00001 class StChain;
00002 class StMuEmcCollection;
00003 StChain *chain=0;
00004 const int mxv=4;
00005 TH1F *hA[8], *hT[4], *hF[4];
00006
00007 TFile *hf=0;
00008
00009 int rdMu2PrimTr(
00010 char* file = "st_physics_6145042_raw_2040013.MuDst.root",
00011 Int_t nFiles = 2,
00012 char* inDir = "out1/",
00013 int nEve=100) {
00014
00015 initHisto();
00016
00017
00018
00019 inDir="/star/data04/sim/balewski/tmp2/";
00020 file="st_physics_6111024_raw_1040001.MuDst.root";
00021
00022 gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
00023 loadSharedLibraries();
00024 cout << " loading done " << endl;
00025
00026
00027
00028 chain = new StChain("StChain");
00029
00030
00031 muMk = new StMuDstMaker(0,0,inDir,file,"MuDst.root",nFiles);
00032 TChain* tree=muMk->chain(); assert(tree);
00033 int nEntries=(int)tree->GetEntries();
00034 printf("total eve in chain =%d\n",nEntries);
00035 printf("in=%s%s=\n",inDir,file);
00036
00037 chain->Init();
00038 chain->ls(3);
00039
00040 int eventCounter=0;
00041 int stat=0;
00042
00043
00044 while ( stat==0 ) {
00045 if(eventCounter>=nEve) break;
00046 chain->Clear();
00047 stat = chain->Make();
00048 if(stat) break;
00049
00050 eventCounter++;
00051
00052
00053 StMuEvent* muEve = muMk->muDst()->event();
00054
00055 StEventInfo &info=muEve->eventInfo();
00056 int nPrimV=muMk->muDst()->numberOfPrimaryVertices();
00057 StMuTriggerIdCollection &tic=muEve->triggerIdCollection();
00058
00059 int trigID=96211;
00060 bool fired=tic.nominal().isTrigger(trigID);
00061
00062 Int_t nPrimTrAll=muMk->muDst()->GetNPrimaryTrack();
00063 Int_t nGlobTrAll=muMk->muDst()->GetNGlobalTrack();
00064
00065 if(eventCounter%1000==0) printf("\n\n ====================%d processing eventID %d nPrimV=%d nPrimTr=%d nGlobTrAll=%d =============\n", eventCounter,info.id(),nPrimV,nPrimTrAll,nGlobTrAll);
00066
00067
00068
00069
00070 hA[0]->Fill(nPrimV);
00071 int iv;
00072 for(iv=0;iv<nPrimV;iv++) {
00073 StMuPrimaryVertex* V= muMk->muDst()->primaryVertex(iv);
00074 assert(V);
00075 StThreeVectorF &r=V->position();
00076 StThreeVectorF &er=V->posError();
00077
00078
00079
00080 float rank1K=V->ranking()/1000.;
00081 hA[5]->Fill(er.z() );
00082 ((TH2F*) hA[6])->Fill(er.z() ,rank1K);
00083
00084 if(nPrimV==1 && iv==0) hA[1]->Fill(rank1K);
00085 if(nPrimV>1){
00086 if(iv==0)hA[2]->Fill(rank1K);
00087 if(iv==1)hA[3]->Fill(rank1K);
00088 if(iv==2)hA[4]->Fill(rank1K);
00089 if(iv==1)hA[7]->Fill(rank1K/muMk->muDst()->primaryVertex(0)->ranking()*1000.);
00090
00091 }
00092 int nPrimTr =0;
00093 int itr;
00094 for(itr=0;itr<nPrimTrAll;itr++) {
00095 StMuTrack *pr_track=muMk->muDst()->primaryTracks(itr);
00096 if(pr_track->vertexIndex()!=iv) continue;
00097 if(pr_track->flag()<=0) continue;
00098 if(iv<mxv) {
00099 if(pr_track->topologyMap().trackTpcOnly()) hT[iv]->Fill(pr_track->nHitsFit());
00100 if(pr_track->topologyMap().trackFtpc()) hF[iv]->Fill(pr_track->nHitsFit());
00101 }
00102 nPrimTr ++;
00103 }
00104 continue;
00105 printf(" nPrimTr=%d , VFid=%d:: ntrVF=%d nCtb=%d nBemc=%d nEEmc=%d nTpc=%d sumPt=%.1f rank=%g\n"
00106 ,nPrimTr, V->vertexFinderId() ,V->nTracksUsed() ,V->nCTBMatch() ,V-> nBEMCMatch() ,V->nEEMCMatch() ,V->nCrossCentralMembrane() ,V->sumTrackPt() ,V->ranking());
00107 }
00108
00109 continue;
00110
00111 for(iv=0;iv<nPrimV;iv++) {
00112 printf(" Prim tracks belonging to %d prim vertex:\n",iv);
00113 int itr;
00114 int ntr=0;
00115 for(itr=0;itr<nPrimTrAll;itr++) {
00116 StMuTrack *pr_track=muMk->muDst()->primaryTracks(itr);
00117 if(pr_track->vertexIndex()!=iv) continue;
00118 if(pr_track->flag()<=0) continue;
00119 ntr++;
00120 cout << "\nPrimary track " << ntr << " momentum " << pr_track->p() << endl; cout << "\t flag=" << pr_track->flag() << " nHits=" << pr_track->nHits()<< " vertID="<< pr_track->vertexIndex()<< endl;
00121 cout << "\t primV("<<iv<<") primDCA=" << pr_track->dca(iv) << endl;
00122 if(pr_track->dca(iv).mag()>5) cout << "^^^^^ 3D DCA magnitude="<<pr_track->dca(iv).mag()<<endl;
00123
00124
00125 }
00126 }
00127
00128 continue;
00129 }
00130 hf->Write();
00131
00132 }
00133
00134
00135
00136 void initHisto() {
00137 hf=new TFile("out.hist.root","recreate");
00138 hA[0]=new TH1F("nV","# of prim vertices per eve; # vertex",10,-0.5,9.5);
00139 int iv;
00140 for(iv=0;iv<mxv;iv++) {
00141 char tt1[100], tt2[300];
00142 sprintf(tt1,"nPTv%d",iv);
00143 sprintf(tt2,"nFitPoint for TPC prim tracks for vertex%d",iv);
00144 hT[iv]=new TH1F(tt1,tt2,50,-0.5,49.5);
00145
00146 sprintf(tt1,"nPFv%d",iv);
00147 sprintf(tt2,"nFitPoint for FTPC prim tracks for vertex%d",iv);
00148 hF[iv]=new TH1F(tt1,tt2,50,-0.5,49.5);
00149 }
00150
00151 hA[1]=new TH1F("L11","vertex rank of 1st vertex if only one found; vertex rank/1000",100,0,20);
00152 hA[2]=new TH1F("Ln1","vertex rank of 1st vertex if more than one found; vertex rank/1000",100,0,20);
00153 hA[3]=new TH1F("Ln2","vertex rank of 2nd vertex if more than one found; vertex rank/1000",100,0,10);
00154 hA[4]=new TH1F("Ln3","vertex rank of 3rd vertex if more than one found; vertex rank/1000",100,0,10);
00155 hA[5]=new TH1F("zEr","vertex Z-error; error Z (cm)",1000,0,1.);
00156 TH2F *h2=new TH2F("zErL","Rank vs. vertex Z-error ; error Z (cm); rank/1000",100,0,0.5,100,0,20.);
00157 hA[6]=(TH1F*) h2;
00158 hA[7]=new TH1F("Ln2d1","ratio of vertex rank 2nd/1st ; rank2/rank1",100,0,1.);
00159 }
00160
00161
00162
00163 pl(int page=1){
00164
00165 fd=new TFile("out.hist.root");
00166 assert(fd->IsOpen());
00167 fd->ls();
00168
00169 TString tt;
00170
00171
00172 switch(page){
00173 case 1:{
00174 tt="nPTv";
00175 } break;
00176
00177 case 2:{
00178 tt="nPFv";
00179 } break;
00180
00181 default:;
00182 }
00183
00184 c=new TCanvas();
00185 c->Divide(2,2);
00186 int iv;
00187 for(iv=0;iv<4;iv++){
00188 TString tt1=tt;
00189 tt1+=iv;
00190 printf("=%s=\n",tt1.Data());
00191
00192 h=(TH1F*)fd->Get(tt1); assert(h);
00193 c->cd(iv+1);
00194 h->Draw();
00195 }
00196
00197 }
