00001 class StChain;
00002 class StMuEmcCollection;
00003 class StMuTriggerIdCollection;
00004 StChain *chain=0;
00005
00006
00007
00008 int rdMuDst2print(
00009 char* file = "st_physics_12033050_raw_4010001.MuDst.root",
00010 int nEve=10,
00011 char* inDir = "./"
00012 )
00013 {
00014 Int_t nFiles = 1;
00015
00016 inDir="/star/data05/scratch/balewski/stW-2012A/data/";
00017 file="st_W_13078014_raw_4360001.MuDst.root";
00018
00019
00020 gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
00021 loadSharedLibraries();
00022
00023 cout << " loading done " << endl;
00024
00025
00026
00027 chain = new StChain("StChain");
00028
00029
00030 muMk = new StMuDstMaker(0,0,inDir,file,"MuDst.root",nFiles);
00031 TChain* tree=muMk->chain(); assert(tree);
00032 int nEntries=(int) tree->GetEntries();
00033 printf("total eve in muDst chain =%d\n",nEntries);
00034 if(nEntries<0) return;
00035
00036
00037 chain->Init();
00038 chain->ls(3);
00039
00040 h1=new TH1F("nP","# of prim tracks per vertex",100,0,200);
00041 h2=new TH1F("vR","PPV Vertex rank; rank", 300, -1.2e6, 1.2e6); h2->SetLineColor(kRed);
00042 h3=new TH1F("vRL","PPV Vertex rank, funny X-axis; X=Log10(rank-1e6)+ offset", 150, -9,21);
00043 h2->GetXaxis()->SetTitleSize(.043);
00044 h2->GetXaxis()->SetTitleSize(.043);
00045
00046 h4=new TH1F("trPhi"," prim tracks phi if PT>1.0 GeV/c; phi (rad)",50,-3.2,3.2);
00047 h5=new TH1F("zVerTrg"," Z trg-vertex , rank>0; Z (cm)", 50,-200,200);
00048 h6=new TH1F("zVerPlp"," Z pileup-vertex , rank<0; Z (cm)", 50,-200,200);
00049 h7=new TH1F("nPrV","# of prim tr used by PPV; rank>0",15,0,30);
00050
00051 int eventCounter=0;
00052 int t1=time(0);
00053
00054 for (Int_t iev=0;iev<nEntries; iev++) {
00055 if(eventCounter>=nEve) break;
00056 chain->Clear();
00057 int stat = chain->Make();
00058 if(stat) break;
00059 eventCounter++;
00060
00061
00062 StMuEvent* muEve = muMk->muDst()->event();
00063
00064 StEventInfo &info=muEve->eventInfo();
00065 int nPrimV=muMk->muDst()->numberOfPrimaryVertices();
00066
00067 StMuTriggerIdCollection *tic=&(muEve->triggerIdCollection());
00068 int trigID=380305;
00069 bool isL2EW=tic->nominal().isTrigger(trigID);
00070
00071 Int_t nGlobTrAll=muMk->muDst()->GetNGlobalTrack();
00072
00073 TArrayI& l2Array = muEve->L2Result();
00074
00075 if(eventCounter%1==0) {
00076 printf("ieve=%d eventID %d nPrimV=%d nGlobTrAll=%d =============\n", eventCounter,info.id(),nPrimV,nGlobTrAll);
00077 printf("TrigID=%d fired=%d\n",trigID,isL2EW);
00078 }
00079 printL0Trig(tic);
00080 printL2Trig(l2Array);
00081
00082 int iv;
00083 if(0)for(iv=0;iv<nPrimV;iv++) {
00084 StMuPrimaryVertex* V= muMk->muDst()->primaryVertex(iv);
00085 assert(V);
00086 muMk->muDst()->setVertexIndex(iv);
00087 StThreeVectorF &r=V->position();
00088 StThreeVectorF &er=V->posError();
00089 printf("iv=%d Vz=%.2f +/-%.2f \n",iv,r.z(),er.z() );
00090
00091 int nPrimTr =0;
00092 int itr;
00093 Int_t nPrimTrAll=muMk->muDst()->GetNPrimaryTrack();
00094 for(itr=0;itr<nPrimTrAll;itr++) {
00095 StMuTrack *pr_track=muMk->muDst()->primaryTracks(itr);
00096 if(pr_track->flag()<=0) continue;
00097 nPrimTr ++;
00098 }
00099 if(nPrimV>0)h1->Fill(nPrimTr);
00100 float rank=V->ranking();
00101 h2->Fill(rank);
00102 if(rank>1e6) h3->Fill(log(rank-1e6)+10);
00103 else if(rank>0) h3->Fill(log(rank));
00104 else h3->Fill(log(rank+1e6)-10);
00105
00106 if(1)printf(" nPrimTr=%d , Z=%.1f VFid=%d:: ntrVF=%d nCtb=%d nBemc=%d nEEmc=%d nTpc=%d sumPt=%.1f rank=%g\n"
00107 ,nPrimTr,r.z(), V->vertexFinderId() ,V->nTracksUsed() ,V->nCTBMatch() ,V-> nBEMCMatch() ,V->nEEMCMatch() ,V->nCrossCentralMembrane() ,V->sumTrackPt() ,V->ranking());
00108 if (rank>0) h5->Fill(r.z());
00109 if (rank>0) h7->Fill(V->nTracksUsed());
00110 if (rank<0) h6->Fill(r.z());
00111 }
00112
00113 if(0)
00114 for(iv=0;iv<nPrimV;iv++) {
00115 muMk->muDst()->setVertexIndex(iv);
00116 Int_t nPrimTrAll=muMk->muDst()->GetNPrimaryTrack();
00117 cout<<"\n Prim "<<nPrimTrAll<<" tracks belonging to "<<iv<<"-th prim vertex"<<endl;
00118
00119 int itr;
00120 int ntr=0;
00121
00122 for(itr=0;itr<nPrimTrAll;itr++) {
00123 StMuTrack *pr_track=muMk->muDst()->primaryTracks(itr);
00124 if(pr_track->flag()<=0) continue;
00125 ntr++;
00126 cout << "\nPrimary track " << ntr << " momentum(P vect) " << pr_track->p() << " PT="<<pr_track->pt()<< " recoCharge="<<pr_track->charge()<< endl; cout << "\t flag=" << pr_track->flag() << " nHits=" << pr_track->nHits()<< " vertID="<< pr_track->vertexIndex()<< endl;
00127 cout << "\t primV("<<iv<<") primDCA=" << pr_track->dca(iv) << endl;
00128 if(pr_track->dca(iv).mag()>5) cout << "^^^^^ 3D DCA magnitude="<<pr_track->dca(iv).mag()<<endl;
00129
00130
00131 if(pr_track->pt()>1.0) h4->Fill(pr_track->phi());
00132
00133 }
00134
00135 }
00136
00137 continue;
00138
00139 StMuEmcCollection* emc = muMk->muDst()->muEmcCollection();
00140 if (!emc) {
00141 printf(" No EMC data for this event\n");
00142 return kStOK;
00143 }
00144
00145
00146
00147 printBEtower(emc);
00148
00149
00150
00151
00152 }
00153 printf("****************************************** \n");
00154
00155 int t2=time(0);
00156 if(t2==t1) t2=t1+1;
00157 float tMnt=(t2-t1)/60.;
00158 float rate=1.*eventCounter/(t2-t1);
00159 printf("sorting done %d of nEve=%d, CPU rate=%.1f Hz, total time %.1f minute(s) \n\n",eventCounter,nEntries,rate,tMnt);
00160 c=new TCanvas(); c->Divide(2,2);
00161 c->cd(1); h5->Fit("gaus");
00162 c->cd(3); h6->Fit("gaus");
00163 c->cd(2); h3->Draw();
00164 c->cd(4); h7->Draw();
00165
00166 return;
00167
00168 }
00169
00170
00171
00172
00173 printEEtower( StMuEmcCollection* emc ) {
00174 int sec,eta,sub,adc;
00175 StMuEmcHit *hit;
00176
00177 int i, nh;
00178
00179 printf("\Total %d hits in Tower (only ADC>0)\n",emc->getNEndcapTowerADC());
00180 nh=0;
00181 for (i=0; i< emc->getNEndcapTowerADC(); i++) {
00182 emc->getEndcapTowerADC(i,adc,sec,sub,eta);
00183 if (adc<=0) continue;
00184 nh++;
00185 printf("i=%d Tower %2.2dT%c%2.2d adc=%4d\n",i,sec,sub+'A'-1,eta,adc );
00186
00187 int adcX=1000+ (eta-1) + (sub-1)*12 +(sec-1)*60;
00188
00189 }
00190 printf(" Total %d towers with ADC>0\n",nh);
00191 }
00192
00193
00194
00195
00196 printEEpre( StMuEmcCollection* emc ) {
00197 int sec,eta,sub,pre,adc;
00198 StMuEmcHit *hit;
00199
00200 int i, nh;
00201 nh= emc->getNEndcapPrsHits();
00202 printf("\nTotal %d hits in pre1+2+post\n",nh);
00203 for (i=0; i<nh; i++) {
00204 hit=emc->getEndcapPrsHit(i,sec,sub,eta,pre);
00205 int ss=sub + 5*(pre-1);
00206 adc=hit->getAdc();
00207 printf("i=%d pre/post(%d) %2.2d%c%c%2.2d : energy=%f adc=%d\n",i,pre,sec,pre+'P'-1,sub+'A'-1,eta,hit->getEnergy(),adc);
00208 int adcX= (eta-1) + (sub-1) *12 +(sec-1)*60 + 1000*pre;
00209
00210
00211
00212 }
00213 }
00214
00215
00216
00217
00218 printEEsmd( StMuEmcCollection* emc ) {
00219 int sec,strip,adc;
00220 char uv='U';
00221
00222 for(uv='U'; uv<='V'; uv++) {
00223 int nh= emc->getNEndcapSmdHits(uv);
00224 printf("\nTotal %d hits in SMD-%c\n",nh,uv);
00225 for (int i=0; i<nh; i++) {
00226 hit=emc->getEndcapSmdHit(uv,i,sec,strip);
00227 adc=hit->getAdc();
00228 printf(" SMD-%c %2.2d%c%3.3d : energy=%f adc=%d\n",uv,sec,uv,strip,hit->getEnergy(),adc);
00229 int adcX= 1000 + strip-1 +(sec-1)*300;
00230
00231 }
00232 }
00233 }
00234
00235
00236
00237 printBEtower( StMuEmcCollection* emc ) {
00238 int sec,eta,sub,adc;
00239 StMuEmcHit *hit;
00240
00241 int i, nh;
00242
00243 printf("\Total hits in Tower (only ADC>0)\n");
00244 nh=0;
00245 for (i=0; i< 4800; i++) {
00246 int adc = emc->getTowerADC(i);
00247 if (adc<=1000) continue;
00248 nh++;
00249 printf(" Tower id=%d adc=%4d\n",i,adc );
00250
00251 }
00252 printf(" --> %d towers with ADC>thr\n",nh);
00253 }
00254
00255
00256
00257
00258
00259
00260 printBEpre( StMuEmcCollection* emc ) {
00261
00262 int i, nh;
00263 nh = emc->getNPrsHits();
00264 printf("\nTotal %d hits in pre1\n",nh);
00265
00266 int n1=0;
00267 for (i=0; i<nh; i++) {
00268 StMuEmcHit * hit = emc->getPrsHit(i);
00269 int adc = hit->getAdc();
00270 int id=hit->getId();
00271 if(adc<4) continue;
00272 n1++;
00273 printf("BPRS i=%d, id=%d adc = %d\n",i,id,adc);
00274 }
00275 printf(" --> %d BPRS hits with ADC>thr\n",n1);
00276 }
00277
00278
00279
00280 printBEsmd( StMuEmcCollection* emc ) {
00281
00282 int n1=0,n2=0;
00283 int nh = emc->getNSmdHits(3);
00284 printf("\nTotal %d hits in SMDE\n",nh);
00285 for (int i=0; i<nh; i++) {
00286 StMuEmcHit * hit = emc->getSmdHit(i,3);
00287 adc = hit->getAdc();
00288 int id=hit->getId();
00289 if(adc<4) continue;
00290 n1++;
00291 printf("BSMDE i=%d, id=%d adc = %d\n",i,id,adc);
00292 }
00293
00294 int nh = emc->getNSmdHits(4);
00295 printf("\nTotal %d hits in SMDP\n",nh);
00296 for (int i=0; i<nh; i++) {
00297 hit = emc->getSmdHit(i,4);
00298 adc = hit->getAdc();
00299 int id=hit->getId();
00300 if(adc<4) continue;
00301 n2++;
00302 printf("BSMDP i=%d,id=%d adc = %d\n",i,id,adc);
00303 }
00304
00305 printf(" --> %d BSMD-E & %d BSMD-P hits with ADC>thr\n",n1,n2);
00306 }
00307
00308
00309 printL2Trig( TArrayI& l2Array) {
00310
00311 printf("AccessL2Decision() from regular muDst: L2Ar-size=%d\n",l2Array.GetSize());
00312 unsigned int *l2res=(unsigned int *)l2Array.GetArray();
00313 for(int i=0;i<l2Array.GetSize();i++){
00314 if(l2res[i]==0) continue;
00315 printf("i=%d val=0x%08x\n",i,l2res[i]);
00316 }
00317 }
00318
00319
00320
00321 void printL0Trig(StMuTriggerIdCollection *tic){
00322
00323 const StTriggerId &l1=tic->l1();
00324 vector<unsigned int> idL=l1.triggerIds();
00325 printf("nTrig=%d, trigID: ",idL.size());
00326 for(unsigned int i=0;i<idL.size(); i++){
00327 printf("%d, ",idL[i]);
00328 }
00329 printf("\n dump L2Array:");
00330
00331
00332
00333
00334
00335
00336 }
00337