General tracking macro by J. Chudoba.
[u/mrichter/AliRoot.git] / TPC / AliTPCTracking.C
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // AliTPCTracking.C 
4 //
5 // date: 22.08.2002
6 // author: Jiri Chudoba based on code of Jourij Belikov
7 // version: 1.0
8 // description: 
9 //      reconstructs of tracks in TPC inthe following steps:
10 //         TPC cluster finding
11 //         TPC track finding
12 // input parameters: 
13 //        Int_t nEvents      ... nr of events to process
14 //        Int_t firstEventNr ... first event number (starts from 0)
15 //        TString fileNameHits ... name of file with hits
16 //        TString fileNameDigits .. name of file with TPC digits
17 //        TString fileNameClusters .. name of file with TPC clusters (output)
18 //        TString fileNameTracks .. name of file with TPC tracks (output)
19 //
20 //        default file names correspond to pp production (2002-04)
21 //
22 // History:
23 //
24 //     22.08.2002 ... first version
25 //
26 ////////////////////////////////////////////////////////////////////////
27
28 #if !defined(__CINT__) || defined(__MAKECINT__)
29 #include "iostream.h"
30 #include "TTree.h"
31 #include "TSystem.h"
32 #include "TArrayF.h"
33 #include "TPC/alles.h"
34 #include "TPC/AliTPCtracker.h"
35 #include "TPC/AliTPCclusterer.h"
36 #include "STEER/AliRun.h"
37 #include "STEER/AliHeader.h"
38 #include "STEER/AliGenEventHeader.h"
39 #include "STEER/AliMagF.h"
40
41 #endif
42
43 Int_t gDEBUG = 2;
44
45 Int_t TPCFindClusters(Int_t nEvents=1, Int_t firstEvent=0,
46                       TString fileNameDigits="rfio:galiceSDR.root", 
47                       TString fileNameClusters="tpc.clusters.root");
48 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
49                       TFile* fileDigits, TFile* fileClusters, 
50                       AliTPCParam* paramTPC=0);
51 Int_t TPCFindTracks(Int_t nEvents=1, Int_t firstEvent=0,
52                     TString fileNameClusters="tpc.clusters.root",
53                     TString fileNameTracks="tpc.tracks.root");
54 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
55                     TFile* fileClusters, TFile* fileTracks,
56                     AliTPCParam* paramTPC=0);
57
58 AliTPCParam* LoadTPCParam(TFile *file);
59 void FindVertex(Int_t iEvent, Double_t *vertex);
60 void PrintVertex(TArrayF &primaryVertex);
61 Int_t SetFieldFactor(TString fileName, Bool_t closeFile = kTRUE);
62 Int_t SetFieldFactor(TFile* file, Bool_t deletegAlice = kTRUE);
63 Int_t SetFieldFactor();
64
65 Int_t AliTPCTracking(Int_t nEvents=1, Int_t firstEvent=0,
66                      TString fileNameHits="rfio:galice.root",
67                      TString fileNameDigits="rfio:galiceSDR.root",
68                      TString fileNameClusters="tpc.clusters.root",
69                      TString fileNameTracks="tpc.tracks.root");
70
71 ////////////////////////////////////////////////////////////////////////
72 Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent,
73                       TString fileNameHits,
74                       TString fileNameDigits,
75                       TString fileNameClusters,
76                       TString fileNameTracks) {
77
78   SetFieldFactor(fileNameHits,kFALSE);
79
80 // ********** Find TPC clusters *********** //
81   if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters)) {
82     cerr<<"Failed to get TPC clusters: !\n";
83     return 1;
84   }      
85
86 // ********** Find TPC tracks *********** //
87   if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks)) {
88     cerr<<"Failed to get TPC tracks !\n";
89     return 2;
90   }
91
92   return 0;
93 }
94
95 ////////////////////////////////////////////////////////////////////////
96 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
97                       TString fileNameDigits, TString fileNameClusters) {
98   
99   Int_t rc;
100   const Char_t *name="TPCFindClusters";
101   if (gDEBUG>1) cout<<name<<" starts...\n";
102   if (gDEBUG>1) gBenchmark->Start(name);
103   TFile *fileClusters = TFile::Open(fileNameClusters,"recreate");
104   TFile *fileDigits = TFile::Open(fileNameDigits);
105   if (!fileDigits->IsOpen()) {
106     cerr<<"Cannnot open "<<fileNameDigits<<" !\n"; 
107     return 1;
108   }
109   if (!fileClusters->IsOpen()) {
110     cerr<<"Cannnot open "<<fileNameClusters<<" !\n"; 
111     return 1;
112   }
113
114   rc = TPCFindClusters(nEvents,firstEvent,fileDigits,fileClusters);
115
116   fileDigits->Close();
117   fileClusters->Close();
118   delete fileDigits;
119   delete fileClusters;
120   if (gDEBUG>1) gBenchmark->Stop(name);
121   if (gDEBUG>1) gBenchmark->Show(name);
122
123   return rc;
124 }
125 ////////////////////////////////////////////////////////////////////////
126 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
127                       TFile* fileDigits, TFile* fileClusters, 
128                       AliTPCParam* paramTPC) {
129
130   fileDigits->cd();
131   if (!paramTPC) paramTPC = LoadTPCParam(fileDigits);
132   if (!paramTPC) return 1;
133
134   for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
135     if (gDEBUG > 2) cout<<"TPCFindClusters: event "<<iEvent<<endl;
136     AliTPCclusterer::Digits2Clusters(paramTPC, fileClusters, iEvent);
137   }
138   return 0;
139 }
140 ////////////////////////////////////////////////////////////////////////
141 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
142                     TString fileNameClusters, TString fileNameTracks) {
143
144   Int_t rc = 0;
145   const Char_t *name="TPCFindTracks";
146   if (gDEBUG>1) cout<<name<<" starts"<<endl;
147   if (gDEBUG>1) gBenchmark->Start(name);
148   TFile *fileTracks = TFile::Open(fileNameTracks,"recreate");
149   TFile *fileClusters =TFile::Open(fileNameClusters);
150
151   rc = TPCFindTracks(nEvents, firstEvent, fileClusters, fileTracks);
152
153   fileClusters->Close();
154   fileTracks->Close();
155   delete fileClusters;
156   delete fileTracks;
157   if (gDEBUG>1) gBenchmark->Stop(name);
158   if (gDEBUG>1) gBenchmark->Show(name);
159   return rc;
160
161 }
162 ////////////////////////////////////////////////////////////////////////
163 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
164                     TFile *fileClusters, TFile * fileTracks,
165                     AliTPCParam* paramTPC) {
166
167   Int_t rc = 0;
168   if (!paramTPC) paramTPC = LoadTPCParam(fileClusters);
169   if (!paramTPC) return 1;
170
171   for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
172     if (gDEBUG > 2) cout<<"TPCFindTracks: event "<<iEvent<<endl;
173     AliTPCtracker *tracker = new AliTPCtracker(paramTPC,iEvent);
174     Double_t vertex[3];
175     FindVertex(iEvent,vertex);
176     tracker->SetVertex(vertex);
177     fileClusters->cd();
178     rc = tracker->Clusters2Tracks(0,fileTracks);
179     delete tracker;
180   }
181   return rc;
182 }
183
184 ////////////////////////////////////////////////////////////////////////
185 Int_t SetFieldFactor() {
186
187    AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
188    if (gDEBUG > 2) cout<<"Magnetic field in kGauss: "<<gAlice->Field()->SolenoidField()<<endl;
189    return 0;
190 }
191 ////////////////////////////////////////////////////////////////////////
192 Int_t SetFieldFactor(TFile *file, Bool_t deletegAlice) {
193
194   if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
195     cerr<<"gAlice has not been found in file "<<file->GetName();
196     return 1;
197   }   
198   Int_t rc = SetFieldFactor();
199   if (deletegAlice) {
200     delete gAlice;  
201     gAlice = 0;
202   }
203   return rc;
204 }
205 ////////////////////////////////////////////////////////////////////////
206 Int_t SetFieldFactor(TString fileName, Bool_t closeFile) {
207
208    TFile *file=TFile::Open(fileName);
209    if (!file->IsOpen()) {cerr<<"Cannnot open "<<fileName<<" !\n"; return 1;}
210    Int_t rc = SetFieldFactor(file, closeFile) ;
211    if (closeFile) file->Close();
212    return rc;
213 }
214 ////////////////////////////////////////////////////////////////////////
215 AliTPCParam* LoadTPCParam(TFile *file) {
216
217   char paramName[50];
218   sprintf(paramName,"75x40_100x60_150x60");
219   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
220   if (paramTPC) {
221     if (gDEBUG > 1) cout<<"TPC parameters "<<paramName<<" found."<<endl;
222   } else {
223     cerr<<"TPC parameters not found. Create new (they may be incorrect)."
224         <<endl;    
225     paramTPC = new AliTPCParamSR;
226   }
227   return paramTPC;
228
229 // the older version of parameters can be accessed with this code.
230 // In some cases, we have old parameters saved in the file but 
231 // digits were created with new parameters, it can be distinguish 
232 // by the name of TPC TreeD. The code here is just for the case 
233 // we would need to compare with old data, uncomment it if needed.
234 //
235 //  char paramName[50];
236 //  sprintf(paramName,"75x40_100x60");
237 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
238 //  if (paramTPC) {
239 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
240 //  } else {
241 //    sprintf(paramName,"75x40_100x60_150x60");
242 //    paramTPC=(AliTPCParam*)in->Get(paramName);
243 //    if (paramTPC) {
244 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
245 //    } else {
246 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
247 //          <<endl;    
248 //      paramTPC = new AliTPCParamSR;
249 //    }
250 //  }
251 //  return paramTPC;
252
253 }
254
255 ////////////////////////////////////////////////////////////////////////
256 void FindVertex(Int_t eventNr, Double_t *vertex) {
257
258   vertex[0] = vertex[1] = vertex[2] = 0.;
259   if (!gAlice) {
260     cerr<<"gAlice was not found! Using vertex position (0,0,0).\n";
261     return;
262   }
263   
264   gAlice->GetEvent(eventNr);
265   AliHeader *header = gAlice->GetHeader();
266   if (!header) {
267     cerr<<"header was not found!\n";
268     return;
269   } 
270   AliGenEventHeader* genEventHeader = header->GenEventHeader();
271   if (!genEventHeader) {
272     cerr<<"AliGenEventHeader was not found!\n";
273     return;
274   } 
275
276   TArrayF primaryVertex(3);
277   genEventHeader->PrimaryVertex(primaryVertex);
278   PrintVertex(primaryVertex);
279   vertex[0] = static_cast<Double_t>(primaryVertex[0]);
280   vertex[1] = static_cast<Double_t>(primaryVertex[1]);
281   vertex[2] = static_cast<Double_t>(primaryVertex[2]);
282 //  delete header;
283   delete genEventHeader;
284   return;
285      
286 }
287 ////////////////////////////////////////////////////////////////////////
288 void PrintVertex(TArrayF &primaryVertex) 
289 {
290   cout <<"Vertex: "
291        <<primaryVertex[0]<<" "
292        <<primaryVertex[1]<<" "
293        <<primaryVertex[2]<<" "<<endl;
294   exit;
295
296 ////////////////////////////////////////////////////////////////////////