StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
rdMu2Print.C
1 #include "StRoot/StMuDSTMaker/COMMON/StMuArrays.h"
2 
3 class StChain;
4 class StMuEmcCollection;
5 class StMuTriggerIdCollection;
6 StChain *chain=0;
7 
8 
9 
10 int rdMu2Print(
11  char* file = "QCDprodMBa.MuDst.root",
12  int nEve=1,
13  char* inDir = "./"
14  )
15 {
16  Int_t nFiles = 1;
17 
18  // inDir="/star/data05/scratch/balewski/fgtAna_Mar17/";
19  file="muonEta2PhiAB_pt20.MuDst.root";
20 
21  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
22  loadSharedLibraries();
23  gSystem->Load("StFgtUtil");
24 
25  cout << " loading done " << endl;
26 
27 
28  // create chain
29  chain = new StChain("StChain");
30 
31  // Now we add Makers to the chain...
32  muMk = new StMuDstMaker(0,0,inDir,file,"MuDst.root",nFiles);
33  TChain* tree=muMk->chain(); assert(tree);
34  int nEntries=(int) tree->GetEntries();
35  printf("total eve in muDst chain =%d\n",nEntries); // return ;
36  if(nEntries<0) return;
37 
38 
39  chain->Init();
40  chain->ls(3);
41 
42  h1=new TH1F("nP","# of prim tracks per vertex",100,0,200);
43  h2=new TH1F("vR","PPV Vertex rank; rank", 300, -1.2e6, 1.2e6); h2->SetLineColor(kRed);
44  h3=new TH1F("vRL","PPV Vertex rank, funny X-axis; X=Log10(rank-1e6)+ offset", 150, -9,21);
45  h2->GetXaxis()->SetTitleSize(.043);
46  h2->GetXaxis()->SetTitleSize(.043);
47 
48  h4=new TH1F("trPhi"," prim tracks phi if PT>1.0 GeV/c; phi (rad)",50,-3.2,3.2);
49  h5=new TH1F("zVerTrg"," Z trg-vertex , rank>0; Z (cm)", 50,-200,200);
50  h6=new TH1F("zVerPlp"," Z pileup-vertex , rank<0; Z (cm)", 50,-200,200);
51  h7=new TH1F("nPrV","# of prim tr used by PPV; rank>0",15,0,30);
52  //---------------------------------------------------
53  int eventCounter=0;
54  int t1=time(0);
55 
56  for (Int_t iev=0;iev<nEntries; iev++) {
57  if(eventCounter>=nEve) break;
58  chain->Clear();
59  int stat = chain->Make();
60  if(stat) break; // EOF or input error
61  //printf("stat=%d\n", stat);
62  eventCounter++;
63  // if(eventCounter<17) continue;
64 
65  // Access to muDst .......................
66  StMuEvent* muEve = muMk->muDst()->event();
67 
68  StEventInfo &info=muEve->eventInfo();
69  int nPrimV=muMk->muDst()->numberOfPrimaryVertices();
70  StMuTriggerIdCollection *tic=&(muEve->triggerIdCollection());
71 
72  int trigID=96211;
73  bool fired=tic->nominal().isTrigger(trigID);
74 
75  Int_t nGlobTrAll=muMk->muDst()->GetNGlobalTrack();
76 
77  //assert(tic.nominal().isTrigger(127271)==0);
78  printTrig(tic);
79 
80  if(eventCounter%1==0) {
81  printf("ieve=%d eventID %d nPrimV=%d nGlobTrAll=%d =============\n", eventCounter,info.id(),nPrimV,nGlobTrAll);
82  // printf("TrigID=%d fired=%d\n",trigID,fired);
83 
84  // if(nPrimV>1) printf("######\n");
85  }
86 
87  int iv;
88  if(1)for(iv=0;iv<nPrimV;iv++) {
89  StMuPrimaryVertex* V= muMk->muDst()->primaryVertex(iv);
90  assert(V);
91  muMk->muDst()->setVertexIndex(iv);
92  StThreeVectorF &r=V->position();
93  StThreeVectorF &er=V->posError();
94  printf("iv=%d Vz=%.2f +/-%.2f \n",iv,r.z(),er.z() );
95  // count prim tracks for this vert
96  int nPrimTr =0;
97  int itr;
98  Int_t nPrimTrAll=muMk->muDst()->GetNPrimaryTrack();
99  for(itr=0;itr<nPrimTrAll;itr++) {
100  StMuTrack *pr_track=muMk->muDst()->primaryTracks(itr);
101  if(pr_track->flag()<=0) continue;
102  nPrimTr ++;
103  }
104  if(nPrimV>0)h1->Fill(nPrimTr);
105  float rank=V->ranking();
106  h2->Fill(rank);
107  if(rank>1e6) h3->Fill(log(rank-1e6)+10);
108  else if(rank>0) h3->Fill(log(rank));
109  else h3->Fill(log(rank+1e6)-10);
110 
111  if(1)printf(" nPrimTr=%d , Z=%.1f VFid=%d:: ntrVF=%d nCtb=%d nBemc=%d nEEmc=%d nTpc=%d sumPt=%.1f rank=%g\n"
112  ,nPrimTr,r.z(), V->vertexFinderId() ,V->nTracksUsed() ,V->nCTBMatch() ,V-> nBEMCMatch() ,V->nEEMCMatch() ,V->nCrossCentralMembrane() ,V->sumTrackPt() ,V->ranking());
113  if (rank>0) h5->Fill(r.z());
114  if (rank>0) h7->Fill(V->nTracksUsed());
115  if (rank<0) h6->Fill(r.z());
116  }
117 
118 
119  // continue; // do NOT print prim tracks for each vertex
120 
121  if(1)
122  for(iv=0;iv<nPrimV;iv++) {
123  muMk->muDst()->setVertexIndex(iv);
124  Int_t nPrimTrAll=muMk->muDst()->GetNPrimaryTrack();
125  cout<<"\n Prim "<<nPrimTrAll<<" tracks belonging to "<<iv<<"-th prim vertex"<<endl;
126 
127  int itr;
128  int ntr=0;
129 
130  for(itr=0;itr<nPrimTrAll;itr++) {
131  StMuTrack *pr_track=muMk->muDst()->primaryTracks(itr);
132  if(pr_track->flag()<=0) continue;
133  ntr++;
134  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;
135  cout << "\t primV("<<iv<<") primDCA=" << pr_track->dca(iv) << endl;
136  if(pr_track->dca(iv).mag()>5) cout << "^^^^^ 3D DCA magnitude="<<pr_track->dca(iv).mag()<<endl;
137  // cout << "\t first point " << pr_track->firstPoint() << endl;
138  // cout << "\t last point " << pr_track->lastPoint() << endl;
139  if(pr_track->pt()>1.0) h4->Fill(pr_track->phi());
140 
141  } // end of loop over tracks
142 
143  }// end of loop over vertices
144 
145  if(1) { // EMC response
146 
147  StMuEmcCollection* emc = muMk->muDst()->muEmcCollection();
148  if (!emc) {
149  printf(" No EMC data for this event\n");
150  return kStOK;
151  }
152  printEEtower(emc);
153  // printEEpre(emc);
154  // printEEsmd(emc);
155  //printBEtower(emc);
156  //printBEpre(emc);
157  //printBEsmd(emc);
158  }
159 
160 
161  if(1) { // fgt hits
162  printFgtStrips(muMk->muDst()->fgtArray( muFgtStrips));
163  printFgtClusters(muMk->muDst()->fgtArray( muFgtClusters ));
164  // printFgtInfo(muMk->muDst()->fgtArray( muFgtInfo));
165  }
166 
167 
168  }
169  printf("****************************************** \n");
170  // return;
171  int t2=time(0);
172  if(t2==t1) t2=t1+1;
173  float tMnt=(t2-t1)/60.;
174  float rate=1.*eventCounter/(t2-t1);
175  printf("sorting done %d of nEve=%d, CPU rate=%.1f Hz, total time %.1f minute(s) \n\n",eventCounter,nEntries,rate,tMnt);
176  c=new TCanvas(); c->Divide(2,2);
177  c->cd(1); h5->Fit("gaus"); //h5->Draw();
178  c->cd(3); h6->Fit("gaus"); //h5->Draw();
179  c->cd(2); h3->Draw();
180  c->cd(4); h7->Draw();
181 
182  return;
183 
184 }
185 
186 
187 //===========================================
188 //===========================================
189  void printFgtStrips( TClonesArray * fgtStrips){
190  assert( fgtStrips );
191 
192  Int_t nStrips = fgtStrips->GetEntriesFast();
193  printf("FGT nStrips=%d\n",nStrips);
194  for( Int_t i = 0; i < nStrips; ++i ){
195  StMuFgtStrip* strip = static_cast< StMuFgtStrip* >( (*fgtStrips)[i] );
196  if(strip==0) continue;
197  Int_t geoId = strip->getGeoId();
198  Short_t seedType=strip->getClusterSeedType();
199  printf("i=%d geoId=%d charge=%.1f +/-%.1f seedType=%d\n",i,geoId, strip->getCharge(), strip->getChargeUncert(),seedType);
200  }
201  }
202 
203 //===========================================
204 //===========================================
205 void printFgtClusters( TClonesArray * fgtClusters){
206  assert( fgtClusters );
207 
208  Int_t nClusters = fgtClusters->GetEntriesFast();
209  printf("FGT nClusters=%d\n",nClusters);
210  for( Int_t i = 0; i < nClusters; ++i ){
211  StMuFgtCluster* clus = static_cast< StMuFgtCluster* >( (*fgtClusters)[i] );
212  if(clus==0) continue;
213  Int_t geoId = clus->getCentralStripGeoId();
214  printf("i=%d cntrGeoId=%d charge=%.1f R=%.1f+/-%.1f phi=%.3f+/-%.3f nStrip=%d\n",i,geoId , clus->getCharge(), clus->getR(),clus->getErrR(), clus->getPhi(), clus->getErrPhi(), clus->getNumStrips());
215 
216  }
217 
218 }
219 
220 //===========================================
221 //===========================================
222 void printFgtInfo( TClonesArray *fgtInfo){
223  assert( fgtInfo );
224  Int_t nInfo = fgtInfo->GetEntriesFast();
225  printf("FGT nInfo=%d\n",nInfo);
226  for( Int_t i = 0; i < nInfo; ++i ){
227  StMuFgtInfo* info = static_cast< StMuFgtInfo* >( (*fgtInfo)[i] );
228  if(info==0) continue;
229  //printf("i=%d info=%s=\n",i,info->mMsg);
230  }
231 
232 }
233 
234 //===========================================
235 //===========================================
236 printEEtower( StMuEmcCollection* emc ) {
237  int sec,eta,sub,adc;
238  StMuEmcHit *hit;
239 
240  int i, nh;
241 
242  printf("\Total %d hits in Tower (only ADC>0)\n",emc->getNEndcapTowerADC());
243  nh=0;
244  for (i=0; i< emc->getNEndcapTowerADC(); i++) {
245  emc->getEndcapTowerADC(i,adc,sec,sub,eta);
246  if (adc<=0) continue; // print only non-zero values
247  nh++;
248  printf("i=%d Tower %2.2dT%c%2.2d adc=%4d\n",i,sec,sub+'A'-1,eta,adc );
249  // printf(" Tower isec=%d ieta=%d isub=%d adc=%4d\n",sec,eta, sub,adc );
250  int adcX=1000+ (eta-1) + (sub-1)*12 +(sec-1)*60;
251  // assert(adc==adcX );
252 }
253  printf(" Total %d E-Towers with ADC>0\n",nh);
254 }
255 
256 
257 //===========================================
258 //===========================================
259 printEEpre( StMuEmcCollection* emc ) {
260  int sec,eta,sub,pre,adc;
261  StMuEmcHit *hit;
262 
263  int i, nh;
264  nh= emc->getNEndcapPrsHits();
265  printf("\nTotal %d hits in pre1+2+post\n",nh);
266  for (i=0; i<nh; i++) {
267  hit=emc->getEndcapPrsHit(i,sec,sub,eta,pre);
268  int ss=sub + 5*(pre-1);
269  adc=hit->getAdc();
270  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);
271  int adcX= (eta-1) + (sub-1) *12 +(sec-1)*60 + 1000*pre;
272 
273  // assert(adc==adcX );
274 
275  }
276 }
277 
278 
279 //===========================================
280 //===========================================
281 printEEsmd( StMuEmcCollection* emc ) {
282  int sec,strip,adc;
283  char uv='U';
284 
285  for(uv='U'; uv<='V'; uv++) {
286  int nh= emc->getNEndcapSmdHits(uv);
287  printf("\nTotal %d hits in SMD-%c\n",nh,uv);
288  for (int i=0; i<nh; i++) {
289  hit=emc->getEndcapSmdHit(uv,i,sec,strip);
290  adc=hit->getAdc();
291  printf(" SMD-%c %2.2d%c%3.3d : energy=%f adc=%d\n",uv,sec,uv,strip,hit->getEnergy(),adc);
292  int adcX= 1000 + strip-1 +(sec-1)*300;
293  // assert(adc==adcX );
294  }
295  }
296 }
297 
298 //===========================================
299 //===========================================
300 printBEtower( StMuEmcCollection* emc ) {
301  int sec,eta,sub,adc;
302  StMuEmcHit *hit;
303 
304  int i, nh;
305 
306  printf("\Total hits in BTower (only ADC>0)\n");
307  nh=0;
308  for (i=0; i< 4800; i++) {
309  int adc = emc->getTowerADC(i);
310  if (adc<=1000) continue; // print only non-zero values
311  nh++;
312  printf(" Tower id=%d adc=%4d\n",i,adc );
313  // assert(adc==adcX );
314 }
315  printf(" --> %d towers with ADC>thr\n",nh);
316 }
317 
318 
319 
320 
321 //===========================================
322 //===========================================
323 printBEpre( StMuEmcCollection* emc ) {
324 
325  int i, nh;
326  nh = emc->getNPrsHits();
327  printf("\nTotal %d hits in pre1\n",nh);
328 
329  int n1=0;
330  for (i=0; i<nh; i++) {
331  StMuEmcHit * hit = emc->getPrsHit(i);
332  int adc = hit->getAdc();
333  int id=hit->getId();
334  if(adc<4) continue;
335  n1++;
336  printf("BPRS i=%d, id=%d adc = %d\n",i,id,adc);
337  }
338  printf(" --> %d BPRS hits with ADC>thr\n",n1);
339 }
340 
341 //===========================================
342 //===========================================
343 printBEsmd( StMuEmcCollection* emc ) {
344 
345  int n1=0,n2=0;
346  int nh = emc->getNSmdHits(3);
347  printf("\nTotal %d hits in SMDE\n",nh);
348  for (int i=0; i<nh; i++) {
349  StMuEmcHit * hit = emc->getSmdHit(i,3);
350  adc = hit->getAdc();
351  int id=hit->getId();
352  if(adc<4) continue;
353  n1++;
354  printf("BSMDE i=%d, id=%d adc = %d\n",i,id,adc);
355  }
356 
357  int nh = emc->getNSmdHits(4);
358  printf("\nTotal %d hits in SMDP\n",nh);
359  for (int i=0; i<nh; i++) {
360  hit = emc->getSmdHit(i,4);
361  adc = hit->getAdc();
362  int id=hit->getId();
363  if(adc<4) continue;
364  n2++;
365  printf("BSMDP i=%d,id=%d adc = %d\n",i,id,adc);
366  }
367 
368  printf(" --> %d BSMD-E & %d BSMD-P hits with ADC>thr\n",n1,n2);
369 }
370 
371 //--------------------------------------
372 void printTrig(StMuTriggerIdCollection *tic){
373 
374  const StTriggerId &l1=tic->l1();
375  vector<unsigned int> idL=l1.triggerIds();
376  printf("nTrig=%d, trigID: ",idL.size());
377  for(unsigned int i=0;i<idL.size(); i++){
378  printf("%d, ",idL[i]);
379  }
380  printf("\n");
381 }
382 
StMuDst::primaryVertex
static StMuPrimaryVertex * primaryVertex()
return pointer to current primary vertex
Definition: StMuDst.h:322
StMuEmcHit::getId
int getId() const
Return Module number.
Definition: StMuEmcHit.h:18
StMuTrack::vertexIndex
Int_t vertexIndex() const
Returns index of associated primary vertex.
Definition: StMuTrack.cxx:600
StMuEmcHit::getEnergy
float getEnergy() const
Return Hit energy.
Definition: StMuEmcHit.h:21
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
StMuTrack::pt
Double_t pt() const
Returns pT at point of dca to primary vertex.
Definition: StMuTrack.h:256
StChain::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StChain.cxx:77
StMuEmcHit::getAdc
int getAdc() const
Return ADC value.
Definition: StMuEmcHit.h:19
StMuTrack::flag
short flag() const
Returns flag, (see StEvent manual for type information)
Definition: StMuTrack.h:230
StMuTrack::charge
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
StMuTrack::p
const StThreeVectorF & p() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:259
TDataSet::ls
virtual void ls(Option_t *option="") const
Definition: TDataSet.cxx:495
StMuDst::primaryTracks
static TObjArray * primaryTracks()
returns pointer to a list of tracks belonging to the selected primary vertex
Definition: StMuDst.h:301
StChain::Make
virtual Int_t Make()
Definition: StChain.cxx:110
StMuDst::muEmcCollection
static StMuEmcCollection * muEmcCollection()
returns pointer to current StMuEmcCollection
Definition: StMuDst.h:389
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
kStOK
Definition: Stypes.h:40
StMuTrack::nHits
UShort_t nHits() const
Bingchu.
Definition: StMuTrack.h:237
StMuTrack::dca
StThreeVectorF dca(Int_t vtx_id=-1) const
Returns 3D distance of closest approach to primary vertex.
Definition: StMuTrack.cxx:359
StMuDstMaker::chain
TChain * chain()
In read mode, returns pointer to the chain of .MuDst.root files that where selected.
Definition: StMuDstMaker.h:426
StMuTrack::phi
Double_t phi() const
Returns phi at point of dca to primary vertex.
Definition: StMuTrack.h:258
StMuDst::fgtArray
static TClonesArray * fgtArray(int type)
returns pointer to the n-th TClonesArray from the fgt arrays
Definition: StMuDst.h:293
StMuDst::setVertexIndex
static void setVertexIndex(Int_t vtx_id)
Set the index number of the current primary vertex (used by both primaryTracks() functions and for St...
Definition: StMuDst.cxx:273
StMuTriggerIdCollection
Collection of trigger ids as stored in MuDst.
Definition: StMuTriggerIdCollection.h:31
Generated by   doxygen 1.8.5