3ddf801c487e07a7e03278b7d410e2cd52385069
[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
95 /**** The ITS corner ********************/
96
97    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
98    if (itsl == 0x0) {
99       cerr<<"AliESDtest.C : Can not get the ITS loader"<<endl;
100       return 3;
101    }
102    itsl->LoadRecPoints("read");
103
104    AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
105    if (!dITS) {
106       cerr<<"AliESDtest.C : Can not find the ITS detector !"<<endl;
107       return 4;
108    }
109    AliITSgeom *geom = dITS->GetITSgeom();
110
111    //An instance of the ITS tracker
112    AliITStrackerV2 itsTracker(geom);
113    
114    //An instance of the ITS PID maker
115    Double_t parITS[]={34.,0.15,10.};
116    AliITSpidESD itsPID(parITS);
117
118    //An instance of the V0 finder
119    Double_t cuts[]={33,  // max. allowed chi2
120                     0.16,// min. allowed negative daughter's impact parameter 
121                     0.05,// min. allowed positive daughter's impact parameter 
122                     0.080,// max. allowed DCA between the daughter tracks
123                     0.998,// max. allowed cosine of V0's pointing angle
124                     0.9,  // min. radius of the fiducial volume
125                     2.9   // max. radius of the fiducial volume
126                    };
127    AliV0vertexer vtxer(cuts);
128
129    Double_t cts[]={33.,    // max. allowed chi2
130                     0.05,   // min. allowed V0 impact parameter 
131                     0.008,  // window around the Lambda mass 
132                     0.035,  // min. allowed bachelor's impact parameter 
133                     0.10,   // max. allowed DCA between a V0 and a track
134                     0.9985, //max. allowed cosine of the cascade pointing angle
135                     0.9,    // min. radius of the fiducial volume
136                     2.9     // max. radius of the fiducial volume
137                    };
138    AliCascadeVertexer cvtxer=AliCascadeVertexer(cts);
139
140 /**** The TPC corner ********************/
141
142    AliTPCLoader* tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
143    if (tpcl == 0x0) {
144       cerr<<"AliESDtest.C : can not get the TPC loader"<<endl;
145       return 5;
146    }
147    tpcl->LoadRecPoints("read");
148
149    rl->CdGAFile();
150    AliTPCParam *par=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
151    if (!par) { 
152       cerr<<"TPC parameters have not been found !\n";
153       return 6;
154    }
155    
156    //An instance of the TPC tracker
157    AliTPCtracker tpcTracker(par);
158
159    //An instance of the TPC PID maker
160    Double_t parTPC[]={47.,0.10,10.};
161    AliTPCpidESD tpcPID(parTPC);
162
163
164 /**** The TRD corner ********************/
165
166    AliLoader* trdl = rl->GetLoader("TRDLoader");
167    if (trdl == 0x0) {
168       cerr<<"AliESDtest.C : can not get the TRD loader"<<endl;
169       return 5;
170    }
171    trdl->LoadRecPoints("read");
172
173    //An instance of the TRD tracker
174    rl->CdGAFile();
175    AliTRDtracker trdTracker(gFile);  //galice.root file
176
177 /*
178    //An instance of the TRD PID maker
179    TFile* pidFile = TFile::Open("pid.root");
180    if (!pidFile->IsOpen()) {
181      cerr << "Can't get pid.root !\n";
182      return 7;
183    }
184    AliTRDPartID* trdPID = (AliTRDPartID*) pidFile->Get("AliTRDPartID");
185    if (!trdPID) {
186      cerr << "Can't get PID object !\n";
187      return 8;
188    }
189 */
190
191
192 /**** The TOF corner ********************/
193
194    AliLoader* tofl = rl->GetLoader("TOFLoader");
195    if (tofl == 0x0) {
196       cerr<<"AliESDtest.C : can not get the TOF loader"<<endl;
197       return 5;
198    }
199    tofl->LoadDigits("read");
200
201
202    //Instance of the TOF PID maker
203    Double_t parTOF[]={130.,5.};
204    AliTOFpidESD tofPID(parTOF);
205
206
207    //rl->UnloadgAlice();
208
209
210    TFile *ef=TFile::Open("AliESDs.root","RECREATE");
211    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
212
213    TStopwatch timer;
214    Int_t rc=0;
215    if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents();
216    //The loop over events
217    for (Int_t i=0; i<nev; i++) {
218      cerr<<"\n\nProcessing event number : "<<i<<endl;
219      AliESD *event=new AliESD(); 
220
221      rl->GetEvent(i);
222  
223      TTree *tpcTree=tpcl->TreeR();
224      if (!tpcTree) {
225         cerr<<"Can't get the TPC cluster tree !\n";
226         return 4;
227      }     
228      tpcTracker.LoadClusters(tpcTree);
229      rc+=tpcTracker.Clusters2Tracks(event);
230
231
232      TTree *itsTree=itsl->TreeR();
233      if (!itsTree) {
234         cerr<<"Can't get the TPC cluster tree !\n";
235         return 4;
236      }     
237      itsTracker.LoadClusters(itsTree);
238      rc+=itsTracker.Clusters2Tracks(event);
239
240      rc+=vtxer.Tracks2V0vertices(event);            // V0 finding
241      rc+=cvtxer.V0sTracks2CascadeVertices(event);   // cascade finding
242
243      rc+=itsTracker.PropagateBack(event); 
244      itsTracker.UnloadClusters();
245      //itsPID.MakePID(event);
246      
247      rc+=tpcTracker.PropagateBack(event);
248      tpcTracker.UnloadClusters();
249      tpcPID.MakePID(event);
250
251
252      TTree *trdTree=trdl->TreeR();
253      if (!trdTree) {
254         cerr<<"Can't get the TPC cluster tree !\n";
255         return 4;
256      } 
257      trdTracker.LoadClusters(trdTree);
258      rc+=trdTracker.PropagateBack(event);
259      trdTracker.UnloadClusters();
260
261 /*
262      for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
263        AliESDtrack* track = event->GetTrack(iTrack);
264        trdPID->MakePID(track);
265      }
266 */
267
268      TTree *tofTree=tofl->TreeD();
269      if (!tofTree) {
270         cerr<<"Can't get the TOF cluster tree !\n";
271         return 4;
272      } 
273      tofPID.LoadClusters(tofTree);
274      tofPID.MakePID(event);
275      tofPID.UnloadClusters();
276
277
278     //Here is the combined PID
279      AliESDpid::MakePID(event);
280
281      if (rc==0) {
282         Char_t ename[100]; 
283         sprintf(ename,"%d",i);
284         ef->cd();
285         if (!event->Write(ename)) rc++;
286      } 
287      if (rc) {
288         cerr<<"Something bad happened...\n";
289      }
290      delete event;
291    }
292    timer.Stop(); timer.Print();
293
294    //pidFile->Close();
295
296    delete par;
297
298    ef->Close();
299
300    delete rl;
301
302    return rc;
303 }