Changes due to the coding conventions
[u/mrichter/AliRoot.git] / TPC / AliBarrelRec_TPCparam.C
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.     * 
4  *                                                                          *
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...)   *
15  *                                                                          *
16  * If mode='A' all 5 steps are executed                                     *
17  * If mode='B' only steps 4-5 are executed                                  *
18  *                                                                          *  
19  *  Origin: A.Dainese, Padova,   andrea.dainese@pd.infn.it                  * 
20  *  (from AliTPCtest.C & AliITStestV2.C by I.Belikov)                       *
21  ****************************************************************************/
22
23 // structure for track references
24 typedef struct {
25   Int_t lab;
26   Int_t pdg;
27   Int_t mumlab;
28   Int_t mumpdg;
29   Float_t Vx,Vy,Vz;
30   Float_t Px,Py,Pz;
31 } RECTRACK;
32
33 //===== Functions definition ================================================= 
34
35 void CopyVtx(const Char_t *inName,const Char_t *outName);
36
37 void ITSFindClustersV2(Char_t SlowOrFast);
38
39 void ITSFindTracksV2(Int_t *skipEvt);
40
41 void ITSMakeRefFile(Int_t *skipEvt);
42
43 void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt);
44
45 void PrimaryVertex(const Char_t *outName,Char_t vtxMode);
46
47 void TPCParamTracks(Int_t coll,Double_t Bfield);
48
49 Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName);
50
51 void VtxFromHeader(const Char_t *outName,Bool_t smear);
52
53 void VtxFromTracks(const Char_t *outName);
54
55 void ZvtxFromSPD(const Char_t *outName);
56
57 //=============================================================================
58
59 // number of events to be processed
60 Int_t    gNevents;
61 // magnetic field
62 Double_t gBfieldValue;
63
64 void AliBarrelRec_TPCparam(Int_t n=-1,Char_t mode='A') {
65
66   //---------------------------------------------------------------------
67   //                    CONFIGURATION
68   //
69   // _Magnetic_field_
70   gBfieldValue = 0.4;
71   //
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
76   Int_t collcode = 1;  
77   // 
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';
83   //
84   // _Primary_vertex_for_ITS_tracking_
85   // Available choices:
86   //    Vtx4Tracking = 'H'   from event Header
87   //    --- for Pb-Pb ---
88   //    Vtx4Tracking = 'S'   from event header + Smearing 
89   //                                           (x=15,y=15,z=10) micron
90   //    --- for pp/pA ---
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)
94   // Available choices:
95   //    Vtx4Analysis = 'C'   Copy the same used for tracking
96   //    --- for pp/pA ---
97   //    Vtx4Analysis = 'T'   x,y,z from Tracks
98   Char_t Vtx4Analysis = 'T';
99   //
100   //                  END CONFIGURATION
101   //---------------------------------------------------------------------
102
103   const Char_t *name=" AliBarrelRec_TPCparam";
104   printf("\n %s\n",name);
105   gBenchmark->Start(name);
106
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();
111     delete gAlice; 
112     gAlice=0;
113     f->Close(); 
114     delete f;
115     printf(" All %d events in file will be processed\n",n);
116   }
117   gNevents = n;
118
119
120   // conversion constant for kalman tracks 
121   AliKalmanTrack::SetConvConst(100/0.299792458/gBfieldValue);
122
123   // Selection of execution mode
124   switch(mode) {
125   case 'A':
126     // Build TPC tracks with parameterization
127     TPCParamTracks(collcode,gBfieldValue);
128   
129     // ITS clusters
130     ITSFindClustersV2(SlowOrFast);
131
132     // Vertex for ITS tracking
133     PrimaryVertex("Vtx4Tracking.root",Vtx4Tracking);
134
135     break;
136   
137   case 'B':
138     printf("       ---> only tracking in ITS <---\n");
139
140     // Update list of events to be skipped
141     if(!UpdateEvtsToSkip("itstracking.log","evtsToSkip.dat")) return;
142
143     break;
144   }
145
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);
151     
152   // Tracking in ITS
153   ITSFindTracksV2(skipEvt); 
154
155   // Vertex for analysis 
156   PrimaryVertex("AliITStracksV2.root",Vtx4Analysis);
157
158   // Create ITS tracks reference file
159   ITSMakeRefFile(skipEvt);
160   delete [] skipEvt;
161
162
163
164   gBenchmark->Stop(name);
165   gBenchmark->Show(name);
166
167   return;
168 }
169 //-----------------------------------------------------------------------------
170 void CopyVtx(const Char_t *inName,const Char_t *outName) {
171
172   // Open input and output files
173   TFile *inFile = new TFile(inName);
174   TFile *outFile = new TFile(outName,"update");
175
176   TDirectory *curdir;
177   Char_t vname[20];
178
179
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;
184     curdir = gDirectory;
185     outFile->cd();
186     vertex->Write();
187     curdir->cd();
188     vertex = 0;
189   }
190
191   inFile->Close();
192   outFile->Close();
193   delete inFile;
194   delete outFile;
195
196   return;
197 }
198 //-----------------------------------------------------------------------------
199 void ITSFindClustersV2(Char_t SlowOrFast) {
200
201   printf("\n------------------------------------\n");
202
203   const Char_t *name="ITSFindClustersV2";
204   printf("\n %s\n",name);
205   gBenchmark->Start(name);
206
207   //---  taken from AliITStestV2.C--------------------------------------
208   //
209   if (SlowOrFast=='f') {
210     //cerr<<"Fast AliITSRecPoint(s) !\n";
211     //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
212     //AliITSHits2FastRecPoints();
213   } else {
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();
220   }
221   gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersV2.C");
222   AliITSFindClustersV2(SlowOrFast,gNevents);
223   //
224   //--------------------------------------------------------------------
225
226
227   gBenchmark->Stop(name);
228   gBenchmark->Show(name);
229
230    return;
231 }
232 //-----------------------------------------------------------------------------
233 Int_t ITSFindTracksV2(Int_t *skipEvt) {
234
235   printf("\n------------------------------------\n");
236   
237   const Char_t *name="ITSFindTracksV2";
238   printf("\n %s\n",name);
239   gBenchmark->Start(name);
240
241  
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");
246   
247   AliITSgeom *geom=(AliITSgeom*)inClusters->Get("AliITSgeom");
248   if(!geom) { printf("can't get ITS geometry !\n"); return;}
249   
250   Double_t vtx[3];
251   Int_t flag1stPass,flag2ndPass;
252   Char_t vname[20];
253
254   // open logfile for done events
255   FILE *logfile = fopen("itstracking.log","w");
256
257   // Instantiate AliITStrackerV2
258   AliITStrackerV2 tracker(geom);
259
260   // loop on events
261   for(Int_t ev=0; ev<gNevents; ev++){
262     // write to logfile of begun events
263     fprintf(logfile,"%d\n",ev);
264
265     if(skipEvt[ev]) continue;
266     printf(" --- Processing event %d ---\n",ev);
267
268     // pass event number to the tracker
269     tracker.SetEventNumber(ev);
270
271     // set position of primary vertex
272     sprintf(vname,"Vertex_%d",ev);
273     AliESDVertex *vertex = (AliESDVertex*)inVertex->Get(vname);
274     if(vertex) {
275       vertex->GetXYZ(vtx);
276       delete vertex;
277     } else {
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.;
281     }
282
283     flag1stPass=1; // vtx constraint
284     flag2ndPass=0; // no vtx constraint
285
286     // no vtx constraint if vertex not found
287     if(vtx[2]<-999.) {
288       flag1stPass=0;
289       vtx[2]=0.;
290     }
291
292     tracker.SetVertex(vtx);  
293
294     // setup vertex constraint in the two tracking passes
295     Int_t flags[2];
296     flags[0]=flag1stPass;
297     tracker.SetupFirstPass(flags);
298     flags[0]=flag2ndPass;
299     tracker.SetupSecondPass(flags);
300     
301     // find the tracks
302     tracker.Clusters2Tracks(inTPCtrks,outFile);
303
304   } // loop on events
305
306   fprintf(logfile,"%d\n",gNevents); //this means all evts are successfully completed
307   fclose(logfile);
308
309   delete geom;
310
311   inTPCtrks->Close();
312   inClusters->Close();
313   inVertex->Close();
314   outFile->Close();
315  
316
317   gBenchmark->Stop(name);
318   gBenchmark->Show(name);
319   
320   return;
321 }
322 //-----------------------------------------------------------------------------
323 void ITSMakeRefFile(Int_t *skipEvt) {
324
325   printf("\n------------------------------------\n");
326
327   const Char_t *name="ITSMakeRefFile";
328   printf("\n %s\n",name);
329   gBenchmark->Start(name);
330   
331   
332   TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
333   TFile *trk = TFile::Open("AliITStracksV2.root");
334   TFile *kin = TFile::Open("galice.root");
335
336   
337   // Get gAlice object from file
338   gAlice=(AliRun*)kin->Get("gAlice");
339   
340   Int_t label;
341   TParticle *Part;  
342   TParticle *Mum;
343   RECTRACK rectrk;
344   
345
346   for(Int_t ev=0; ev<gNevents; ev++){
347     if(skipEvt[ev]) continue;
348     printf(" --- Processing event %d ---\n",ev);
349
350     gAlice->GetEvent(ev);  
351
352     trk->cd();
353
354     // Tree with ITS tracks
355     char tname[100];
356     sprintf(tname,"TreeT_ITS_%d",ev);
357
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();
363
364     // Tree for true track parameters
365     char ttname[100];
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");
369
370     for(Int_t i=0; i<nentr; i++) {
371       tracktree->GetEvent(i);
372       label = TMath::Abs(itstrack->GetLabel());
373
374       Part = (TParticle*)gAlice->Particle(label);
375       rectrk.lab=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();
381       } else {
382         rectrk.mumpdg=-1;
383       }
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();
390       
391       reftree->Fill();
392     } // loop on tracks   
393
394     out->cd();
395     reftree->Write();
396
397     delete itstrack;
398     delete reftree;
399   } // loop on events
400
401   trk->Close();
402   kin->Close();
403   out->Close();
404   
405   gBenchmark->Stop(name);
406   gBenchmark->Show(name);
407   
408
409   return;
410 }
411 //-----------------------------------------------------------------------------
412 void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt) {
413
414   printf("\n------------------------------------\n");
415   printf("\nChecking for events to skip...\n");
416
417   Int_t evt,ncol;
418
419   FILE *f = fopen(evtsName,"r");
420   while(1) {
421     ncol = fscanf(f,"%d",&evt);
422     if(ncol<1) break;
423     skipEvt[evt] = 1;
424     printf(" event %d will be skipped\n",evt);
425   }
426   fclose(f);
427
428   return;
429 }
430 //-----------------------------------------------------------------------------
431 void PrimaryVertex(const Char_t *outName,Char_t vtxMode) {
432
433   printf("\n------------------------------------\n");
434
435   const Char_t *name="PrimaryVertex";
436   printf("\n %s\n",name);
437   gBenchmark->Start(name);
438
439   switch(vtxMode) {
440   case 'H':
441     printf(" ... from event header\n");
442     VtxFromHeader(outName,kFALSE);
443     break;
444   case 'S':
445     printf(" ... from event header + smearing\n");
446     VtxFromHeader(outName,kTRUE);
447     break;
448   case 'P':
449     printf(" ... z from pixels for pp/pA\n");
450     ZvtxFromSPD(outName);
451     break;
452   case 'T':
453     printf(" ... from tracks for pp/pA\n");
454     VtxFromTracks(outName);
455     break;
456   case 'C':
457     printf(" ... copied from Vtx4Tracking.root to AliITStracksV2.root\n");
458     CopyVtx("Vtx4Tracking.root",outName);
459     break;
460   }
461
462   gBenchmark->Stop(name);
463   gBenchmark->Show(name);
464
465   return;
466 }
467 //-----------------------------------------------------------------------------
468 void TPCParamTracks(Int_t coll,Double_t Bfield) {
469
470   printf("\n------------------------------------\n");
471
472   const Char_t *name="TPCParamTracks";
473   printf("\n %s\n",name);
474   gBenchmark->Start(name);
475
476   TFile *outFile=TFile::Open("AliTPCtracksParam.root","recreate");
477   TFile *inFile =TFile::Open("galice.root");
478  
479   AliTPCtrackerParam tracker(coll,Bfield,gNevents);
480   tracker.BuildTPCtracks(inFile,outFile);
481
482   delete gAlice; gAlice=0;
483
484   inFile->Close();
485   outFile->Close();
486
487   gBenchmark->Stop(name);
488   gBenchmark->Show(name);
489
490   return;
491 }
492 //-----------------------------------------------------------------------------
493 Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName) {
494
495     if(!gSystem->AccessPathName(logName,kFileExists)) { 
496       FILE *ifile = fopen(logName,"r");
497       Int_t lEvt=0,nCol=1;
498       while(nCol>0) {
499         nCol = fscanf(ifile,"%d",&lEvt);
500       }
501       fclose(ifile);
502       if(lEvt==gNevents) { 
503         printf(" All events already reconstructed\n"); 
504         return 0;  
505       } else {
506         FILE *ofile = fopen("evtsToSkip.dat","a");
507         fprintf(ofile,"%d\n",lEvt);
508         fclose(ofile);
509       }
510     } else { 
511       printf("File itstracking.log not found\n");
512     }
513
514     return 1;
515 }
516 //-----------------------------------------------------------------------------
517 void VtxFromHeader(const Char_t *outName,Bool_t smear) {
518
519   TDatime t;
520   UInt_t seed = t.Get();
521   gRandom->SetSeed(seed);
522
523   TFile *galice  = new TFile("galice.root");  
524   TFile *outFile = new TFile(outName,"update");
525
526   TDirectory *curdir;
527   Double_t pos[3],sigma[3];
528   if(smear) {
529     sigma[0]=15.e-4;
530     sigma[1]=15.e-4;
531     sigma[2]=10.e-4;
532   } else {
533     sigma[0]=0.;
534     sigma[1]=0.;
535     sigma[2]=0.;
536   }
537   Char_t vname[20];
538
539   galice->cd();
540
541   for(Int_t ev=0; ev<gNevents; ev++){
542     printf(" event %d\n",ev);
543     sprintf(vname,"Vertex_%d",ev);
544     TArrayF o = 0;
545     o.Set(3);
546     AliHeader* header = 0;
547     TTree* treeE = (TTree*)gDirectory->Get("TE");
548     treeE->SetBranchAddress("Header",&header);
549     treeE->GetEntry(ev);
550     AliGenEventHeader* genHeader = header->GenEventHeader();
551     if(genHeader) {
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];
557       if(smear) {
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]);
561       }
562       // create AliESDVertex
563       AliESDVertex *vertex = new AliESDVertex(pos,sigma,vname);
564     } else {
565       printf(" ! event header not found : setting vertex to (0,0,0) !");
566       pos[0] = 0.;
567       pos[1] = 0.;
568       pos[2] = 0.;
569       // create AliESDVertex
570       AliESDVertex *vertex = new AliESDVertex(pos,sigma,vname);
571     }    
572     delete header;
573     // write AliESDVertex to file
574     curdir = gDirectory;
575     outFile->cd();
576     if(smear) {
577       vertex->SetTitle("vertex from header, smeared");
578     } else {
579       vertex->SetTitle("vertex from header");
580     }
581     vertex->Write();
582     curdir->cd();
583     vertex = 0;
584   }
585
586   outFile->Close();
587   galice->Close();
588
589   delete outFile;
590   delete galice;
591
592   return;
593 }
594 //-----------------------------------------------------------------------------
595 void VtxFromTracks(const Char_t *outName) {
596
597   // Open input and output files
598   TFile *inFile  = new TFile("AliITStracksV2.root");
599   TFile *outFile = new TFile(outName,"update");
600
601   // set AliRun object to 0
602   if(gAlice) gAlice = 0;
603
604   // Create vertexer
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();
611   // Find vertices
612   vertexer->FindVertices();
613
614   delete vertexer;
615
616   inFile->Close();
617   outFile->Close();
618   delete inFile;
619   delete outFile;
620
621   return;
622 }
623 //-----------------------------------------------------------------------------
624 void ZvtxFromSPD(const Char_t *outName) {
625
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);
630
631   // delphi ---> azimuthal range to accept tracklets
632   // window ---> window in Z around the peak of tracklets proj. in mm
633   Float_t delphi=0.05;
634   Float_t window=3.;
635   Float_t initx=0.;
636   Float_t inity=0.;
637
638   TFile *infile = new TFile("galice.root");
639   TFile *outfile = new TFile(outName,"update");
640
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();
649   delete vertexer;
650   vertexer=0;
651
652   outfile->Close();
653   infile->Close();
654   delete infile;
655   delete outfile;
656
657
658   return;
659 }
660 //-----------------------------------------------------------------------------
661
662
663
664
665
666
667
668
669