StRoot: StSpinPool/StWalgo2011/St2011W_acessMuDst.cxx Source File

00001 // $Id: St2011W_acessMuDst.cxx,v 1.6 2011/02/25 06:03:45 stevens4 Exp $
00002 //
00003 //*-- Author : Jan Balewski, MIT
00004 //*-- Author for Endcap: Justin Stevens, IUCF
00005 
00006 //MuDst
00007 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
00008 #include <StMuDSTMaker/COMMON/StMuDst.h>
00009 #include <StMuDSTMaker/COMMON/StMuTriggerIdCollection.h>
00010 #include <StMuDSTMaker/COMMON/StMuEvent.h>
00011 #include <StMuDSTMaker/COMMON/StMuTrack.h>
00012 #include <StMuDSTMaker/COMMON/StMuPrimaryVertex.h>
00013 
00014 #include "StEmcRawMaker/defines.h"
00015 #include "StEmcUtil/database/StBemcTables.h"
00016 
00017 #include "StEmcUtil/geometry/StEmcGeom.h"
00018 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00019 
00020 #include "St2011WMaker.h"
00021 //--------------------------------------
00022 //--------------------------------------
00023 int  
00024 St2011WMaker::accessBarrelTrig(){ // return non-zero on abort 
00025 
00026   if (isMC){
00027     
00028     /*
00029       When the trigger emulator is ready, this should hook into that
00030       instead of the two functions used below.  For now, check that is passes both
00031       L0 and L2, and set the l2bitET flag to true if so.
00032     */
00033     
00034     //if (!passes_L0()) return -1;
00035     if(!passes_L2()) return -2;
00036     hA[0]->Fill("L2bwET",1.);
00037 
00038     wEve->l2bitET=true;
00039     return 0; // we haven't set everything, but it should be good enough for simu.
00040   }
00041 
00042   StMuEvent* muEve = mMuDstMaker->muDst()->event();
00043 
00044   //collect info for the luminosity monitor
00045   int highestT=0;
00046   int highestM=0;  
00047   for (int m=0;m<300;m++)
00048     {
00049       int myT=muEve->emcTriggerDetector().highTower(m);
00050       if  (myT>highestT)
00051         {
00052           highestT=myT;
00053           highestM=m;
00054         }
00055     }
00056   int highestPhi, tempPhi, tempEta;
00057   int awaySum[16];
00058   int totalSum=0;
00059   for (int i=0;i<16;i++) awaySum[i]=0;
00060   
00061   patchToEtaPhi(highestM,&tempEta,&highestPhi);
00062   
00063   for (int m=0;m<300;m++)
00064     {
00065       int myT=muEve->emcTriggerDetector().highTower(m);
00066       patchToEtaPhi(m,&tempEta,&tempPhi);
00067       for (int away_width=0;away_width<16;away_width++)
00068         if ((highestPhi+30-tempPhi)%30>(15-away_width) && (highestPhi+30-tempPhi)%30<(15+away_width))
00069           {
00070             //printf("==> adding %d to awaySum",myT);
00071             awaySum[away_width]+=myT;
00072           }
00073       totalSum+=myT;
00074     }
00075   for (int i=0;i<16;i++)  wEve->trigAwaySum[i]=awaySum[i];
00076   wEve->trigTotalSum=totalSum;
00077     
00078   StMuTriggerIdCollection *tic=&(muEve->triggerIdCollection());
00079   assert(tic);
00080   const StTriggerId &l1=tic->l1();
00081   vector<unsigned int> idL=l1.triggerIds();
00082 
00083   //printf("nTrig=%d, trigID: ",idL.size());
00084   for(unsigned int i=0;i<idL.size(); i++){
00085     char txt[100];
00086     sprintf(txt,"%d",idL[i]);
00087     //printf("%d, ",idL[i]);
00088     hA[1]->Fill(txt,1.);
00089   }
00090   //printf("\n isTrg=%d trgId=%d\n",tic->nominal().isTrigger(par_l2bwTrgID),par_l2bwTrgID);
00091 
00092   //get bX info
00093   StL0Trigger *trig=&(muEve->l0Trigger());
00094   wEve->bx48=trig->bunchCrossingId();
00095   wEve->bx7=trig->bunchCrossingId7bit(mRunNo);
00096 
00097   //check trigger ID
00098   if(!tic->nominal().isTrigger(par_l2bwTrgID)) return -2;
00099   hA[0]->Fill("L2bwId",1.);
00100 
00101   TArrayI& l2Array = muEve->L2Result();
00102   LOG_DEBUG <<Form("AccessL2Decision() from regular muDst: L2Ar-size=%d",l2Array.GetSize())<<endm;
00103   unsigned int *l2res=(unsigned int *)l2Array.GetArray();
00104   const int BEMCW_off=20; // valid only for 2009 & 2011 run
00105   L2wResult2009 *l2algo= ( L2wResult2009 *) &l2res[BEMCW_off];
00106 
00107 #if 0
00108   if(l2algo->trigger==0) return -3;
00109   printf(" L2-jet online results below:\n");
00110   for (int k=0;k<64;k++) 
00111     if(l2res[k]) printf("k=%2d  val=0x%04x\n",k,l2res[k]);
00112   L2wResult2009_print(l2algo);
00113 #endif
00114 
00115   wEve->l2bitET=(l2algo->trigger&2)>0;  // bit1=ET>thr
00116   wEve->l2bitRnd=(l2algo->trigger&1)>0; // bit0=rnd,
00117 
00118   if( (wEve->l2bitRnd || wEve->l2bitET)==0) return -3; // L2W-algo did not accept this event
00119   hA[0]->Fill("L2bwBits",1.); // confirmation bits were set properly
00120 
00121   if(wEve->l2bitRnd) {
00122     hA[0]->Fill("L2bwRnd",1.);
00123     for (int m=0;m<300;m++){
00124       int val=muEve->emcTriggerDetector().highTower(m);
00125       hA[7]->Fill(val);
00126     }
00127     hA[61]->Fill(wEve->bx7);
00128   }
00129 
00130   if(!wEve->l2bitET)  return -3; // drop L2W-random accepts
00131   if(wEve->l2bitET) hA[0]->Fill("L2bwET",1.);
00132 
00133   //.... only monitor below ....
00134   hA[2]->Fill(wEve->bx48);
00135   hA[3]->Fill(wEve->bx7);
00136 
00137   // access L0-HT data
00138   for (int m=0;m<300;m++)       {
00139     int val=muEve->emcTriggerDetector().highTower(m);
00140     if(wEve->l2bitET) hA[6]->Fill(val);
00141     if(val<par_DsmThres) continue;
00142     if(wEve->l2bitET) hA[8]->Fill(m);
00143     //printf("Fired L0 HT m=%d val=%d\n",m,val);
00144   }
00145   return 0;
00146 }
00147 
00148 
00149 //--------------------------------------
00150 //--------------------------------------
00151 int
00152 St2011WMaker::accessVertex(){ // return non-zero on abort 
00153   int nInpPrimV=mMuDstMaker->muDst()->numberOfPrimaryVertices();
00154   
00155   if(nInpPrimV <par_minPileupVert) return -1;
00156   //separate histos for barrel and endcap triggers
00157   if(wEve->l2bitET) hA[0]->Fill("tpcOn",1.);
00158   if(wEve->l2EbitET) hE[0]->Fill("tpcOn",1.);
00159 
00160   int nVer=0; int nVerR=0;
00161   for(int iv=0;iv<nInpPrimV;iv++) {
00162     StMuPrimaryVertex* V= mMuDstMaker->muDst()->primaryVertex(iv);
00163     assert(V);
00164     mMuDstMaker->muDst()->setVertexIndex(iv);
00165     float rank=V->ranking(), funnyR=999;
00166     if(rank>1e6)  funnyR=log(rank-1e6)+10;
00167     else if(rank>0)   funnyR=log(rank);
00168     else   funnyR=log(rank+1e6)-10;
00169     if(wEve->l2bitET) hA[10]->Fill(funnyR);
00170     if(wEve->l2EbitET) hE[10]->Fill(funnyR);
00171     //keep some neg. rank vertices for endcap if matched to ETOW
00172     if (rank<=0 && V->nEEMCMatch()<=0) continue; 
00173     const StThreeVectorF &r=V->position();
00174     //   StThreeVectorF &er=V->posError();
00175     if(wEve->l2bitET && rank>0) hA[11]->Fill(r.z());
00176     if(wEve->l2EbitET) hE[11]->Fill(r.z());
00177     nVer++; // count valid vertices
00178     if(rank>0) nVerR++; //count vertices with rank>0
00179     if(fabs(r.z()) > par_vertexZ) continue;
00180     WeveVertex wv;
00181     wv.id=iv;
00182     wv.z=r.z();
00183     wv.rank=rank;
00184     wv.funnyRank=funnyR;
00185     wv.nEEMCMatch=V->nEEMCMatch();
00186     wEve->vertex.push_back(wv);
00187   }
00188   if(nVer<=0) return -2;
00189   if(wEve->l2bitET && nVerR>0) {
00190     hA[0]->Fill("primVert",1.);
00191     hA[4]->Fill(wEve->bx48);
00192     hA[5]->Fill(wEve->bx7);
00193   }
00194   if(wEve->l2EbitET) {
00195     hE[0]->Fill("primVert",1.);
00196     hE[4]->Fill(wEve->bx48);
00197     hE[5]->Fill(wEve->bx7);
00198   }
00199 
00200   // access L0-HT data
00201   StMuEvent* muEve = mMuDstMaker->muDst()->event();
00202   for (int m=0;m<300;m++)       {
00203     int val=muEve->emcTriggerDetector().highTower(m);
00204     if(val<par_DsmThres) continue;
00205     if(wEve->l2bitET && nVerR>0) hA[9]->Fill(m);
00206   }
00207   for (int m=0;m<90;m++)        {
00208     int val=muEve->emcTriggerDetector().highTowerEndcap(m);
00209     if(val<parE_DsmThres) continue;
00210     if(wEve->l2EbitET) hE[9]->Fill(m);
00211   }
00212 
00213   if(wEve->l2bitET) hA[12]->Fill(nVerR);
00214   if(wEve->l2EbitET) hE[12]->Fill(nVer);
00215 
00216   if(wEve->vertex.size()<=0) return -3;
00217   if(wEve->l2bitET && nVerR>0) hA[0]->Fill("vertZ",1.);
00218   if(wEve->l2EbitET) hE[0]->Fill("vertZ",1.);
00219   return 0;
00220 }
00221 
00222 
00223 
00224 //--------------------------------------
00225 //--------------------------------------
00226 int  
00227 St2011WMaker::accessTracks(){ // return non-zero on abort 
00228   int nTrOK=0;
00229   // printf("\n nInp=%d eveID=%d nPVer=%d nAnyV= %d\n",nInpEve,mMuDstMaker->muDst()->event()->eventId(),wEve->vertex.size(),mMuDstMaker->muDst()->numberOfPrimaryVertices());
00230   for(uint iv=0;iv<wEve->vertex.size(); iv++) {
00231     uint vertID=wEve->vertex[iv].id;
00232     assert(vertID<mMuDstMaker->muDst()->numberOfPrimaryVertices());
00233     assert(vertID>=0);
00234     StMuPrimaryVertex* V= mMuDstMaker->muDst()->primaryVertex(vertID);
00235     assert(V);
00236     mMuDstMaker->muDst()->setVertexIndex(vertID);
00237     float rank=V->ranking();
00238     assert(rank>0 || (rank<0 && V->nEEMCMatch()));
00239     Int_t nPrimTrAll=mMuDstMaker->muDst()->GetNPrimaryTrack();
00240     for(int itr=0;itr<nPrimTrAll;itr++) {
00241       StMuTrack *prTr=mMuDstMaker->muDst()->primaryTracks(itr);
00242       if(prTr->flag()<=0) continue;
00243       const StMuTrack *glTr=prTr->globalTrack();
00244       if(glTr==0) continue; // see the reason at the end of this method
00245       //keep list of all tracks for TPC cone sum in tree ana
00246       wEve->vertex[iv].prTrList.push_back(prTr); 
00247       StThreeVectorF ro=glTr->lastPoint();
00248 
00249       // TPC+prim vertex tracks and short EEMC tracks
00250       if(prTr->flag()!=301 && prTr->flag()!=311) continue;
00251       if(wEve->l2bitET && rank>0 && prTr->flag()==301) 
00252         hA[20]->Fill("flag",1.); 
00253       if(wEve->l2EbitET && ro.pseudoRapidity()>parE_trackEtaMin)
00254         hE[20]->Fill("flag",1.); 
00255       
00256       float pt=prTr->pt();
00257       if(pt<1.0) continue;
00258       if(wEve->l2bitET && rank>0 && prTr->flag()==301) 
00259         hA[20]->Fill("pt1",1.);
00260       if(wEve->l2EbitET && ro.pseudoRapidity()>parE_trackEtaMin)
00261         hE[20]->Fill("pt1",1.);
00262   
00263       //accepted tracks......
00264       float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
00265       StThreeVectorF ri=glTr->firstPoint();
00266       /* Victor: in reality mChiSqXY is a normal Xi2 for track and mChiSqZ is Xi2 of fit to  primary vertex   */
00267       float globChi2dof=glTr->chi2();
00268       float dedx=prTr->dEdx()*1e6;
00269      
00270       //barrel algo track monitors
00271       if(wEve->l2bitET && rank>0 && prTr->flag()==301){ 
00272         hA[21]->Fill(prTr->nHitsFit());
00273         hA[22]->Fill(hitFrac);
00274         hA[23]->Fill(ri.perp());
00275         hA[24]->Fill(ro.perp());
00276 
00277         //TPC sector dependent filter 
00278         int secID=WtpcFilter::getTpcSec(ro.phi(),ro.pseudoRapidity());
00279         if (mTpcFilter[secID-1].accept(prTr)==false) continue;
00280         //if (secID==20) continue; //poorly calibrated sector?
00281 
00282         hA[25]->Fill(glTr->p().perp());
00283         if(glTr->charge()<0) hA[27]->Fill(glTr->p().perp());
00284         hA[29]->Fill(pt);
00285         if(prTr->charge()<0)hA[30]->Fill(pt);
00286         hA[26]->Fill(ro.pseudoRapidity(),ro.phi());
00287         if(pt>5) //estimate TPC inefficiency in data
00288           hA[57]->Fill(ro.pseudoRapidity(),ro.phi());
00289         hA[35]->Fill(globChi2dof);
00290         // monitor chi2 for east/west TPC separately
00291         if(ri.z()>0 && ro.z()>0)  hA[58]->Fill(globChi2dof);
00292         if(ri.z()<0 && ro.z()<0)  hA[59]->Fill(globChi2dof);
00293         hA[36]->Fill(globChi2dof,ro.pseudoRapidity());
00294         hA[28]->Fill(prTr->p().mag(),dedx);
00295 
00296         if(pt>10) 
00297           hA[197]->Fill(ro.pseudoRapidity(),ro.phi());
00298         hA[198]->Fill(ro.pseudoRapidity(),prTr->pt());
00299       }
00300 
00301       //endcap algo track monitors
00302       if(wEve->l2EbitET && ro.pseudoRapidity()>parE_trackEtaMin){ 
00303         hE[20]->Fill("#eta>0.7",1.);
00304         hE[21]->Fill(prTr->nHitsFit());
00305         hE[22]->Fill(hitFrac);
00306         hE[23]->Fill(ri.perp());
00307         hE[24]->Fill(ro.perp());
00308 
00309         // TPC sector dependent filter 
00310         int secID=WtpcFilter::getTpcSec(ro.phi(),ro.pseudoRapidity());
00311         if ( mTpcFilterE[secID-1].accept(prTr)==false) continue;
00312         //if (secID==20) continue; //poorly calibrated sector?
00313 
00314         hE[25]->Fill(glTr->p().perp());
00315         if(glTr->charge()<0)hE[27]->Fill(glTr->p().perp());
00316         hE[29]->Fill(pt);
00317         if(prTr->charge()<0)hE[30]->Fill(pt);
00318         
00319         hE[26]->Fill(ro.pseudoRapidity(),ro.phi());
00320         if(pt>5) //estimate TPC inefficiency in data
00321           hE[57]->Fill(ro.pseudoRapidity(),ro.phi());
00322         hE[35]->Fill(globChi2dof);
00323         hE[36]->Fill(globChi2dof,ro.pseudoRapidity());
00324         hE[28]->Fill(prTr->p().mag(),dedx);
00325       }
00326 
00327       
00328       bool barrelTrack=(wEve->l2bitET && rank>0 && prTr->flag()==301 && pt>par_trackPt); 
00329       if(barrelTrack) hA[20]->Fill("ptOK",1.);//good barrel candidate
00330       bool endcapTrack=(wEve->l2EbitET && ro.pseudoRapidity()>parE_trackEtaMin && pt>parE_trackPt); 
00331       if(endcapTrack) hE[20]->Fill("ptOK",1.);//good endcap candidate
00332 
00333       if(!barrelTrack && !endcapTrack) continue;
00334       
00335       //keep all tracks in one container
00336       nTrOK++;
00337       WeveEleTrack wTr;
00338       
00339       wTr.prMuTrack=prTr;
00340       wTr.glMuTrack=glTr;
00341       StThreeVectorF prPvect=prTr->p();
00342       wTr.primP=TVector3(prPvect.x(),prPvect.y(),prPvect.z());
00343 
00344       wEve->vertex[iv].eleTrack.push_back(wTr);
00345     }// loop over tracks
00346   }// loop over vertices
00347   if(nTrOK<=0) return -1;
00348   if(wEve->l2bitET) hA[0]->Fill("Pt10",1.);
00349   if(wEve->l2EbitET) hE[0]->Fill("Pt10",1.);
00350   return 0;
00351 }
00352 
00353 /* from Pibero:
00354    It looks like your global track is null. See this post:
00355 
00356    http://www.star.bnl.gov/HyperNews-star/get/mudst/53.html
00357 
00358    My reading of this hypernews says its just the way ITTF/MuDst
00359    works. You can get a good primary track, but its global track
00360    fails the chi2 fit. So the primary track is kept in the MuDst
00361    but the global track is dropped. I would suggest you skip those
00362    rare primary tracks that have no global tracks, that way you
00363    still use most of the tracks in the MuDst. You don't need to
00364    skip the entire event, just that track. I guess the down side
00365    is you couldn't make a global DCA cut on those rare tracks, right?
00366    I guess you could also request S&C to change ITTF/MuDst not to drop
00367    the global track for every good primary track regardless of chi2.
00368 */
00369 
00370   
00371 
00372 /* $STAR/StRoot/StEvent/StTrack.h
00373  *  mFlag=zxyy, where  z = 1 for pile up track in TPC (otherwise 0) 
00374  *                     x indicates the detectors included in the fit and 
00375  *                    yy indicates the status of the fit. 
00376  *  Positive mFlag values are good fits, negative values are bad fits. 
00377  *
00378  *  The first digit indicates which detectors were used in the refit: 
00379  *
00380  *      x=1 -> TPC only 
00381  *      x=3 -> TPC       + primary vertex 
00382  *      x=5 -> SVT + TPC 
00383  *      x=6 -> SVT + TPC + primary vertex 
00384  *      x=7 -> FTPC only 
00385  *      x=8 -> FTPC      + primary 
00386  *      x=9 -> TPC beam background tracks            
00387  *
00388  *  The last two digits indicate the status of the refit: 
00389  *       = +x01 -> good track 
00390  *
00391  *      = -x01 -> Bad fit, outlier removal eliminated too many points 
00392  *      = -x02 -> Bad fit, not enough points to fit 
00393  *      = -x03 -> Bad fit, too many fit iterations 
00394  *      = -x04 -> Bad Fit, too many outlier removal iterations 
00395  *      = -x06 -> Bad fit, outlier could not be identified 
00396  *      = -x10 -> Bad fit, not enough points to start 
00397  *
00398  *      = +x11 -> Short track pointing to EEMC
00399 */
00400 
00401 
00402 //________________________________________________
00403 //________________________________________________
00404 int
00405 St2011WMaker::accessBTOW(){
00406   
00407   StMuEmcCollection* emc = mMuDstMaker->muDst()->muEmcCollection();
00408   if (!emc) {
00409     gMessMgr->Warning() <<"No EMC data for this event"<<endm;    return -4;
00410   }
00411     
00412   int ibp=kBTow; // my index for tower & preshower set to BTOW
00413   int jBP=BTOW; // official BTOW detector ID
00414   int n5=0,n0=0,n1=0,n2=0,n3=0,n4=0;
00415   int maxID=0;
00416   double maxADC=0,adcSum=0;
00417   for (int softID=1; softID <=mxBtow ; softID++) {
00418     float rawAdc= emc->getTowerADC(softID);
00419     if(rawAdc==0) n0++;
00420 
00421     int statPed,statOfl,statGain;
00422     mBarrelTables->getStatus(jBP, softID, statPed,"pedestal");
00423     mBarrelTables->getStatus(jBP, softID, statOfl);
00424     mBarrelTables->getStatus(jBP, softID, statGain,"calib");
00425     
00426     if(statPed!=1) {
00427       wEve->bemc.statTile[ibp][softID-1]=1;
00428       n1++; continue;}
00429     if(statOfl!=1) {
00430       wEve->bemc.statTile[ibp][softID-1]=2;
00431       n2++; continue;} 
00432     if(statGain!=1) {
00433       wEve->bemc.statTile[ibp][softID-1]=4;
00434       n3++; continue;} 
00435 
00436     wEve->bemc.statTile[ibp][softID-1]=0 ; 
00437 
00438 
00439     float ped,sigPed,gain;
00440     int capID=0;// just one value for btow
00441     mBarrelTables->getPedestal(jBP,softID,capID,ped,sigPed); 
00442     mBarrelTables->getCalib(jBP, softID, 1, gain);
00443     if (use_gains_file == 1) {
00444       gain = gains_BTOW[softID];
00445     }
00446     //printf("id=%d gain=%f\n",softID,gain);
00447 
00448     //method for shifting energy scale 
00449     gain=gain*par_btowScale;//(default is par_btowScale=1)
00450 
00451     float adc=rawAdc-ped;
00452     if(adc>0) n4++;
00453     if(adc<par_kSigPed*sigPed) continue;
00454     if(adc<par_AdcThres) continue;
00455     n5++;
00456     wEve->bemc.adcTile[ibp][softID-1]=adc;
00457     wEve->bemc.eneTile[ibp][softID-1]=adc*gain;
00458 
00459     if(maxADC<adc) { maxID=softID; maxADC=adc;}
00460     adcSum+=adc;
00461   }
00462 
00463   if(isMC>0 && isMC<20) assert(n1==0);//prevent using real peds for private MC
00464   //printf("NNN %d %d %d %d %d %d id=%d\n",n0,n1,n2,n3,n4,n5,maxID);
00465   if(n0==mxBtow) return -1 ;  // BTOW was not present in this events
00466 
00467   wEve->bemc.tileIn[ibp]=1; //tag usable data
00468   
00469   if(nInpEve%5000==1) { 
00470     LOG_INFO << Form("unpackMuBTOW() dataIn=%d, nBbad: ped=%d stat=%d gain=%d ; nAdc: %d>0, %d>thres\n    maxADC=%.0f softID=%d adcSum=%.0f",
00471                      wEve->bemc.tileIn[ibp],n1,n2,n3,n4,n5,
00472                      maxADC,maxID,adcSum
00473                      )<<endm;
00474   }
00475   hA[31]->Fill(maxADC);
00476   hA[32]->Fill(adcSum);
00477   wEve->bemc.maxAdc=maxADC;
00478 
00479   if(maxID<=2400) hA[195]->Fill(maxADC);
00480   else hA[196]->Fill(maxADC);
00481 
00482   if(maxADC<par_maxADC)  return -2 ;  // not enough energy
00483   
00484   return 0;
00485 }
00486 
00487 
00488 //________________________________________________
00489 //________________________________________________
00490 void
00491 St2011WMaker::fillTowHit(bool vert){
00492   if(!wEve->l2bitET) return; //only barrel triggers
00493 
00494   //find highest rank vertex
00495   float maxRank=0; uint maxRankId=0;
00496   for(uint iv=0;iv<wEve->vertex.size(); iv++) {
00497     float rank=wEve->vertex[iv].rank;
00498     if(rank<0) continue;
00499     if(rank > maxRank){
00500       maxRank=rank;
00501       maxRankId=iv;
00502     }
00503   }
00504 
00505   int bx7=wEve->bx7; int bxBin=-1;
00506   if(bx7>=0 && bx7<30)
00507     bxBin=0;
00508   else if(bx7<40)
00509     bxBin=1;
00510   else if(bx7<110)
00511     bxBin=2;
00512   else if(bx7<120)
00513     bxBin=3;
00514   
00515   float Rcylinder= mBtowGeom->Radius(), Rcylinder2=Rcylinder*Rcylinder;
00516   //loop barrel towers and fill histo
00517   for(int i=0; i<mxBtow; i++){
00518     float adc=wEve->bemc.adcTile[kBTow][i];
00519     bool fillAdc=false;
00520     if(adc > 10) fillAdc=true; //~150 MeV threshold for tower firing
00521     if(vert){
00522       if(fillAdc) hA[215+bxBin]->Fill(positionBtow[i].Eta(),positionBtow[i].Phi()); 
00523       float ene=wEve->bemc.eneTile[kBTow][i];
00524       float delZ=positionBtow[i].z()-wEve->vertex[maxRankId].z;
00525       float e2et=Rcylinder/sqrt(Rcylinder2+delZ*delZ);
00526       float ET=ene*e2et;
00527       if(ET > 2.0) hA[219+bxBin]->Fill(positionBtow[i].Eta(),positionBtow[i].Phi());
00528     }
00529     else if(fillAdc) hA[223+bxBin]->Fill(positionBtow[i].Eta(),positionBtow[i].Phi());
00530   }
00531   
00532   for(int isec=0;isec<mxEtowSec;isec++){
00533     for(int isub=0;isub<mxEtowSub;isub++){
00534       for(int ieta=0;ieta<mxEtowEta;ieta++){
00535         int iPhi=isec*mxEtowSub+isub;
00536         float adc=wEve->etow.adc[iPhi][ieta];
00537         bool fillAdc=false;
00538         if(adc > 10) fillAdc=true; //~150 MeV threshold for tower firing
00539         
00540         if(vert){
00541           if(fillAdc) hA[227+bxBin]->Fill(ieta,iPhi); 
00542           float ene=wEve->etow.ene[iPhi][ieta];
00543           float delZ=positionEtow[iPhi][ieta].z()-wEve->vertex[maxRankId].z;
00544           float Rxy=positionEtow[iPhi][ieta].Perp();
00545           float e2et=Rxy/sqrt(Rxy*Rxy+delZ*delZ);
00546           float ET=ene*e2et;
00547           if(ET > 2.0) hA[231+bxBin]->Fill(ieta,iPhi);
00548         }
00549         else if(fillAdc) hA[235+bxBin]->Fill(ieta,iPhi);
00550       }
00551     }
00552   }
00553   
00554   return;
00555 }
00556 
00557 
00558 //________________________________________________
00559 //________________________________________________
00560 float
00561 St2011WMaker::sumTpcCone(int vertID, TVector3 refAxis, int flag, int pointTowId){ 
00562 
00563   // flag=2 use 2D cut, 1= only delta phi
00564 
00565   // printf("******* sumTpcCone, flag=%d eveId=%d vertID=%d  eta0=%.2f phi0/rad=%.2f  \n",flag,wEve->id,vertID,refAxis.PseudoRapidity() ,refAxis.Phi());
00566 
00567   assert(vertID>=0);
00568   assert(vertID<(int)mMuDstMaker->muDst()->numberOfPrimaryVertices());
00569   
00570   StMuPrimaryVertex* V= mMuDstMaker->muDst()->primaryVertex(vertID);
00571   assert(V);
00572   mMuDstMaker->muDst()->setVertexIndex(vertID);
00573   float rank=V->ranking();
00574   assert(rank>0 || (rank<0 && V->nEEMCMatch()));
00575   double ptSum=0;
00576   Int_t nPrimTrAll=mMuDstMaker->muDst()->GetNPrimaryTrack();
00577   for(int itr=0;itr<nPrimTrAll;itr++) {
00578     StMuTrack *prTr=mMuDstMaker->muDst()->primaryTracks(itr);
00579     if(prTr->flag()<=0) continue;
00580     if(prTr->flag()!=301 && pointTowId>0) continue;// TPC-only regular tracks for barrel candidate
00581     if(prTr->flag()!=301 && prTr->flag()!=311 && pointTowId<0) continue;// TPC regular and short EEMC tracks for endcap candidate
00582     float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
00583     if(hitFrac<par_nHitFrac) continue;
00584     StThreeVectorF prPvect=prTr->p();
00585     TVector3 primP=TVector3(prPvect.x(),prPvect.y(),prPvect.z());
00586     // printf(" prTrID=%4d  prTrEta=%.3f prTrPhi/deg=%.1f prPT=%.1f  nFitPts=%d\n", prTr->id(),prTr->eta(),prTr->phi()/3.1416*180.,prTr->pt(),prTr->nHitsFit());
00587     if(flag==1) {
00588       float deltaPhi=refAxis.DeltaPhi(primP);
00589       if(fabs(deltaPhi)> par_awayDeltaPhi) continue;
00590     }
00591     if(flag==2) {
00592       float deltaR=refAxis.DeltaR(primP);
00593       if(deltaR>par_nearDeltaR) continue;
00594     }
00595     float pT=prTr->pt();
00596     //    printf(" passed pt=%.1f\n",pT);
00597 
00598     //separate quench for barrel and endcap candidates
00599     if(pT>par_trackPt && pointTowId>0) ptSum+=par_trackPt;
00600     else if(pT>parE_trackPt && pointTowId<0) ptSum+=parE_trackPt;
00601     else  ptSum+=pT;
00602   }
00603   return ptSum;
00604 }
00605 
00606 
00607 
00608 //________________________________________________
00609 void
00610 St2011WMaker::accessBSMD(){
00611   const char cPlane[ mxBSmd]={'E','P'};
00612   // Access to muDst .......................
00613   StMuEmcCollection* emc = mMuDstMaker->muDst()->muEmcCollection();
00614   if (!emc) {
00615     gMessMgr->Warning() <<"No EMC data for this muDst event"<<endm;    return;
00616   }
00617 
00618   //....................... B S M D .........................
00619   for(int iEP=bsmde; iEP<=bsmdp;iEP++) { // official BSMD plane IDs
00620     int iep=iEP-3; assert(bsmde==3);// what a hack
00621     int nh= emc->getNSmdHits(iEP);
00622     //printf("muDst BSMD-%c nHit=%d\n",cPlane[iep],nh);
00623     int n5=0,n1=0,n2=0,n3=0,n4=0;
00624     for (int i=0; i < nh; i++) {
00625       StMuEmcHit *hit=emc->getSmdHit(i,iEP);
00626       float  adc=hit->getAdc();
00627       int softID=hit->getId();
00628 
00629       int statPed,statOfl,statGain;
00630       mBarrelTables->getStatus(iEP, softID, statPed,"pedestal");
00631       mBarrelTables->getStatus(iEP, softID, statOfl);
00632       mBarrelTables->getStatus(iEP, softID, statGain,"calib");
00633       
00634       if(statPed!=1) {
00635         wEve->bemc.statBsmd[iep][softID-1]=1;
00636         n1++; continue;}
00637       if(statOfl!=1) {
00638         wEve->bemc.statBsmd[iep][softID-1]=2;
00639         n2++; continue;} 
00640       if(statGain<1 || statGain>19) {
00641         wEve->bemc.statBsmd[iep][softID-1]=4;
00642         n3++; continue;} 
00643       
00644       float pedRes,sigPed,gain;
00645       int capID=0;// just one value for ped residua in pp500, 2009 run
00646       mBarrelTables->getPedestal(iEP,softID,capID,pedRes,sigPed); 
00647       mBarrelTables->getCalib(iEP, softID, 1, gain);
00648       
00649       if(isMC) { // overwrite it based on genat DE & private calibration
00650         float par_bsmdAbsGain=6e6;// tmp arbitrary absolute calib of bsmd, was 3e6
00651         float  de = hit->getEnergy();// Geant energy deposit (GeV)
00652         adc=de*par_bsmdAbsGain;
00653       } else { // correct for pedestal residua
00654         adc-=pedRes;
00655         if(adc>0) n4++;
00656         if(adc<par_kSigPed*sigPed) continue;
00657       }
00658       
00659       n5++;
00660       assert(softID>=1);      assert(softID<=mxBStrips);
00661       int id0=softID-1;
00662       wEve->bemc.adcBsmd[ iep][id0]=adc;
00663       hA[70+10*iep]->Fill(adc);
00664       
00665       //if(nInpEve<3 || i <20 )printf("  i=%d, smd%c id=%d, m=%d adc=%.3f pedRes=%.1f, sigP=%.1f stat: O=%d P=%d G=%d  gain=%.2f\n",i,cPlane[iep],softID,1+id0/150,adc,pedRes,sigPed, statOfl,statPed,statGain, gain);
00666     }// end of hit list
00667     if(nTrigEve%5000==1) { 
00668       LOG_INFO << Form("unpackMuBSMD-%c() nBbad: ped=%d stat=%d gain=%d ; nAdc: %d>0, %d>thres", cPlane[iep],n1,n2,n3,n4,n5)<<endm;
00669     }
00670   }// end of E-, P-plane loop
00671 
00672 }
00673 
00674 
00675 //$Log: St2011W_acessMuDst.cxx,v $
00676 //Revision 1.6  2011/02/25 06:03:45  stevens4
00677 //addes some histos and enabled running on MC
00678 //
00679 //Revision 1.5  2011/02/17 04:16:20  stevens4
00680 //move sector dependent track QA cuts before track pt>10 cut and lower par_clustET and par_ptBalance thresholds to 14 GeV
00681 //
00682 //Revision 1.4  2011/02/16 15:25:16  stevens4
00683 //remove relative gain correction for BSMD adc
00684 //
00685 //Revision 1.3  2011/02/15 17:39:12  stevens4
00686 //remove accept-all for L2btowW bits
00687 //
00688 //Revision 1.2  2011/02/14 01:36:17  stevens4
00689 //*** empty log message ***
00690 //
00691 //Revision 1.1  2011/02/10 20:33:22  balewski
00692 //start
00693 //