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