]>
Commit | Line | Data |
---|---|---|
1b446896 | 1 | #include "AliHBTReaderTPC.h" |
2 | ||
61ab8ea8 | 3 | #include <Riostream.h> |
4 | //#include <Riostream.h> | |
1b446896 | 5 | #include <TTree.h> |
6 | #include <TFile.h> | |
16701d1b | 7 | #include <TParticle.h> |
1b446896 | 8 | |
16701d1b | 9 | #include <AliRun.h> |
10 | #include <AliMagF.h> | |
1b446896 | 11 | #include <AliTPCtrack.h> |
12 | #include <AliTPCParam.h> | |
13 | #include <AliTPCtracker.h> | |
14 | ||
1b446896 | 15 | #include "AliHBTRun.h" |
16 | #include "AliHBTEvent.h" | |
17 | #include "AliHBTParticle.h" | |
18 | #include "AliHBTParticleCut.h" | |
19 | ||
20 | ||
21 | ClassImp(AliHBTReaderTPC) | |
22 | //reader for TPC tracking | |
61ab8ea8 | 23 | //needs galice.root, AliTPCtracks.root, AliTPCclusters.root, good_tracks_tpc |
1b446896 | 24 | // |
25 | //more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html | |
26 | //Piotr.Skowronski@cern.ch | |
27 | ||
61ab8ea8 | 28 | AliHBTReaderTPC:: |
29 | AliHBTReaderTPC(const Char_t* trackfilename,const Char_t* clusterfilename, | |
30 | const Char_t* galicefilename): | |
31 | fTrackFileName(trackfilename),fClusterFileName(clusterfilename), | |
32 | fGAliceFileName(galicefilename) | |
1b446896 | 33 | { |
16701d1b | 34 | //constructor, |
1b446896 | 35 | //Defaults: |
36 | // trackfilename = "AliTPCtracks.root" | |
37 | // clusterfilename = "AliTPCclusters.root" | |
1b446896 | 38 | // galicefilename = "" - this means: Do not open gAlice file - |
39 | // just leave the global pointer untached | |
40 | ||
41 | fParticles = new AliHBTRun(); | |
42 | fTracks = new AliHBTRun(); | |
16701d1b | 43 | fIsRead = kFALSE; |
44 | } | |
45 | /********************************************************************/ | |
61ab8ea8 | 46 | AliHBTReaderTPC:: |
47 | AliHBTReaderTPC(TObjArray* dirs, | |
df9a1768 | 48 | const Char_t* trackfilename, const Char_t* clusterfilename, |
49 | const Char_t* galicefilename): | |
61ab8ea8 | 50 | AliHBTReader(dirs), fTrackFileName(trackfilename), |
51 | fClusterFileName(clusterfilename),fGAliceFileName(galicefilename) | |
52 | ||
16701d1b | 53 | { |
54 | //constructor, | |
55 | //Defaults: | |
56 | // trackfilename = "AliTPCtracks.root" | |
57 | // clusterfilename = "AliTPCclusters.root" | |
58 | // galicefilename = "" - this means: Do not open gAlice file - | |
59 | // just leave the global pointer untached | |
60 | ||
16701d1b | 61 | fParticles = new AliHBTRun(); |
62 | fTracks = new AliHBTRun(); | |
1b446896 | 63 | fIsRead = kFALSE; |
61ab8ea8 | 64 | |
1b446896 | 65 | } |
66 | /********************************************************************/ | |
67 | ||
68 | AliHBTReaderTPC::~AliHBTReaderTPC() | |
69 | { | |
70 | //desctructor | |
71 | delete fParticles; | |
72 | delete fTracks; | |
73 | } | |
74 | /********************************************************************/ | |
75 | ||
76 | AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n) | |
77 | { | |
78 | //returns Nth event with simulated particles | |
16701d1b | 79 | if (!fIsRead) |
80 | if(Read(fParticles,fTracks)) | |
81 | { | |
82 | Error("GetParticleEvent","Error in reading"); | |
83 | return 0x0; | |
84 | } | |
1b446896 | 85 | return fParticles->GetEvent(n); |
86 | } | |
87 | /********************************************************************/ | |
88 | AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n) | |
89 | { | |
90 | //returns Nth event with reconstructed tracks | |
16701d1b | 91 | if (!fIsRead) |
92 | if(Read(fParticles,fTracks)) | |
93 | { | |
94 | Error("GetTrackEvent","Error in reading"); | |
95 | return 0x0; | |
96 | } | |
1b446896 | 97 | return fTracks->GetEvent(n); |
98 | } | |
99 | /********************************************************************/ | |
100 | ||
101 | Int_t AliHBTReaderTPC::GetNumberOfPartEvents() | |
102 | { | |
103 | //returns number of events of particles | |
16701d1b | 104 | if (!fIsRead) |
105 | if ( Read(fParticles,fTracks)) | |
106 | { | |
107 | Error("GetNumberOfPartEvents","Error in reading"); | |
108 | return 0; | |
109 | } | |
1b446896 | 110 | return fParticles->GetNumberOfEvents(); |
111 | } | |
112 | ||
113 | /********************************************************************/ | |
114 | Int_t AliHBTReaderTPC::GetNumberOfTrackEvents() | |
115 | { | |
116 | //returns number of events of tracks | |
16701d1b | 117 | if (!fIsRead) |
118 | if(Read(fParticles,fTracks)) | |
119 | { | |
120 | Error("GetNumberOfTrackEvents","Error in reading"); | |
121 | return 0; | |
122 | } | |
1b446896 | 123 | return fTracks->GetNumberOfEvents(); |
124 | } | |
125 | /********************************************************************/ | |
126 | ||
127 | ||
128 | Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks) | |
129 | { | |
130 | //reads data and puts put to the particles and tracks objects | |
131 | //reurns 0 if everything is OK | |
132 | // | |
61ab8ea8 | 133 | cout<<"AliHBTReaderTPC::Read()"<<endl; |
1b446896 | 134 | Int_t i; //iterator and some temprary values |
16701d1b | 135 | Int_t Nevents = 0; |
136 | Int_t totalNevents = 0; | |
137 | TFile *aTracksFile;//file with tracks | |
138 | TFile *aClustersFile;//file with clusters | |
139 | TFile *aGAliceFile;//!ile name with galice | |
140 | ||
1b446896 | 141 | if (!particles) //check if an object is instatiated |
142 | { | |
7da9bdd2 | 143 | Error("Read"," particles object must instatiated before passing it to the reader"); |
1b446896 | 144 | } |
145 | if (!tracks) //check if an object is instatiated | |
146 | { | |
7da9bdd2 | 147 | Error("Read"," tracks object must instatiated before passing it to the reader"); |
1b446896 | 148 | } |
149 | particles->Reset();//clear runs == delete all old events | |
150 | tracks->Reset(); | |
16701d1b | 151 | |
152 | TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks | |
153 | tarray->SetOwner(); //set the ownership of the objects it contains | |
154 | //when array is is deleted or cleared all objects | |
155 | //that it contains are deleted | |
156 | Int_t currentdir = 0; | |
b71a5edb | 157 | |
158 | Int_t Ndirs; | |
159 | if (fDirs) //if array with directories is supplied by user | |
160 | { | |
161 | Ndirs = fDirs->GetEntries(); //get the number if directories | |
162 | } | |
163 | else | |
164 | { | |
165 | Ndirs = 0; //if the array is not supplied read only from current directory | |
166 | } | |
167 | ||
16701d1b | 168 | do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./" |
1b446896 | 169 | { |
16701d1b | 170 | |
171 | if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) ) | |
172 | { | |
7da9bdd2 | 173 | Error("Read","Exiting due to problems with opening files. Errorcode %d",i); |
5ee8b891 | 174 | currentdir++; |
175 | continue; | |
16701d1b | 176 | } |
1b446896 | 177 | |
178 | ||
16701d1b | 179 | if (gAlice->TreeE())//check if tree E exists |
180 | { | |
181 | Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice | |
61ab8ea8 | 182 | cout<<"________________________________________________________\n"; |
183 | cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl; | |
184 | cout<<"Setting Magnetic Field: B="<<gAlice->Field()->SolenoidField()<<"T"<<endl; | |
5ee8b891 | 185 | AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField()); |
16701d1b | 186 | } |
187 | else | |
188 | {//if not return an error | |
7da9bdd2 | 189 | Error("Read","Can not find Header tree (TreeE) in gAlice"); |
5ee8b891 | 190 | currentdir++; |
191 | continue; | |
16701d1b | 192 | } |
1b446896 | 193 | |
16701d1b | 194 | aClustersFile->cd();//set cluster file active |
2f1287d2 | 195 | AliTPCParam *TPCParam= (AliTPCParam*)aClustersFile->Get("75x40_100x60_150x60"); |
16701d1b | 196 | if (!TPCParam) |
197 | { | |
2f1287d2 | 198 | TPCParam= (AliTPCParam*)aClustersFile->Get("75x40_100x60"); |
199 | if (!TPCParam) | |
200 | { | |
201 | Error("Read","TPC parameters have not been found !\n"); | |
202 | currentdir++; | |
203 | continue; | |
204 | } | |
16701d1b | 205 | } |
1b446896 | 206 | |
1b446896 | 207 | |
16701d1b | 208 | for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events |
209 | { | |
61ab8ea8 | 210 | cout<<"Reading Event "<<currentEvent<<endl; |
16701d1b | 211 | /**************************************/ |
212 | /**************************************/ | |
213 | /**************************************/ | |
1b446896 | 214 | |
16701d1b | 215 | aTracksFile->cd();//set track file active |
216 | ||
1b446896 | 217 | Char_t treename[100]; |
218 | sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree | |
219 | ||
220 | TTree *tracktree=0; | |
221 | ||
16701d1b | 222 | tracktree=(TTree*)aTracksFile->Get(treename);//get the tree |
1b446896 | 223 | if (!tracktree) //check if we got the tree |
224 | {//if not return with error | |
7da9bdd2 | 225 | Error("Read","Can't get a tree with TPC tracks !\n"); |
5ee8b891 | 226 | continue; |
1b446896 | 227 | } |
228 | ||
229 | TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks | |
230 | if (!trackbranch) ////check if we got the branch | |
231 | {//if not return with error | |
7da9bdd2 | 232 | Error("Read","Can't get a branch with TPC tracks !\n"); |
5ee8b891 | 233 | continue; |
1b446896 | 234 | } |
235 | Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks | |
61ab8ea8 | 236 | cout<<"Found "<<NTPCtracks<<" TPC tracks.\n"; |
1b446896 | 237 | //Copy tracks to array |
238 | ||
239 | AliTPCtrack *iotrack=0; | |
240 | ||
16701d1b | 241 | aClustersFile->cd();//set cluster file active |
61ab8ea8 | 242 | //AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);//create the tacker for this event |
243 | AliTPCtracker *tracker = new AliTPCtracker(TPCParam); //I.B. | |
244 | tracker->SetEventNumber(currentEvent); //I.B. | |
1b446896 | 245 | if (!tracker) //check if it has created succeffuly |
246 | {//if not return with error | |
7da9bdd2 | 247 | Error("Read","Can't get a tracker !\n"); |
5ee8b891 | 248 | continue; |
1b446896 | 249 | } |
61ab8ea8 | 250 | //tracker->LoadInnerSectors(); //I.B. |
251 | //tracker->LoadOuterSectors(); //I.B. | |
252 | tracker->LoadClusters(); //I.B. | |
1b446896 | 253 | |
254 | for (i=0; i<NTPCtracks; i++) //loop over all tpc tracks | |
255 | { | |
256 | iotrack=new AliTPCtrack; //create new tracks | |
257 | trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file) | |
258 | tracktree->GetEvent(i); //stream track i to the iotrack | |
259 | tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track | |
260 | //which is the label of corresponding simulated particle | |
261 | tarray->AddLast(iotrack); //put the track in the array | |
262 | } | |
263 | ||
16701d1b | 264 | aTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more |
265 | aTracksFile->Delete("tracks");//delete branch from memmory | |
1b446896 | 266 | delete tracker; //delete tracker |
267 | ||
268 | tracker = 0x0; | |
269 | trackbranch = 0x0; | |
270 | tracktree = 0x0; | |
1b446896 | 271 | |
272 | Double_t xk; | |
273 | Double_t par[5]; | |
274 | Float_t phi, lam, pt;//angles and transverse momentum | |
275 | Int_t label; //label of the current track | |
16701d1b | 276 | |
277 | aGAliceFile->cd(); | |
278 | gAlice->GetEvent(currentEvent); | |
279 | ||
280 | gAlice->Particles(); | |
281 | ||
282 | for (i=0; i<NTPCtracks; i++) //loop over all good tracks | |
1b446896 | 283 | { |
16701d1b | 284 | iotrack = (AliTPCtrack*)tarray->At(i); |
285 | label = iotrack->GetLabel(); | |
286 | ||
287 | if (label < 0) continue; | |
1b446896 | 288 | |
16701d1b | 289 | TParticle *p = (TParticle*)gAlice->Particle(label); |
1b446896 | 290 | |
5ee8b891 | 291 | if(p == 0x0) continue; //if returned pointer is NULL |
292 | if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database) | |
f1e2da22 | 293 | |
16701d1b | 294 | if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type |
295 | //if not take next partilce | |
1b446896 | 296 | |
16701d1b | 297 | AliHBTParticle* part = new AliHBTParticle(*p); |
1b446896 | 298 | if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts |
299 | //if it does not delete it and take next good track | |
300 | ||
16701d1b | 301 | iotrack->PropagateToVertex(); |
7da9bdd2 | 302 | |
1b446896 | 303 | iotrack->GetExternalParameters(xk,par); //get properties of the track |
304 | phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); | |
305 | if (phi<-TMath::Pi()) phi+=2*TMath::Pi(); | |
306 | if (phi>=TMath::Pi()) phi-=2*TMath::Pi(); | |
307 | lam=par[3]; | |
308 | pt=1.0/TMath::Abs(par[4]); | |
309 | ||
310 | Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum | |
311 | Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum | |
312 | Double_t tpz = pt * lam; //track z coordinate of momentum | |
313 | ||
16701d1b | 314 | Double_t mass = p->GetMass(); |
1b446896 | 315 | Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track |
316 | ||
16701d1b | 317 | AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.); |
5ee8b891 | 318 | if(Pass(track))//check if meets all criteria of any of our cuts |
7da9bdd2 | 319 | //if it does not delete it and take next good track |
320 | { | |
321 | delete track; | |
322 | delete part; | |
323 | continue; | |
324 | } | |
16701d1b | 325 | particles->AddParticle(totalNevents,part);//put track and particle on the run |
326 | tracks->AddParticle(totalNevents,track); | |
1b446896 | 327 | |
328 | } | |
329 | tarray->Clear(); //clear the array | |
16701d1b | 330 | |
331 | /**************************************/ | |
1b446896 | 332 | /**************************************/ |
16701d1b | 333 | /**************************************/ |
334 | totalNevents++; | |
335 | } | |
1b446896 | 336 | |
16701d1b | 337 | //save environment (resouces) -- |
338 | //clean your place after the work | |
339 | CloseFiles(aTracksFile,aClustersFile,aGAliceFile); | |
340 | currentdir++; | |
b71a5edb | 341 | }while(currentdir < Ndirs); |
16701d1b | 342 | |
343 | ||
1b446896 | 344 | delete tarray; |
1b446896 | 345 | fIsRead = kTRUE; |
346 | return 0; | |
347 | } | |
348 | ||
349 | /********************************************************************/ | |
16701d1b | 350 | Int_t AliHBTReaderTPC::OpenFiles |
351 | (TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event) | |
1b446896 | 352 | { |
353 | //opens all the files | |
16701d1b | 354 | |
355 | ||
356 | const TString& dirname = GetDirName(event); | |
357 | if (dirname == "") | |
358 | { | |
359 | Error("OpenFiles","Can not get directory name"); | |
360 | return 4; | |
361 | } | |
362 | ||
363 | TString filename = dirname +"/"+ fTrackFileName; | |
364 | aTracksFile = TFile::Open(filename.Data()); | |
365 | if ( aTracksFile == 0x0 ) | |
1b446896 | 366 | { |
16701d1b | 367 | Error("OpenFiles","Can't open file with tacks named %s",filename.Data()); |
1b446896 | 368 | return 1; |
369 | } | |
16701d1b | 370 | if (!aTracksFile->IsOpen()) |
371 | { | |
372 | Error("OpenFiles","Can't open file with tacks named %s",filename.Data()); | |
373 | return 1; | |
374 | } | |
375 | ||
376 | filename = dirname +"/"+ fClusterFileName; | |
377 | aClustersFile = TFile::Open(filename.Data()); | |
378 | if ( aClustersFile == 0x0 ) | |
1b446896 | 379 | { |
16701d1b | 380 | Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data()); |
1b446896 | 381 | return 2; |
382 | } | |
16701d1b | 383 | if (!aClustersFile->IsOpen()) |
384 | { | |
385 | Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data()); | |
386 | return 2; | |
387 | } | |
388 | ||
389 | filename = dirname +"/"+ fGAliceFileName; | |
390 | agAliceFile = TFile::Open(filename.Data()); | |
391 | if ( agAliceFile== 0x0) | |
392 | { | |
393 | Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data()); | |
394 | return 3; | |
395 | } | |
396 | if (!agAliceFile->IsOpen()) | |
397 | { | |
398 | Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data()); | |
399 | return 3; | |
400 | } | |
401 | ||
402 | if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice"))) | |
403 | { | |
404 | Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data()); | |
405 | return 5; | |
406 | } | |
1b446896 | 407 | |
16701d1b | 408 | return 0; |
1b446896 | 409 | } |
16701d1b | 410 | /********************************************************************/ |
1b446896 | 411 | |
412 | /********************************************************************/ | |
413 | ||
16701d1b | 414 | void AliHBTReaderTPC::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile) |
1b446896 | 415 | { |
416 | //closes the files | |
16701d1b | 417 | tracksFile->Close(); |
7da9bdd2 | 418 | delete tracksFile; |
16701d1b | 419 | tracksFile = 0x0; |
420 | clustersFile->Close(); | |
7da9bdd2 | 421 | delete clustersFile; |
16701d1b | 422 | clustersFile = 0x0; |
7da9bdd2 | 423 | |
424 | delete gAlice; | |
425 | gAlice = 0; | |
426 | ||
427 | if (gAliceFile) | |
428 | { | |
429 | gAliceFile->Close(); | |
430 | delete gAliceFile; | |
431 | gAliceFile = 0x0; | |
432 | } | |
1b446896 | 433 | } |
434 | ||
435 | /********************************************************************/ | |
436 | ||
437 | /********************************************************************/ | |
438 | /********************************************************************/ | |
439 | /********************************************************************/ | |
440 |