Primary vertex included in ESD (Yu.Belikov)
[u/mrichter/AliRoot.git] / STEER / AliESDtest.C
CommitLineData
ae982df3 1//********************************************************************
2// Example of the reconstruction that generates the ESD
3// Input files:
c630aafd 4// a) file containing the ITS clusters
ae982df3 5// (the AliITSFindClustersV2.C macro can be used to generate it)
c630aafd 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)
ae982df3 12// Ouput file:
13// AliESDs.root containing the ESD events
14//
15// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
16//********************************************************************
17
c630aafd 18#if !defined(__CINT__) || defined(__MAKECINT__)
ae982df3 19 #include <Riostream.h>
20 #include "TFile.h"
c630aafd 21 #include "TSystem.h"
ae982df3 22 #include "TStopwatch.h"
c630aafd 23 #include "TGeant3.h"
873f1f73 24 #include "TArrayF.h"
c630aafd 25
26 #include "AliMagF.h"
27 #include "AliRun.h"
28 #include "AliRunLoader.h"
29 #include "AliLoader.h"
873f1f73 30 #include "AliHeader.h"
31 #include "AliGenEventHeader.h"
ae982df3 32
33 #include "AliESD.h"
8c6a71ab 34 #include "AliESDpid.h"
c630aafd 35
36 #include "AliITS.h"
ae982df3 37 #include "AliITSgeom.h"
38 #include "AliITStrackerV2.h"
e23730c7 39 #include "AliV0vertexer.h"
40 #include "AliCascadeVertexer.h"
c630aafd 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
79e94bf8 49 #include "AliTRDtracker.h"
50 #include "AliTRDPartID.h"
c630aafd 51
52 #include "AliTOFpidESD.h"
ae982df3 53#endif
54
c630aafd 55extern TSystem *gSystem;
56extern AliRun *gAlice;
57extern TFile *gFile;
58
59Int_t AliESDtest(Int_t nev=1) {
60
61/**** Initialization of the NewIO *******/
79e94bf8 62
c630aafd 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;
ae982df3 87 return 2;
88 }
c630aafd 89 gAlice=rl->GetAliRun();
90
ae982df3 91
c630aafd 92 AliKalmanTrack::SetConvConst(
93 1000/0.299792458/gAlice->Field()->SolenoidField()
94 );
ae982df3 95
96
c630aafd 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;
ae982df3 109 return 4;
110 }
c630aafd 111 AliITSgeom *geom = dITS->GetITSgeom();
ae982df3 112
113 //An instance of the ITS tracker
114 AliITStrackerV2 itsTracker(geom);
8c6a71ab 115
116 //An instance of the ITS PID maker
c630aafd 117 Double_t parITS[]={34.,0.15,10.};
8c6a71ab 118 AliITSpidESD itsPID(parITS);
79e94bf8 119
e23730c7 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);
c630aafd 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";
79e94bf8 155 return 6;
156 }
c630aafd 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");
79e94bf8 174
175 //An instance of the TRD tracker
c630aafd 176 rl->CdGAFile();
177 AliTRDtracker trdTracker(gFile); //galice.root file
79e94bf8 178
c630aafd 179/*
79e94bf8 180 //An instance of the TRD PID maker
c630aafd 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
79e94bf8 211
212 TFile *ef=TFile::Open("AliESDs.root","RECREATE");
ae982df3 213 if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
214
215 TStopwatch timer;
216 Int_t rc=0;
c630aafd 217 if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents();
ae982df3 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();
ae982df3 222
c630aafd 223 rl->GetEvent(i);
224
873f1f73 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
bb2ceb1f 236
237//***** Initial path towards the primary vertex
873f1f73 238 tpcTracker.SetVertex(vtx,cvtx);
c630aafd 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);
ae982df3 245 rc+=tpcTracker.Clusters2Tracks(event);
246
873f1f73 247 itsTracker.SetVertex(vtx,cvtx);
c630aafd 248 TTree *itsTree=itsl->TreeR();
249 if (!itsTree) {
bb2ceb1f 250 cerr<<"Can't get the ITS cluster tree !\n";
c630aafd 251 return 4;
252 }
253 itsTracker.LoadClusters(itsTree);
ae982df3 254 rc+=itsTracker.Clusters2Tracks(event);
255
e23730c7 256
bb2ceb1f 257//***** Back propagation towards the outer barrel detectors
ae982df3 258 rc+=itsTracker.PropagateBack(event);
e23730c7 259 //itsPID.MakePID(event);
ae982df3 260
261 rc+=tpcTracker.PropagateBack(event);
8c6a71ab 262 tpcPID.MakePID(event);
263
c630aafd 264 TTree *trdTree=trdl->TreeR();
265 if (!trdTree) {
bb2ceb1f 266 cerr<<"Can't get the TRD cluster tree !\n";
c630aafd 267 return 4;
268 }
269 trdTracker.LoadClusters(trdTree);
79e94bf8 270 rc+=trdTracker.PropagateBack(event);
c630aafd 271/*
79e94bf8 272 for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
273 AliESDtrack* track = event->GetTrack(iTrack);
274 trdPID->MakePID(track);
275 }
c630aafd 276*/
c630aafd 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
79e94bf8 286
bb2ceb1f 287
288//***** Here is the combined PID
8c6a71ab 289 AliESDpid::MakePID(event);
ae982df3 290
bb2ceb1f 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
873f1f73 306 vtxer.SetVertex(vtx);
bb2ceb1f 307 rc+=vtxer.Tracks2V0vertices(event); // V0 finding
308 rc+=cvtxer.V0sTracks2CascadeVertices(event); // cascade finding
309
310
311
312//***** Some final manipulations with this event
ae982df3 313 if (rc==0) {
314 Char_t ename[100];
315 sprintf(ename,"%d",i);
79e94bf8 316 ef->cd();
ae982df3 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
c630aafd 326 //pidFile->Close();
327
ae982df3 328 delete par;
c630aafd 329
ae982df3 330 ef->Close();
331
c630aafd 332 delete rl;
333
ae982df3 334 return rc;
335}