StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
find_vertex.C
1 //
2 // macro: find_vertex.C
3 //
4 // author: Marco van Leeuwen
5 //
6 // notes: placed in CVS by Gene Van Buren
7 // Need to load libs before running this macro,
8 // because of function with StEvent, StPrimaryVertex
9 // When you run this, it will actually add
10 // new vertices to the event, so you have to
11 // keep track of how many were there when
12 // you read the file at first; I think the
13 // new vertices are added after the old one.
14 //
15 
16 //class StEmcDetector;
17 class StChain;
18 //class StEmcRawData;
19 //class EEfeeDataBlock;
20 //class Collection;
21 //class StSPtrVecTrackNodeIterator ;
22 //class StTriggerIdCollection;
23 class StMinuitVertexFinder;
24 class StGenericVertexFinder;
25 class StGenericVertexMaker;
26 class StEvent;
27 class StPrimaryVertex;
28 
29 StChain *chain=0;
30 TH1F *dca_z_h=0;
31 TNtuple *vtx_tuple=0;
32 
33 float eval_print_vertex(StEvent *event, StPrimaryVertex *vtx) {
34  cout << "Vertex at " << vtx->position() << " chisq " << vtx->chiSquared() << endl;
35  cout << "CTB matches: " << vtx->numMatchesWithCTB() << " BEMC: "
36  << vtx->numMatchesWithBEMC() << " cross central membrane "
37  << vtx->numTracksCrossingCentralMembrane()
38  << " rank " << vtx->ranking() << endl;
39  int cnt_trk=0;
40  int cnt_dcaz=0;
41  int cnt_dca=0;
42  Int_t n_node = event->trackNodes().size();
43  Double_t avg_dip = 0;
44  Double_t rms_dip = 0;
45  for (int i_node = 0; i_node < n_node; i_node++) {
46  StTrackNode *node = event->trackNodes()[i_node];
47  StTrack *track = node->track(global);
48  if (track &&
49  track->flag() >= 0 &&
50  track->fitTraits().numberOfFitPoints() >= 15 &&
51  !track->topologyMap().trackFtpc() &&
52  TMath::Finite(track->length()) ) {
53  //cout << "track " << track << endl;
54  cnt_trk++;
55  double pathlength = track->geometry()->helix().pathLength( vtx->position(), false ); // do not scan periods
56  StThreeVectorF dca = track->geometry()->helix().at(pathlength)-vtx->position();
57  dca_z_h->Fill(dca.z());
58  if (fabs(dca.z()) < 3) {
59  cnt_dcaz++;
60  float dip = track->geometry()->helix().dipAngle();
61  avg_dip += dip;
62  rms_dip += dip * dip;
63  }
64  if (dca.mag() < 3) {
65  cnt_dca++;
66  }
67  }
68  }
69  if (cnt_dcaz) {
70  avg_dip /= cnt_dcaz;
71  rms_dip = sqrt(rms_dip - cnt_dcaz*avg_dip*avg_dip) / cnt_dcaz;
72  }
73  cout << "Tracks used in vtx " << vtx->numTracksUsedInFinder() << ", num with dcaz cut " << cnt_dcaz << ", dca < 3 " << cnt_dca << endl;
74  cout << "avg dip " << avg_dip << ", rms dip " << rms_dip << endl;
75  vtx_tuple->Fill(vtx->numTracksUsedInFinder(),vtx->position().z(),vtx->chiSquared(),cnt_dcaz,cnt_dca,vtx->numMatchesWithBEMC(),vtx->numTracksCrossingCentralMembrane(),avg_dip,rms_dip);
76 
77  return vtx->chiSquared();
78 }
79 
80 void find_vertex(char * fname="high_053/st_physics_6053108_raw_2020002.event.root", Int_t nevents=1000){
81 
82  //fname="high_053/st_physics_6053108_raw_2020002.event.root";// daq1
83  //fname="low_mb_49/st_physics_adc_6049129_raw_3080002.event.root";
84 
85  // ROOT libs
86  gSystem->Load("libPhysics");
87  gSystem->Load("libTable");
88  gSystem->Load("libGeom");
89 
90 // STAR libs
91  gSystem->Load("StarMagField");
92  gSystem->Load("St_base");
93  gSystem->Load("StChain");
94  gSystem->Load("St_Tables");
95  gSystem->Load("StUtilities"); // new addition 22jul99
96  gSystem->Load("StTreeMaker");
97  gSystem->Load("StIOMaker");
98  gSystem->Load("StarClassLibrary");
99  gSystem->Load("StTriggerDataMaker"); // new starting from April 2003
100  gSystem->Load("StBichsel");
101  gSystem->Load("StEvent");
102  gSystem->Load("StTpcDb");
103  gSystem->Load("StEventUtilities");
104  gSystem->Load("StEmcUtil");
105  gSystem->Load("StTofUtil");
106  gSystem->Load("StPmdUtil");
107  gSystem->Load("StPreEclMaker");
108  gSystem->Load("StStrangeMuDstMaker");
109  gSystem->Load("StMuDSTMaker");
110  gSystem->Load("StMagF");
111 
112  gSystem->Load("StDbLib.so");
113  gSystem->Load("StDbBroker.so");
114  gSystem->Load("libStDb_Tables.so");
115  gSystem->Load("St_db_Maker.so");
116 
117 
118  cout << " loading of shared libraries done" << endl;
119 
120 
121  //gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
122  //loadSharedLibraries();
123 
124  // Load my makers
125 
126  gSystem->Load("Sti");
127  gSystem->Load("libStEEmcUtil");
128  gSystem->Load("libMinuit");
129  gSystem->Load("StGenericVertexMaker");
130 
131  TFile *fout = new TFile("vtx_tree.root","RECREATE");
132  dca_z_h = new TH1F("dca_z_h","dz_z_h",100,-50,50);
133  vtx_tuple = new TNtuple("vtx_tuple","vertex info","num_trk:vtx_z:chisq:n_dcaz:n_dca:n_bemc:n_cross:avg_dip:rms_dip");
134 
135  // book some histos
136  //gSystem->Load("libStGenericVertexMaker");
137  // create chain
138  chain = new StChain("bfc");
139  //chain->SetDebug();
140 
141  // Now we add Makers to the chain...
142 
143  // StIOMaker - to read files ...
144  StIOMaker* ioMaker = new StIOMaker();
145  //ioMaker->SetFile("photon_bemc.event.root");
146  ioMaker->SetFile(fname);
147  //ioMaker->SetDebug();
148  ioMaker->SetIOMode("r");
149  ioMaker->SetBranch("*",0,"0"); //deactivate all branches
150  ioMaker->SetBranch("eventBranch",0,"r"); //activate Event Branch
151  ioMaker->SetIOMode("r");
152 
153  St_db_Maker* dbMk = new St_db_Maker("db","MySQL:StarDb","$STAR/StarDb","StarDb");
154 
155  //StGenericVertexMaker *myfinder=new StGenericVertexMaker("myvertexfinder");
156 
157  StGenericVertexMaker *myfinder=new StGenericVertexMaker("myvertexfinder");
158  //StGenericVertexMaker *myfinder=new StMinuitVertexFinder("myvertexfinder");
159  //StMaker *myfinder=new StMinuitVertexFinder("myvertexfinder");
160  myfinder->SetDebug(1);
161  myfinder->SetMode(1); // Select Minuit finder
162  myfinder->SetInternalFind(); // Activate finding
163  myfinder->Init();
164 
165  chain->PrintInfo();
166  //chain->ls(3);
167  Int_t initStat = chain->Init(); // This should call the Init() method in ALL makers
168  if (initStat) chain->Fatal(initStat, "during Init()");
169 
170  int istat=0,iev=0;
171 
172  // Do the event loop
173  while(1) {
174  if (iev>=nevents) break;
175  chain->Clear();
176  cout << "---------------------- Processing Event : " << iev << " ---------------------- " << endl;
177  istat = chain->Make();
178  iev++;
179  if(istat) break;
180  cout << "istat " << istat<<endl;
181 
182  // if(iev<17) continue;
183  if (istat == kStEOF || istat == kStFatal) break;
184 
185  StEvent* mEvent = (StEvent*)chain->GetInputDS("StEvent");
186  assert(mEvent);// fix your chain or open the right event file
187 
188  //StMinuitVertexFinder* test_finder=(StMinuitVertexFinder*) myfinder->GetGenericFinder();
189  //myfinder->GetGenericFinder().printInfo();
190  //StMinuitVertexFinder test_finder;//=new StMinuitVertexFinder();
191  //test_finder->setPrintLevel(3);
192  //StThreeVectorD myvertex;
193  //if (myfinder->fit(mEvent)) {
194  //myvertex = myfinder->result();
196  //}
197  //else
198  // cout << "Error: vertex fit failed, no vertex." << endl;
199 
200  //delete myfinder;
201 
202  //myfinder->Finish();
203  Int_t nV=mEvent->numberOfPrimaryVertices();
204  if (nV == 0) continue;
205  int iv;
206  Float_t best_rank=1e9;
207  StPrimaryVertex *best_vtx=0;
208  for(iv=0;iv<nV;iv++) {
209  StPrimaryVertex *V=mEvent->primaryVertex(iv);
210  assert(V);
211  StThreeVectorF &r=V->position();
212  StThreeVectorF &er=V->positionError();
213  Float_t rank=eval_print_vertex(mEvent, V);
214  if (best_vtx==0 || rank < best_rank) {
215  best_vtx = V;
216  best_rank = rank;
217  }
218  }
219  if (best_vtx)
220  cout << "Best vertex: " << *best_vtx << endl;
221 // printf(" nPrimTr=%d , VFid=%d:: ntrVF=%d nCtb=%d nBemc=%d nEEmc=%d nTpc=%d sumPt=%.1f rank=%g\n"
222 // ,V->numberOfDaughters(), V->vertexFinderId() ,V->numTracksUsedInFinder() ,
223 // V->numMatchesWithCTB() ,V-> numMatchesWithBEMC() ,V->numMatchesWithEEMC() ,
224 // V->numTracksCrossingCentralMembrane() ,V->sumOfTrackPt() ,V->ranking());
225 
226 // continue;
227 // int nPrTr=0;
228 // //.... access prim tracks for given vertex
229 // int itr;
230 // for(itr=0; itr<V->numberOfDaughters(); itr++) {
231 // StTrack *track=V-> daughter(itr);
232 // if(track==0) continue;
233 // if (track->flag() <0 ) continue;
234 // printf("itr=%d pT=%.1f eta=%.2f nFitP=%d DCA=%.1f\n",itr,
235 // track->geometry()->momentum().mag(),
236 // track->geometry()->momentum().pseudoRapidity(),
237 // track->fitTraits().numberOfFitPoints(),
238 // track->geometry()->helix().distance(V->position()));
239 // nPrTr++;
240 // }
241 
242 // printf(" counted nPrimTr=%d \n",nPrTr);
243 // } // end of loop over vertices
244 
245 
246  } // Event Loop
247  //chain->Finish();
248  fout->Write();
249 }
250 
251 
252 //
253 // $Id: find_vertex.C,v 3.3 2021/03/20 02:38:13 genevb Exp $
254 // $Log: find_vertex.C,v $
255 // Revision 3.3 2021/03/20 02:38:13 genevb
256 // Remove StTpcDb dependence for StEventUtilities
257 //
258 // Revision 3.2 2021/03/19 14:49:29 genevb
259 // Add StTpcDb dependency for StEvenntUtilities
260 //
261 // Revision 3.1 2008/03/19 19:39:08 genevb
262 // Introduction of macro
263 //
264 //
St_db_Maker
Definition: St_db_Maker.h:141
StIOInterFace::SetIOMode
virtual void SetIOMode(Option_t *iomode="w")
number of transactions
Definition: StIOInterFace.h:35
StChain::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StChain.cxx:77
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
kStEOF
Definition: Stypes.h:43
StChain::Make
virtual Int_t Make()
Definition: StChain.cxx:110
Generated by   doxygen 1.8.5