8683b0355fd4dc731a7fa26937ff344c680834d1
[u/mrichter/AliRoot.git] / STEER / AliESDtest.C
1 //********************************************************************
2 //     Example of the reconstruction that generates the ESD
3 // Input files: 
4 //   a) file containing the ITS clusters
5 //      (the AliITSFindClustersV2.C macro can be used to generate it)
6 //   b) file containing the TPC clusters
7 //      (the AliTPCFindClusters.C macro can be used to generate it)
8 //   c) file containing the TRD clusters
9 //      (the AliTRDdigits2cluster.C macro can be used to generate it)
10 //   d) file containing the TOF digits
11 //      (the AliTOFSDigits2Digits.C macro can be used to generate it)
12 // Ouput file:
13 //      AliESDs.root containing the ESD events 
14 //
15 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
16 //********************************************************************
17
18 #if !defined(__CINT__) || defined(__MAKECINT__)
19   #include <Riostream.h>
20   #include "TFile.h"
21   #include "TSystem.h"
22   #include "TStopwatch.h"
23   #include "TGeant3.h"
24
25   #include "AliMagF.h"
26   #include "AliRun.h"
27   #include "AliRunLoader.h"
28   #include "AliLoader.h"
29
30   #include "AliESD.h"
31   #include "AliESDpid.h"
32
33   #include "AliITS.h"
34   #include "AliITSgeom.h"
35   #include "AliITStrackerV2.h"
36   #include "AliV0vertexer.h"
37   #include "AliCascadeVertexer.h"
38   #include "AliITSpidESD.h"
39   #include "AliITSLoader.h"
40
41   #include "AliTPCParam.h"
42   #include "AliTPCtracker.h"
43   #include "AliTPCpidESD.h"
44   #include "AliTPCLoader.h"
45
46   #include "AliTRDtracker.h"
47   #include "AliTRDPartID.h"
48
49   #include "AliTOFpidESD.h"
50 #endif
51
52 extern TSystem *gSystem;
53 extern AliRun *gAlice;
54 extern TFile *gFile;
55
56 Int_t AliESDtest(Int_t nev=1) {
57
58 /**** Initialization of the NewIO *******/
59
60    if (gAlice) {
61       delete gAlice->GetRunLoader();
62       delete gAlice; 
63       gAlice=0;
64    }
65
66    gSystem->Load("libgeant321");     // needed for the PID in TOF 
67    new TGeant3("");                  // must be re-done !
68
69    AliRunLoader *rl = AliRunLoader::Open("galice.root");
70    if (rl == 0x0) {
71       cerr<<"Can not open session"<<endl;
72       return 1;
73    }
74    Int_t retval = rl->LoadgAlice();
75    if (retval) {
76       cerr<<"AliESDtest.C : LoadgAlice returned error"<<endl;
77       delete rl;
78       return 1;
79    }
80    retval = rl->LoadHeader();
81    if (retval) {
82       cerr<<"AliESDtest.C : LoadHeader returned error"<<endl;
83       delete rl;
84       return 2;
85    }
86    gAlice=rl->GetAliRun();
87        
88
89    AliKalmanTrack::SetConvConst(
90       1000/0.299792458/gAlice->Field()->SolenoidField()
91    );
92
93
94 /**** The ITS corner ********************/
95
96    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
97    if (itsl == 0x0) {
98       cerr<<"AliESDtest.C : Can not get the ITS loader"<<endl;
99       return 3;
100    }
101    itsl->LoadRecPoints("read");
102
103    AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
104    if (!dITS) {
105       cerr<<"AliESDtest.C : Can not find the ITS detector !"<<endl;
106       return 4;
107    }
108    AliITSgeom *geom = dITS->GetITSgeom();
109
110    //An instance of the ITS tracker
111    AliITStrackerV2 itsTracker(geom);
112    
113    //An instance of the ITS PID maker
114    Double_t parITS[]={34.,0.15,10.};
115    AliITSpidESD itsPID(parITS);
116
117    //An instance of the V0 finder
118    Double_t cuts[]={33,  // max. allowed chi2
119                     0.16,// min. allowed negative daughter's impact parameter 
120                     0.05,// min. allowed positive daughter's impact parameter 
121                     0.080,// max. allowed DCA between the daughter tracks
122                     0.998,// max. allowed cosine of V0's pointing angle
123                     0.9,  // min. radius of the fiducial volume
124                     2.9   // max. radius of the fiducial volume
125                    };
126    AliV0vertexer vtxer(cuts);
127
128    Double_t cts[]={33.,    // max. allowed chi2
129                     0.05,   // min. allowed V0 impact parameter 
130                     0.008,  // window around the Lambda mass 
131                     0.035,  // min. allowed bachelor's impact parameter 
132                     0.10,   // max. allowed DCA between a V0 and a track
133                     0.9985, //max. allowed cosine of the cascade pointing angle
134                     0.9,    // min. radius of the fiducial volume
135                     2.9     // max. radius of the fiducial volume
136                    };
137    AliCascadeVertexer cvtxer=AliCascadeVertexer(cts);
138
139 /**** The TPC corner ********************/
140
141    AliTPCLoader* tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
142    if (tpcl == 0x0) {
143       cerr<<"AliESDtest.C : can not get the TPC loader"<<endl;
144       return 5;
145    }
146    tpcl->LoadRecPoints("read");
147
148    rl->CdGAFile();
149    AliTPCParam *par=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
150    if (!par) { 
151       cerr<<"TPC parameters have not been found !\n";
152       return 6;
153    }
154    
155    //An instance of the TPC tracker
156    AliTPCtracker tpcTracker(par);
157
158    //An instance of the TPC PID maker
159    Double_t parTPC[]={47.,0.10,10.};
160    AliTPCpidESD tpcPID(parTPC);
161
162
163 /**** The TRD corner ********************/
164
165    AliLoader* trdl = rl->GetLoader("TRDLoader");
166    if (trdl == 0x0) {
167       cerr<<"AliESDtest.C : can not get the TRD loader"<<endl;
168       return 5;
169    }
170    trdl->LoadRecPoints("read");
171
172    //An instance of the TRD tracker
173    rl->CdGAFile();
174    AliTRDtracker trdTracker(gFile);  //galice.root file
175
176 /*
177    //An instance of the TRD PID maker
178    TFile* pidFile = TFile::Open("pid.root");
179    if (!pidFile->IsOpen()) {
180      cerr << "Can't get pid.root !\n";
181      return 7;
182    }
183    AliTRDPartID* trdPID = (AliTRDPartID*) pidFile->Get("AliTRDPartID");
184    if (!trdPID) {
185      cerr << "Can't get PID object !\n";
186      return 8;
187    }
188 */
189
190
191 /**** The TOF corner ********************/
192
193    AliLoader* tofl = rl->GetLoader("TOFLoader");
194    if (tofl == 0x0) {
195       cerr<<"AliESDtest.C : can not get the TOF loader"<<endl;
196       return 5;
197    }
198    tofl->LoadDigits("read");
199
200
201    //Instance of the TOF PID maker
202    Double_t parTOF[]={130.,5.};
203    AliTOFpidESD tofPID(parTOF);
204
205
206    //rl->UnloadgAlice();
207
208
209    TFile *ef=TFile::Open("AliESDs.root","RECREATE");
210    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
211
212    TStopwatch timer;
213    Int_t rc=0;
214    if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents();
215    //The loop over events
216    for (Int_t i=0; i<nev; i++) {
217      cerr<<"\n\nProcessing event number : "<<i<<endl;
218      AliESD *event=new AliESD(); 
219
220      rl->GetEvent(i);
221  
222
223 //***** Initial path towards the primary vertex
224      TTree *tpcTree=tpcl->TreeR();
225      if (!tpcTree) {
226         cerr<<"Can't get the TPC cluster tree !\n";
227         return 4;
228      }     
229      tpcTracker.LoadClusters(tpcTree);
230      rc+=tpcTracker.Clusters2Tracks(event);
231
232      TTree *itsTree=itsl->TreeR();
233      if (!itsTree) {
234         cerr<<"Can't get the ITS cluster tree !\n";
235         return 4;
236      }     
237      itsTracker.LoadClusters(itsTree);
238      rc+=itsTracker.Clusters2Tracks(event);
239
240
241 //***** Back propagation towards the outer barrel detectors
242      rc+=itsTracker.PropagateBack(event); 
243      //itsPID.MakePID(event);
244      
245      rc+=tpcTracker.PropagateBack(event);
246      tpcPID.MakePID(event);
247
248      TTree *trdTree=trdl->TreeR();
249      if (!trdTree) {
250         cerr<<"Can't get the TRD cluster tree !\n";
251         return 4;
252      } 
253      trdTracker.LoadClusters(trdTree);
254      rc+=trdTracker.PropagateBack(event);
255 /*
256      for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
257        AliESDtrack* track = event->GetTrack(iTrack);
258        trdPID->MakePID(track);
259      }
260 */
261      TTree *tofTree=tofl->TreeD();
262      if (!tofTree) {
263         cerr<<"Can't get the TOF cluster tree !\n";
264         return 4;
265      } 
266      tofPID.LoadClusters(tofTree);
267      tofPID.MakePID(event);
268      tofPID.UnloadClusters();
269
270
271
272 //***** Here is the combined PID
273      AliESDpid::MakePID(event);
274
275
276
277 //***** Now the final refit at the primary vertex...
278      rc+=trdTracker.RefitInward(event);
279      trdTracker.UnloadClusters();
280
281      rc+=tpcTracker.RefitInward(event);
282      tpcTracker.UnloadClusters();
283
284      rc+=itsTracker.RefitInward(event); 
285      itsTracker.UnloadClusters();
286
287
288
289 //***** Hyperon reconstruction 
290      rc+=vtxer.Tracks2V0vertices(event);            // V0 finding
291      rc+=cvtxer.V0sTracks2CascadeVertices(event);   // cascade finding
292
293
294
295 //***** Some final manipulations with this event 
296      if (rc==0) {
297         Char_t ename[100]; 
298         sprintf(ename,"%d",i);
299         ef->cd();
300         if (!event->Write(ename)) rc++;
301      } 
302      if (rc) {
303         cerr<<"Something bad happened...\n";
304      }
305      delete event;
306    }
307    timer.Stop(); timer.Print();
308
309    //pidFile->Close();
310
311    delete par;
312
313    ef->Close();
314
315    delete rl;
316
317    return rc;
318 }