StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMuLcp2TreeMaker.cxx
1 
2 // *-- Author : Jan Balewski
3 //
4 // $Id: StMuLcp2TreeMaker.cxx,v 1.7 2009/12/02 16:35:58 fine Exp $
5 
6 #include <TFile.h>
7 #include <TH2.h>
8 
9 #include "StMuLcp2TreeMaker.h"
10 
11 #include "TChain.h"
12 #include "TClonesArray.h"
13 #include "StL0Trigger.h"
14 #include "StEventInfo.h"
15 #include "StEventSummary.h"
16 #include "StarClassLibrary/StThreeVectorF.hh"
17 #include "StMuDSTMaker/COMMON/StMuEvent.h"
18 #include "StMuDSTMaker/COMMON/StMuTrack.h"
19 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
20 //#include "MikesRejector/ExampleUsage.h"
21 
22 #include "CtbMatching.h"
23 
24 #include "StThreeVector.hh"
25 
26 ClassImp(StMuLcp2TreeMaker)
27 
28 
29 //________________________________________________
30 //________________________________________________
31 StMuLcp2TreeMaker::StMuLcp2TreeMaker(const char* self ,const char* muDstMakerName) : StMaker(self){
32  mMuDstMaker = (StMuDstMaker*)GetMaker(muDstMakerName);
33  assert(mMuDstMaker);
34 
35  // rejector=new ExampleUsage; // Mike's cosmic/laser rejector
36  primTrA=0;
37 
38  runID=101;// default not initialized value
39  treeName="./";
40  eve_cosm=-8;
41  SetMinNFitPoint(15);
42  SetMaxDCAxy(3.);// cm
43  SetMaxEta(1.4);
44  SetMinPt(0.4); // GeV/c
45  SetMinFitPfrac(0.55);
46  SetMaxZvert(100.); // (cm)
47  ctb=new CtbMatching();
48  eve_nGlob=-1;
49 }
50 
51 
52 //___________________ _____________________________
53 //________________________________________________
54 StMuLcp2TreeMaker::~StMuLcp2TreeMaker(){
55  if(tree) tree->Print();
56 
57  // Save all objects in this file
58  hfile->Write();
59 
60 }
61 
62 
63 //________________________________________________
64 //________________________________________________
65 Int_t StMuLcp2TreeMaker::InitRunFromMake (int runNumber){
66  float Pi=3.14159;
67  if(runID==runNumber) return kStOK;
68  if(runID==888999) return kStOK;// special run number
69 
70  if(runID !=101) {
71  printf("\n\n FATAL %s::InitRunFromMake() runID=%d changed to %d, (not implemented),STOP \n",GetName(),runID,runNumber);
72  assert(runID==runNumber);
73  }
74 
75  runID=runNumber;
76 
77  treeName+="/R";
78  treeName+=runID;
79  treeName+=".tree.root";
80  printf("%s::InitRunFromMake() \n runID=%d Setup output tree & histo to '%s' , use BXoff48=%d , CUTS: minNFitPoint=%d maxDCAxy/cm=%.2f maxEta=%.2f minPt/GeV/c=%.2f, minFitPfrac=%.2f maxZvert=%f\n",GetName(),runID,treeName.Data(),off48,C_minNFitPoint,C_maxDCAxy,C_maxEta,C_minPt,C_minFitPfrac,C_maxZvertex);
81 
82  hfile=new TFile(treeName,"RECREATE"," histograms & trees with LCP ");
83  assert(hfile->IsOpen());
84 
85  // Create some histograms and a profile histogram
86  h[0]=new TH1F("bX7","rates vs. 7bit bXing",129,-0.5,128.5);
87  h[1]=new TH1F("bX48","rates vs. 48bit bXing +offset%120 ",129,-0.5,128.5);
88  h[2]=new TH1F("bX120","rates vs. true bXing [0-119]",129,-0.5,128.5);
89  h[3]=(TH1F *)new TH2F("sBit120","spinBits vs. bXing120",129,-0.5,128.5,17,-0.5,16.5);
90  h[4] = new TH1F("Vz","Vertex Z/cm",100,-250,250.);
91  h[5] = new TH1F("nPrim","No. of prim TPC tracks",50,-0.5,49.5);
92  h[6] =new TH1F("eta","eta of LCP",100,-2.,2.);
93 
94  h[7] = new TH1F("phi"," phi/rad of LCP", 240,-Pi,Pi);
95  h[8] = new TH1F("pT","pT/GeV of LCP",100,0,10.);
96  h[9] = new TH1F("pTm1","pT/GeV of LCP, CTB match any",100,0,10.);
97  h[10] = new TH1F("pTm2","pT/GeV of LCP, CTB match LCP",100,0,10.);
98 
99 
100  if (runNumber>100) {
101  // Create tree
102  tree = new TTree("T-LCP","LCP selected from MinB trigger");
103  tree->Branch("bx48", &eve_bx48,"bx48/I");
104  tree->Branch("bx120",&eve_bx120,"bx120/I");
105  tree->Branch("id", &eve_id,"id/I");
106  tree->Branch("sb", &eve_sb,"sb/I");
107  tree->Branch("nPrim",&eve_nPrim,"nPrim/I");
108  tree->Branch("nGlob",&eve_nGlob,"nGlob/I");
109  tree->Branch("vz", &eve_vz ,"vz/F");
110  tree->Branch("cosm", &eve_cosm ,"cosm/F");
111 
112  tree->Branch("pt",&lcp_pt ,"pt/F");
113  tree->Branch("phi",&lcp_phi ,"phi/F");
114  tree->Branch("eta",&lcp_eta ,"eta/F");
115  tree->Branch("q",&lcp_q ,"q/I");
116  tree->Branch("nFit",&lcp_nFit ,"nFit/I");
117  } else {
118  printf(" %s::InitRunFromMake(RunNumber=%d) WARN : TTRee NOT created\n",GetName(),runNumber);
119  tree=0;
120  }
121 
122 #if 0
123  // cuts study
124  int nd=0; float x1=1,x2=2; TString titCore="fixMe2";
125 
126  // titCore="nFitPoint"; nd=15; x1=9.5; x2=39.5;
127  // titCore="maxDCAxy/cm"; nd=15; x1=0.; x2=3.0;
128  // titCore="maxEta"; nd=15; x1=0.35; x2=1.85;
129  titCore="minFitPfrac"; nd=10; x1=0.325; x2=0.825;
130 
131  hc[0]=new TH1F("CnFP-A","Yield vs."+titCore+", for LCP #Delta#phi <1/8 #pi ",nd,x1,x2);
132  hc[1]=new TH1F("CnFP-B","Yield vs."+titCore+", for LCP #Delta#phi #in [ 1/8 #pi, 3/8 #pi ] ",nd,x1,x2);
133  hc[2]=new TH1F("CnFP-C","Yield vs."+titCore+", for LCP #Delta#phi #in [ 3/8 #pi, 5/8 #pi ] ",nd,x1,x2);
134  hc[3]=new TH1F("CnFP-D","Yield vs. "+titCore+", for LCP #Delta#phi #in [ 5/8 #pi, 7/8 #pi ] ",nd,x1,x2);
135  hc[4]=new TH1F("CnFP-E","Yield vs."+titCore+", for LCP #Delta#phi > 7/8 #pi ",nd,x1,x2);
136  hc[5]=new TH1F("CnFP-L","Yield vs."+titCore+", lost LCP",nd,x1,x2);
137  hc[6]=new TH1F("CnFP-W","Yield vs."+titCore+", won LCP",nd,x1,x2);
138 #endif
139 
140 
141  return kStOK;
142 }
143 
144 //________________________________________________
145 //________________________________________________
146 Int_t StMuLcp2TreeMaker::Init(){
147  return StMaker::Init();
148 }
149 
150 
151 //________________________________________________
152 //________________________________________________
153 void StMuLcp2TreeMaker::clearLCP(){
154  lcp_eta=lcp_phi=lcp_pt=0;
155  lcp_q=0;
156  lcp_nFit=0;
157 }
158 
159 
160 //________________________________________________
161 //________________________________________________
163 
164  static int kEve=0;
165  kEve++;
166 
167  printf("%s::Make() is called ..........\n",GetName());
168 
169  StMuDst* dst = mMuDstMaker->muDst();
170  primTrA=dst->primaryTracks();
171 
172  StMuEvent* muEve = dst->event();
173  StL0Trigger &trig=muEve->l0Trigger();
174  StEventInfo &info=muEve->eventInfo();
175  StEventSummary &smry=muEve->eventSummary();
176  ctb->loadHits(muEve);
177  StThreeVectorF ver=smry.primaryVertexPosition();
178  if(runID>100) InitRunFromMake(info.runId());
179 
180  // printf("trgWrd=%d nPrim=%d, Vz=%f \n", trig.triggerWord(),primTrA->GetEntries(),ver.z());
181  // use only minB trigger events
182  if(trig.triggerWord()!=0x2000) return kStOK;
183 
184  // reject events without primtracks
185  if(primTrA->GetEntries()<=0) return kStOK;
186 
187  // fill eve-related tree values
188  eve_vz=ver.z();
189  if(fabs(eve_vz)>C_maxZvertex) return kStOK;
190 
191  int bx7=trig.bunchCrossingId7bit(info.runId());
192  eve_id=info.id();
193  eve_bx48 = trig.bunchCrossingId();
194  eve_bx120 =(trig.bunchCrossingId()+off48 )%120;
195  eve_sb=trig.spinBits(info.runId());
196 
197  // test for false vertex
198  //eve_cosm=rejector->acceptEvent(dst);
199  eve_cosm=0;
200 
201  // search for LCP
202  clearLCP();
203  const StMuTrack* lcp= findLCP(C_minPt,C_minNFitPoint, C_maxDCAxy, C_maxEta,C_minFitPfrac);
204 
205  // fill histo 1:1
206  h[0]->Fill(bx7);
207  h[1]->Fill(eve_bx48);
208  h[2]->Fill(eve_bx120);
209  ((TH2F*)h[3])->Fill(eve_bx120,eve_sb);
210  h[4]->Fill(eve_vz);
211  h[5]->Fill(eve_nPrim); // <<-- assigned by findLCP()
212 
213  if(lcp) { // LCP was found
214  lcp_eta=lcp->eta();
215  lcp_phi=lcp->phi();
216  lcp_pt= lcp->pt();
217  lcp_q = lcp->charge();
218  lcp_nFit=lcp->nHitsFit();
219 
220  // fill those histo only for LCP
221  h[6]->Fill(lcp_eta);
222  h[7]->Fill(lcp_phi);
223  h[8]->Fill(lcp_pt);
224  if(eve_cosm) h[9]->Fill(lcp_pt);
225  if(eve_cosm>1000) h[10]->Fill(lcp_pt);
226  }
227 
228  // scan for good globals, to match Olg's analysis
229  {
230  // TClonesArray* globTrA=dst->globalTracks();
231  TObjArray* globTrA=dst->primaryTracks();
232  //printf("tot Glob=%d\n", globTrA->GetEntries());
233  int nTr=0;
234  for (int i=0; i<globTrA->GetEntries();i++) {
235  TObject &oo=globTrA[i];
236  const StMuTrack* gTr = (StMuTrack*)(&oo);// perhaps it works? JB
237  if(gTr->flag() <=0) continue;
238  StThreeVectorF dca=gTr->dcaGlobal();
239  if(!gTr->topologyMap().trackTpcOnly()) continue;
240  if(gTr->nHitsFit()<15 )continue;
241  if(dca.perp()>3.) continue;
242  if(fabs(dca.z())>3.) continue;
243  //printf("nFit=%d,DCA x=%f y=%f z=%f zV=%f\n",gTr->nHitsFit(),dca.x(),dca.y(),dca.z(),ver.z());
244  nTr++;
245  }
246  eve_nGlob=nTr;
247  }
248 
249 
250  if(tree)tree->Fill();
251  //if(kEve%500==0)
252  printf("#%s::Make(%d) nPrim=%d nGlob=%d lcp: pt=%f nFit=%d eta=%f\n",GetName(),kEve,eve_nPrim, eve_nGlob, lcp_pt,lcp_nFit,lcp_eta);
253 
254  // extra tests of cuts
255  // examinCut(lcp);
256 
257  return kStOK;
258 }
259 
260 //________________________________________________
261 //________________________________________________
262 const StMuTrack* StMuLcp2TreeMaker::findLCP( float XminPt,int XminNFitP , float XmaxDCAxy, float XmaxEta, float XminFitPfrac) {
263 
264  eve_cosm=0;
265  const StMuTrack* lp=0;
266  int nTr=primTrA->GetEntries();
267  int nTr0=0,nTr1=0; // counts valid prim TPC tracks
268  int lpMatch=0;
269  float maxPt=XminPt;// init pT
270  for (int i=0; i<nTr; i++) {
271  TObject &oo=primTrA[i];
272  const StMuTrack* pTr = (StMuTrack*)(&oo);// perhaps it works? JB
273  if(pTr->flag() <=0) continue;
274  if(!pTr->topologyMap().trackTpcOnly()) continue;
275  if(pTr->globalTrack()->nHitsFit()>15 ) nTr0++; // was just: nTr0++;
276  // printf("itr=%d pT=%f nFit=%d\n",i, pTr->pt(),pTr->nHitsFit());
277 
278  StThreeVectorF dca=pTr->dcaGlobal();
279  //printf("xyz=%f %f %f, pt=%f\n",dca.x(),dca.y(),dca.z(),dca.perp());
280  if(dca.perp()>XmaxDCAxy) continue;
281 
282  // switch to oryginal glob-track, since vertex constrain biases pt
283  const StMuTrack* gTr=pTr->globalTrack();
284  assert(gTr);// every prim must have its global
285 
286  if ( gTr->nHitsFit() < XminNFitP) continue;
287  if ( fabs(gTr->eta()) > XmaxEta ) continue;
288  float frac=1.*gTr->nHitsFit()/gTr->nHitsPoss();
289  if ( frac < XminFitPfrac) continue;
290  nTr1++;
291  // printf("itr=%d pT=%f phi/deg=%f eta=%f \n",i, pTr->pt(),pTr->phi()/3.1416*180,pTr->eta());
292 
293  int match=ctb->match(gTr);
294  if(match>0) eve_cosm++;
295  // printf("itr=%d pT=%f nFit=%d match=%d\n",i, pTr->pt(),pTr->nHitsFit(),match);
296  // now search for high pT
297  if ( gTr->pt() < maxPt ) continue;
298  maxPt=gTr->pt();
299  lp=gTr;
300  if(match>0)lpMatch++;
301  }
302  // tmp - apply pT cut on LCP
303  //if(maxPt<1.5 || maxPt>3) lp=0;
304 
305  eve_nPrim=nTr0;
306  if(lpMatch) eve_cosm+=1000;
307  // printf("nTr0=%d nTr1=%d lp=%p match=%d\n",nTr0,nTr1,lp,eve_cosm);
308  return lp;
309 }
310 
311 
312 //________________________________________________
313 //________________________________________________
314 void StMuLcp2TreeMaker::examinCut(const StMuTrack*lcp0){
315 
316  const float Pi=3.1416;
317 
318  StThreeVectorF p0;
319  float phi0=0;
320 
321  if(lcp0) {
322  p0=lcp0->p();
323  phi0=Pi/2+ atan2(p0.y(),p0.x());
324  }
325 
326  int nb=hc[0]->GetNbinsX();
327 
328  int i;
329  for(i=1;i<=nb;i++) {// vary nFitP ......................
330  float x=hc[0]->GetBinCenter(i);
331  const StMuTrack* lcp1=0;
332 
333  // lcp1=findLCP(C_minPt,(int)x,C_maxDCAxy,C_maxEta,C_minFitPfrac); // test nFitPoint
334  //lcp1=findLCP(C_minPt,C_minNFitPoint,x,C_maxEta,C_minFitPfrac); // test maxDCAxy
335  // lcp1=findLCP(C_minPt,C_minNFitPoint,C_maxDCAxy,x,C_minFitPfrac); // test maxEta
336  lcp1=findLCP(C_minPt,C_minNFitPoint,C_maxDCAxy,C_maxEta,x); // test fitPfrac
337  //printf("i=%d x=%f %p %p\n",i,x,lcp0,lcp1);
338 
339  float delPhi=Pi*1.1;
340  if(lcp0 && lcp1==0) { // lost LCP
341  hc[5]->Fill(x);
342  } else if(lcp1 && lcp0==0) { // won LCP
343  hc[6]->Fill(x);
344  } else if(lcp1 && lcp0) { // phi could change
345  StThreeVectorF p1=lcp1->p();
346  float phi1=Pi/2+ atan2(p1.y(),p1.x());
347  delPhi=phi1-phi0;
348  if(delPhi <0) delPhi*=-1;
349  if(delPhi >Pi) delPhi=2*Pi-delPhi;
350  // printf("(deg) phi0=%.1f, phi1=%.1f, del=%.1f\n",phi0/Pi*180, phi1/Pi*180, delPhi/Pi*180);
351  if( delPhi <1./8*Pi) {
352  hc[0]->Fill(x);
353  } else if( delPhi <3./8*Pi) {
354  hc[1]->Fill(x);
355  } else if( delPhi <5./8*Pi) {
356  hc[2]->Fill(x);
357  } else if( delPhi <7./8*Pi) {
358  hc[3]->Fill(x);
359  } else {
360  hc[4]->Fill(x);
361  }
362  }
363  }
364 }
365 
366 
367 
368 
369 // $Log: StMuLcp2TreeMaker.cxx,v $
370 // Revision 1.7 2009/12/02 16:35:58 fine
371 // Fix StMuTrack interface
372 //
373 // Revision 1.6 2005/10/11 17:28:30 balewski
374 // fix related L0trig code change
375 //
376 // Revision 1.5 2005/08/23 21:09:23 balewski
377 // fix to follow muDst evolution
378 //
379 // Revision 1.4 2004/01/06 17:25:26 balewski
380 // get LCP from Geant info
381 //
382 // Revision 1.3 2003/11/12 18:43:41 balewski
383 // final for LCP paper
384 //
385 // Revision 1.2 2003/10/20 17:04:39 balewski
386 // LCP analysis code
387 //
388 // Revision 1.1 2003/09/16 19:18:35 balewski
389 // extraction of LCP from muDst
390 //
StMuDst * muDst()
Definition: StMuDstMaker.h:425
Double_t pt() const
Returns pT at point of dca to primary vertex.
Definition: StMuTrack.h:256
virtual Int_t Make()
UShort_t nHitsFit() const
Return total number of hits used in fit.
Definition: StMuTrack.h:239
short flag() const
Returns flag, (see StEvent manual for type information)
Definition: StMuTrack.h:230
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
const StThreeVectorF & p() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:259
Double_t eta() const
Returns pseudo rapidity at point of dca to primary vertex.
Definition: StMuTrack.h:257
static TObjArray * primaryTracks()
returns pointer to a list of tracks belonging to the selected primary vertex
Definition: StMuDst.h:301
UShort_t nHitsPoss() const
Return number of possible hits on track.
Definition: StMuTrack.cxx:242
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
Definition: Stypes.h:40
This commented block at the top ...
Double_t phi() const
Returns phi at point of dca to primary vertex.
Definition: StMuTrack.h:258
const StMuTrack * globalTrack() const
Returns pointer to associated global track. Null pointer if no global track available.
Definition: StMuTrack.h:272
StThreeVectorF dcaGlobal(Int_t vtx_id=-1) const
Returns 3D distance of closest approach to primary vertex of associated global track.
Definition: StMuTrack.cxx:371
StTrackTopologyMap topologyMap() const
Returns topology map.
Definition: StMuTrack.h:254