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