1 /****************************************************************************
2 * This macro performs track and vertex reconstruction in TPC and ITS. *
3 * The ITS Kalman tracker V2 is feeded "with" parameterized TPC tracks. *
5 * Reconstruction is performed in the following steps: *
6 * 1) TPC tracking parameterization *
7 * 2) ITS clusters: slow or fast *
8 * 3) Primary vertex reconstruction *
9 * - read from event header for Pb-Pb events *
10 * - determined using points in pixels for pp/pA events *
11 * 4) ITS track finding V2 *
12 * - in pp/pA, redetermine the position of primary vertex *
13 * using the reconstructed tracks *
14 * 5) Create a reference file with simulation info (p,PDG...) *
16 * If mode='A' all 5 steps are executed *
17 * If mode='B' only steps 4-5 are executed *
19 * Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it *
20 * (from AliTPCtest.C & AliITStestV2.C by I.Belikov) *
21 ****************************************************************************/
23 // structure for track references
33 //===== Functions definition =================================================
35 void CopyVtx(const Char_t *inName,const Char_t *outName);
37 void ITSFindClustersV2(Char_t SlowOrFast);
39 void ITSFindTracksV2(Int_t *skipEvt);
41 void ITSMakeRefFile(Int_t *skipEvt);
43 void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt);
45 void PrimaryVertex(const Char_t *outName,Char_t vtxMode);
47 void TPCParamTracks(Int_t coll,Double_t Bfield);
49 Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName);
51 void VtxFromHeader(const Char_t *outName,Bool_t smear);
53 void VtxFromTracks(const Char_t *outName);
55 void ZvtxFromSPD(const Char_t *outName);
57 //=============================================================================
59 // number of events to be processed
62 Double_t gBfieldValue;
64 void AliBarrelRec_TPCparam(Int_t n=-1,Char_t mode='A') {
66 //---------------------------------------------------------------------
72 // _Type_of_collision_ (needed for TPC tracking parameterization)
73 // Available choices: !!! ONLY B = 0.4 TESLA !!!
74 // collcode = 0 -> PbPb6000 (HIJING with b<2fm)
75 // collcode = 1 -> low multiplicity: pp or pA
78 // _ITS_clusters_reconstruction_
79 // Available choices: (from AliITStestV2.C)
80 // SlowOrFast = 's' slow points
81 // SlowOrFast = 'f' fast points
82 Char_t SlowOrFast = 'f';
84 // _Primary_vertex_for_ITS_tracking_
86 // Vtx4Tracking = 'H' from event Header
88 // Vtx4Tracking = 'S' from event header + Smearing
89 // (x=15,y=15,z=10) micron
91 // Vtx4Tracking = 'P' z from pixels, x,y in(0,0)
92 Char_t Vtx4Tracking = 'P';
93 // _Primary_vertex_for_analysis_ (AliESDVertex stored in tracks file)
95 // Vtx4Analysis = 'C' Copy the same used for tracking
97 // Vtx4Analysis = 'T' x,y,z from Tracks
98 Char_t Vtx4Analysis = 'T';
101 //---------------------------------------------------------------------
103 const Char_t *name=" AliBarrelRec_TPCparam";
104 printf("\n %s\n",name);
105 gBenchmark->Start(name);
107 if(n==-1) { // read number of events to be processed from file
108 TFile *f = new TFile("galice.root");
109 gAlice = (AliRun*)f->Get("gAlice");
110 n = gAlice->GetEventsPerRun();
115 printf(" All %d events in file will be processed\n",n);
120 // conversion constant for kalman tracks
121 AliKalmanTrack::SetConvConst(100/0.299792458/gBfieldValue);
123 // Selection of execution mode
126 // Build TPC tracks with parameterization
127 TPCParamTracks(collcode,gBfieldValue);
130 ITSFindClustersV2(SlowOrFast);
132 // Vertex for ITS tracking
133 PrimaryVertex("Vtx4Tracking.root",Vtx4Tracking);
138 printf(" ---> only tracking in ITS <---\n");
140 // Update list of events to be skipped
141 if(!UpdateEvtsToSkip("itstracking.log","evtsToSkip.dat")) return;
146 // Mark events that have to be skipped (if any)
147 Int_t *skipEvt = new Int_t[gNevents];
148 for(Int_t i=0; i<gNevents; i++) skipEvt[i] = 0;
149 if(!gSystem->AccessPathName("evtsToSkip.dat",kFileExists))
150 MarkEvtsToSkip("evtsToSkip.dat",skipEvt);
153 ITSFindTracksV2(skipEvt);
155 // Vertex for analysis
156 PrimaryVertex("AliITStracksV2.root",Vtx4Analysis);
158 // Create ITS tracks reference file
159 ITSMakeRefFile(skipEvt);
164 gBenchmark->Stop(name);
165 gBenchmark->Show(name);
169 //-----------------------------------------------------------------------------
170 void CopyVtx(const Char_t *inName,const Char_t *outName) {
172 // Open input and output files
173 TFile *inFile = new TFile(inName);
174 TFile *outFile = new TFile(outName,"update");
180 for(Int_t ev=0; ev<gNevents; ev++) {
181 sprintf(vname,"Vertex_%d",ev);
182 AliESDVertex *vertex = (AliESDVertex*)inFile->Get(vname);
183 if(!vertex) continue;
198 //-----------------------------------------------------------------------------
199 void ITSFindClustersV2(Char_t SlowOrFast) {
201 printf("\n------------------------------------\n");
203 const Char_t *name="ITSFindClustersV2";
204 printf("\n %s\n",name);
205 gBenchmark->Start(name);
207 //--- taken from AliITStestV2.C--------------------------------------
209 if (SlowOrFast=='f') {
210 //cerr<<"Fast AliITSRecPoint(s) !\n";
211 //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
212 //AliITSHits2FastRecPoints();
214 gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2SDigits.C");
215 AliITSHits2SDigits();
216 gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSSDigits2Digits.C");
217 AliITSSDigits2Digits();
218 //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSDigits2RecPoints.C");
219 //AliITSDigits2RecPoints();
221 gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersV2.C");
222 AliITSFindClustersV2(SlowOrFast,gNevents);
224 //--------------------------------------------------------------------
227 gBenchmark->Stop(name);
228 gBenchmark->Show(name);
232 //-----------------------------------------------------------------------------
233 Int_t ITSFindTracksV2(Int_t *skipEvt) {
235 printf("\n------------------------------------\n");
237 const Char_t *name="ITSFindTracksV2";
238 printf("\n %s\n",name);
239 gBenchmark->Start(name);
242 TFile *outFile = new TFile("AliITStracksV2.root","recreate");
243 TFile *inTPCtrks = new TFile("AliTPCtracksParam.root");
244 TFile *inVertex = new TFile("Vtx4Tracking.root");
245 TFile *inClusters = new TFile("AliITSclustersV2.root");
247 AliITSgeom *geom=(AliITSgeom*)inClusters->Get("AliITSgeom");
248 if(!geom) { printf("can't get ITS geometry !\n"); return;}
251 Int_t flag1stPass,flag2ndPass;
254 // open logfile for done events
255 FILE *logfile = fopen("itstracking.log","w");
257 // Instantiate AliITStrackerV2
258 AliITStrackerV2 tracker(geom);
261 for(Int_t ev=0; ev<gNevents; ev++){
262 // write to logfile of begun events
263 fprintf(logfile,"%d\n",ev);
265 if(skipEvt[ev]) continue;
266 printf(" --- Processing event %d ---\n",ev);
268 // pass event number to the tracker
269 tracker.SetEventNumber(ev);
271 // set position of primary vertex
272 sprintf(vname,"Vertex_%d",ev);
273 AliESDVertex *vertex = (AliESDVertex*)inVertex->Get(vname);
278 printf(" AliESDVertex not found for event %d\n",ev);
279 printf(" Using (0,0,0) for ITS tracking\n");
280 vtx[0] = vtx[1] = vtx[2] = 0.;
283 flag1stPass=1; // vtx constraint
284 flag2ndPass=0; // no vtx constraint
286 // no vtx constraint if vertex not found
292 tracker.SetVertex(vtx);
294 // setup vertex constraint in the two tracking passes
296 flags[0]=flag1stPass;
297 tracker.SetupFirstPass(flags);
298 flags[0]=flag2ndPass;
299 tracker.SetupSecondPass(flags);
302 tracker.Clusters2Tracks(inTPCtrks,outFile);
306 fprintf(logfile,"%d\n",gNevents); //this means all evts are successfully completed
317 gBenchmark->Stop(name);
318 gBenchmark->Show(name);
322 //-----------------------------------------------------------------------------
323 void ITSMakeRefFile(Int_t *skipEvt) {
325 printf("\n------------------------------------\n");
327 const Char_t *name="ITSMakeRefFile";
328 printf("\n %s\n",name);
329 gBenchmark->Start(name);
332 TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
333 TFile *trk = TFile::Open("AliITStracksV2.root");
334 TFile *kin = TFile::Open("galice.root");
337 // Get gAlice object from file
338 gAlice=(AliRun*)kin->Get("gAlice");
346 for(Int_t ev=0; ev<gNevents; ev++){
347 if(skipEvt[ev]) continue;
348 printf(" --- Processing event %d ---\n",ev);
350 gAlice->GetEvent(ev);
354 // Tree with ITS tracks
356 sprintf(tname,"TreeT_ITS_%d",ev);
358 TTree *tracktree=(TTree*)trk->Get(tname);
359 if(!tracktree) continue;
360 AliITStrackV2 *itstrack=new AliITStrackV2;
361 tracktree->SetBranchAddress("tracks",&itstrack);
362 Int_t nentr=(Int_t)tracktree->GetEntries();
364 // Tree for true track parameters
366 sprintf(ttname,"Tree_Ref_%d",ev);
367 TTree *reftree = new TTree(ttname,"Tree with true track params");
368 reftree->Branch("rectracks",&rectrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
370 for(Int_t i=0; i<nentr; i++) {
371 tracktree->GetEvent(i);
372 label = TMath::Abs(itstrack->GetLabel());
374 Part = (TParticle*)gAlice->Particle(label);
376 rectrk.pdg=Part->GetPdgCode();
377 rectrk.mumlab = Part->GetFirstMother();
378 if(Part->GetFirstMother()>=0) {
379 Mum = (TParticle*)gAlice->Particle(Part->GetFirstMother());
380 rectrk.mumpdg=Mum->GetPdgCode();
384 rectrk.Vx=Part->Vx();
385 rectrk.Vy=Part->Vy();
386 rectrk.Vz=Part->Vz();
387 rectrk.Px=Part->Px();
388 rectrk.Py=Part->Py();
389 rectrk.Pz=Part->Pz();
405 gBenchmark->Stop(name);
406 gBenchmark->Show(name);
411 //-----------------------------------------------------------------------------
412 void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt) {
414 printf("\n------------------------------------\n");
415 printf("\nChecking for events to skip...\n");
419 FILE *f = fopen(evtsName,"r");
421 ncol = fscanf(f,"%d",&evt);
424 printf(" event %d will be skipped\n",evt);
430 //-----------------------------------------------------------------------------
431 void PrimaryVertex(const Char_t *outName,Char_t vtxMode) {
433 printf("\n------------------------------------\n");
435 const Char_t *name="PrimaryVertex";
436 printf("\n %s\n",name);
437 gBenchmark->Start(name);
441 printf(" ... from event header\n");
442 VtxFromHeader(outName,kFALSE);
445 printf(" ... from event header + smearing\n");
446 VtxFromHeader(outName,kTRUE);
449 printf(" ... z from pixels for pp/pA\n");
450 ZvtxFromSPD(outName);
453 printf(" ... from tracks for pp/pA\n");
454 VtxFromTracks(outName);
457 printf(" ... copied from Vtx4Tracking.root to AliITStracksV2.root\n");
458 CopyVtx("Vtx4Tracking.root",outName);
462 gBenchmark->Stop(name);
463 gBenchmark->Show(name);
467 //-----------------------------------------------------------------------------
468 void TPCParamTracks(Int_t coll,Double_t Bfield) {
470 printf("\n------------------------------------\n");
472 const Char_t *name="TPCParamTracks";
473 printf("\n %s\n",name);
474 gBenchmark->Start(name);
476 TFile *outFile=TFile::Open("AliTPCtracksParam.root","recreate");
477 TFile *inFile =TFile::Open("galice.root");
479 AliTPCtrackerParam tracker(coll,Bfield,gNevents);
480 tracker.BuildTPCtracks(inFile,outFile);
482 delete gAlice; gAlice=0;
487 gBenchmark->Stop(name);
488 gBenchmark->Show(name);
492 //-----------------------------------------------------------------------------
493 Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName) {
495 if(!gSystem->AccessPathName(logName,kFileExists)) {
496 FILE *ifile = fopen(logName,"r");
499 nCol = fscanf(ifile,"%d",&lEvt);
503 printf(" All events already reconstructed\n");
506 FILE *ofile = fopen("evtsToSkip.dat","a");
507 fprintf(ofile,"%d\n",lEvt);
511 printf("File itstracking.log not found\n");
516 //-----------------------------------------------------------------------------
517 void VtxFromHeader(const Char_t *outName,Bool_t smear) {
520 UInt_t seed = t.Get();
521 gRandom->SetSeed(seed);
523 TFile *galice = new TFile("galice.root");
524 TFile *outFile = new TFile(outName,"update");
527 Double_t pos[3],sigma[3];
541 for(Int_t ev=0; ev<gNevents; ev++){
542 printf(" event %d\n",ev);
543 sprintf(vname,"Vertex_%d",ev);
546 AliHeader* header = 0;
547 TTree* treeE = (TTree*)gDirectory->Get("TE");
548 treeE->SetBranchAddress("Header",&header);
550 AliGenEventHeader* genHeader = header->GenEventHeader();
552 // get primary vertex position
553 genHeader->PrimaryVertex(o);
554 pos[0] = (Double_t)o[0];
555 pos[1] = (Double_t)o[1];
556 pos[2] = (Double_t)o[2];
558 pos[0] = gRandom->Gaus(pos[0],sigma[0]);
559 pos[1] = gRandom->Gaus(pos[1],sigma[1]);
560 pos[2] = gRandom->Gaus(pos[2],sigma[2]);
562 // create AliESDVertex
563 AliESDVertex *vertex = new AliESDVertex(pos,sigma,vname);
565 printf(" ! event header not found : setting vertex to (0,0,0) !");
569 // create AliESDVertex
570 AliESDVertex *vertex = new AliESDVertex(pos,sigma,vname);
573 // write AliESDVertex to file
577 vertex->SetTitle("vertex from header, smeared");
579 vertex->SetTitle("vertex from header");
594 //-----------------------------------------------------------------------------
595 void VtxFromTracks(const Char_t *outName) {
597 // Open input and output files
598 TFile *inFile = new TFile("AliITStracksV2.root");
599 TFile *outFile = new TFile(outName,"update");
601 // set AliRun object to 0
602 if(gAlice) gAlice = 0;
605 AliITSVertexerTracks *vertexer =
606 new AliITSVertexerTracks(inFile,outFile,gBfieldValue);
607 vertexer->SetFirstEvent(0);
608 vertexer->SetLastEvent(gNevents-1);
609 vertexer->SetDebug(0);
610 vertexer->PrintStatus();
612 vertexer->FindVertices();
623 //-----------------------------------------------------------------------------
624 void ZvtxFromSPD(const Char_t *outName) {
626 // create fast RecPoints, which are used for vertex finding
627 cerr<<"Fast AliITSRecPoint(s) !\n";
628 gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
629 AliITSHits2FastRecPoints(0,gNevents-1);
631 // delphi ---> azimuthal range to accept tracklets
632 // window ---> window in Z around the peak of tracklets proj. in mm
638 TFile *infile = new TFile("galice.root");
639 TFile *outfile = new TFile(outName,"update");
641 AliITSVertexerPPZ *vertexer = new AliITSVertexerPPZ(infile,outfile,initx,inity);
642 vertexer->SetFirstEvent(0);
643 vertexer->SetLastEvent(gNevents-1);
644 vertexer->SetDebug(0);
645 vertexer->SetDiffPhiMax(delphi);
646 vertexer->SetWindow(window);
647 vertexer->PrintStatus();
648 vertexer->FindVertices();
660 //-----------------------------------------------------------------------------