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