Compiler Warning removed
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderInternal.cxx
CommitLineData
91923a38 1#include "AliHBTReaderInternal.h"
08d8c545 2//_________________________________________________________________________
3///////////////////////////////////////////////////////////////////////////
4// //
5// class AliHBTReaderInternal //
6// //
7// Multi file reader for Internal Data Format //
8// //
9// This reader reads data created by itself //
10// (method AliHBTReaderInternal::Write) //
11// Data are stored in form of tree of TClonesArray of AliHBTParticle's //
12// //
13// Piotr.Skowronski@cern.ch //
ad4ad1d6 14// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html //
08d8c545 15// //
16///////////////////////////////////////////////////////////////////////////
91923a38 17
91923a38 18#include <TTree.h>
19#include <TFile.h>
20#include <TParticle.h>
f701e886 21#include <TError.h>
91923a38 22#include <AliRun.h>
23#include <AliMagF.h>
24
25#include "AliHBTRun.h"
26#include "AliHBTEvent.h"
27#include "AliHBTParticle.h"
28#include "AliHBTParticleCut.h"
29
91923a38 30ClassImp(AliHBTReaderInternal)
31/********************************************************************/
32
bed069a4 33AliHBTReaderInternal::AliHBTReaderInternal():
34 fFileName(),
35 fPartBranch(0x0),
36 fTrackBranch(0x0),
37 fTree(0x0),
2efd465d 38 fFile(0x0),
39 fPartBuffer(0x0),
40 fTrackBuffer(0x0)
91923a38 41{
f701e886 42//Defalut constructor
91923a38 43}
44/********************************************************************/
45
bed069a4 46AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):
47 fFileName(filename),
48 fPartBranch(0x0),
49 fTrackBranch(0x0),
50 fTree(0x0),
2efd465d 51 fFile(0x0),
52 fPartBuffer(0x0),
53 fTrackBuffer(0x0)
91923a38 54{
f701e886 55//constructor
56//filename - name of file to open
91923a38 57}
58/********************************************************************/
f701e886 59
91923a38 60AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
bed069a4 61 AliHBTReader(dirs),
62 fFileName(filename),
63 fPartBranch(0x0),
64 fTrackBranch(0x0),
65 fTree(0x0),
2efd465d 66 fFile(0x0),
67 fPartBuffer(0x0),
68 fTrackBuffer(0x0)
91923a38 69{
f701e886 70//ctor
71//dirs contains strings with directories to look data in
72//filename - name of file to open
91923a38 73}
f701e886 74/********************************************************************/
75
ad4ad1d6 76AliHBTReaderInternal::AliHBTReaderInternal(const AliHBTReaderInternal& in):
77 AliHBTReader(in),
78 fFileName(in.fFileName),
79 fPartBranch(0x0),
80 fTrackBranch(0x0),
81 fTree(0x0),
82 fFile(0x0),
83 fPartBuffer(0x0),
84 fTrackBuffer(0x0)
85{
86 //cpy constructor
87}
88/********************************************************************/
89
90AliHBTReaderInternal& AliHBTReaderInternal::operator=(const AliHBTReaderInternal& in)
91{
92 //Assigment operator
93 if (this == &in) return *this;
94 Rewind();//close current session
95 AliHBTReader::operator=((const AliHBTReader&)in);
96 fFileName = in.fFileName;
5b28b3e2 97 return *this;
ad4ad1d6 98}
99/********************************************************************/
91923a38 100AliHBTReaderInternal::~AliHBTReaderInternal()
101 {
102 //desctructor
bed069a4 103 delete fTree;
104 delete fFile;
91923a38 105 }
106/********************************************************************/
bed069a4 107
108void AliHBTReaderInternal::Rewind()
109{
110 delete fTree;
111 fTree = 0x0;
112 delete fFile;
113 fFile = 0x0;
114 fCurrentDir = 0;
115 fNEventsRead= 0;
116}
91923a38 117/********************************************************************/
118
bed069a4 119Int_t AliHBTReaderInternal::ReadNext()
91923a38 120 {
121 //reads data and puts put to the particles and tracks objects
122 //reurns 0 if everything is OK
123 //
bed069a4 124 Info("ReadNext","");
91923a38 125 Int_t i; //iterator and some temprary values
f701e886 126 Int_t counter;
9d09c308 127 AliHBTParticle* tpart = 0x0, *ttrack = 0x0;
bed069a4 128
27bb2a54 129 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
130 if (pdgdb == 0x0)
131 {
bed069a4 132 Error("ReadNext","Can not get PDG Particles Data Base");
27bb2a54 133 return 1;
91923a38 134 }
bed069a4 135
136 if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
137 if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
138
139 fParticlesEvent->Reset();
140 fTracksEvent->Reset();
d4c1d7b6 141
91923a38 142 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
143 {
bed069a4 144 if (fTree == 0x0)
145 if( (i=OpenNextFile()) )
146 {
4d91c73a 147 Error("ReadNext","Skipping directory due to problems with opening files. Errorcode %d",i);
bed069a4 148 fCurrentDir++;
149 continue;
150 }
151 if (fCurrentEvent == (Int_t)fTree->GetEntries())
91923a38 152 {
bed069a4 153 delete fTree;
154 fTree = 0x0;
155 delete fFile;
156 fFile = 0x0;
157 fPartBranch = 0x0;
158 fTrackBranch= 0x0;
159 fCurrentDir++;
27bb2a54 160 continue;
91923a38 161 }
162 /***************************/
163 /***************************/
164 /***************************/
165
bed069a4 166
167 Info("ReadNext","Event %d",fCurrentEvent);
168 fTree->GetEvent(fCurrentEvent);
91923a38 169
bed069a4 170 counter = 0;
171 if (fPartBranch && fTrackBranch)
4d1d0fe9 172 {
e4cbdf6b 173 Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());
174 Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
bed069a4 175 for(i = 0; i < fPartBuffer->GetEntries(); i++)
176 {
177 tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
178 ttrack = dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
179
180 if( tpart == 0x0 ) continue; //if returned pointer is NULL
4d1d0fe9 181 if( ttrack == 0x0 ) continue; //if returned pointer is NULL
e4cbdf6b 182
183 if (AliHBTParticle::GetDebug() > 9)
184 {
185 Info("ReadNext","Particle:");
186 tpart->Print();
187 Info("ReadNext","Track:");
188 ttrack->Print();
189 }
bed069a4 190 if (ttrack->GetUID() != tpart->GetUID())
9d09c308 191 {
bed069a4 192 Error("ReadNext","Sth. is wrong: Track and Particle has different UID.");
193 Error("ReadNext","They probobly do not correspond to each other.");
194 }
e4cbdf6b 195
4d1d0fe9 196 for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
bed069a4 197 {
e4cbdf6b 198 //check if we are intersted with particles of this type
199 //if not take next partilce
200 if( Pass(ttrack->GetNthPid(s)) )
201 {
202 if (AliHBTParticle::GetDebug() > 9)
203 Info("ReadNext","Track Incarnation %d did not pass PID cut.",ttrack->GetNthPid(s));
204 continue;
205 }
4d1d0fe9 206 TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
207 if (pdgp == 0x0)//PDG part corresponding to new incarnation
208 {
209 Error("ReadNext","Particle code unknown to PDG DB.");
210 continue;
211 }
7889564f 212
213 AliHBTParticle* track = new AliHBTParticle(*ttrack);
214
215 //apart of setting PDG code of an incarnation
216 //it is necessary tu recalculate energy on the basis of
217 //new PDG code (mass) hypothesis
4d1d0fe9 218 Double_t mass = pdgp->Mass();//mass of new incarnation
219 Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() +
220 ttrack->Py()*ttrack->Py() +
221 ttrack->Pz()*ttrack->Pz() +
222 mass*mass);//total energy of the new incarnation
223 //update Energy and Calc Mass
224 track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
225 track->SetCalcMass(mass);
226 track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
227
228 if( Pass(track) )
27bb2a54 229 {
e4cbdf6b 230 if (AliHBTParticle::GetDebug() > 9)
231 Info("ReadNext","Track Incarnation %d did not pass cut.",ttrack->GetNthPid(s));
4d1d0fe9 232 delete track;
bed069a4 233 continue;
27bb2a54 234 }
4d1d0fe9 235 AliHBTParticle* part = new AliHBTParticle(*tpart);//particle has only 1 incarnation (real)
bed069a4 236
237 fParticlesEvent->AddParticle(part);
238 fTracksEvent->AddParticle(track);
239
240 counter++;
241 }
242 }
243 Info("ReadNext"," Read: %d particles and tracks.",counter);
244 }
245 else
246 {
247 if (fPartBranch)
248 {
249 Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
250 for(i = 0; i < fPartBuffer->GetEntries(); i++)
251 {
252 tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
253 if(tpart == 0x0) continue; //if returned pointer is NULL
254
255 for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
256 {
257 if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
258 if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
bed069a4 259 AliHBTParticle* part = new AliHBTParticle(*tpart);
260 part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
261 if( Pass(part) )
27bb2a54 262 {
bed069a4 263 delete part;
264 continue;
27bb2a54 265 }
bed069a4 266 fParticlesEvent->AddParticle(part);
267 counter++;
268 }
269 }
270 Info("ReadNext"," Read: %d particles.",counter);
271 }
272 else Info("ReadNext"," Read: 0 particles.");
273
274 if (fTrackBranch)
275 {
e4cbdf6b 276 Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());
bed069a4 277 for(i = 0; i < fTrackBuffer->GetEntries(); i++)
278 {
4d1d0fe9 279 ttrack = dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
280 if( ttrack == 0x0 ) continue; //if returned pointer is NULL
bed069a4 281
4d1d0fe9 282 for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
283 {
4d1d0fe9 284 if( Pass(ttrack->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
7889564f 285 //if not take next partilce
4d1d0fe9 286 TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
287 if (pdgp == 0x0)//PDG part corresponding to new incarnation
288 {
289 Error("ReadNext","Particle code unknown to PDG DB.");
290 continue;
291 }
7889564f 292 AliHBTParticle* track = new AliHBTParticle(*ttrack);
293
294 //apart of setting PDG code of an incarnation
295 //it is necessary tu recalculate energy on the basis of
296 //new PDG code (mass) hypothesis
4d1d0fe9 297 Double_t mass = pdgp->Mass();//mass of new incarnation
298 Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() +
299 ttrack->Py()*ttrack->Py() +
300 ttrack->Pz()*ttrack->Pz() +
301 mass*mass);//total energy of the new incarnation
302 //update Energy and Calc Mass
303 track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
304 track->SetCalcMass(mass);
305 track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
306
307 if( Pass(track) )
27bb2a54 308 {
4d1d0fe9 309 delete track;
bed069a4 310 continue;
27bb2a54 311 }
4d1d0fe9 312 fTracksEvent->AddParticle(track);
313
bed069a4 314 counter++;
315 }
316 }
317 Info("ReadNext"," Read: %d tracks",counter);
318 }
319 else Info("ReadNext"," Read: 0 tracks.");
320 }
d7f24a3d 321
322 if (fPartBuffer) fPartBuffer->Delete();
323 if (fTrackBuffer) fTrackBuffer->Delete();
bed069a4 324
325 fCurrentEvent++;
326 fNEventsRead++;
327
328 return 0;
9d09c308 329
4d91c73a 330 }while(fCurrentDir < GetNumberOfDirs());
91923a38 331
bed069a4 332 return 1;//no more directories to read
91923a38 333}
91923a38 334/********************************************************************/
335
bed069a4 336Int_t AliHBTReaderInternal::OpenNextFile()
91923a38 337{
bed069a4 338 //open file in current directory
4d91c73a 339 const TString& dirname = GetDirName(fCurrentDir);
91923a38 340 if (dirname == "")
341 {
bed069a4 342 Error("OpenNextFile","Can not get directory name");
91923a38 343 return 4;
344 }
345
346 TString filename = dirname +"/"+ fFileName;
bed069a4 347 fFile = TFile::Open(filename.Data());
348 if ( fFile == 0x0 )
91923a38 349 {
bed069a4 350 Error("OpenNextFile","Can't open file named %s",filename.Data());
91923a38 351 return 1;
352 }
bed069a4 353 if (fFile->IsOpen() == kFALSE)
91923a38 354 {
bed069a4 355 Error("OpenNextFile","Can't open filenamed %s",filename.Data());
91923a38 356 return 1;
357 }
bed069a4 358
359 fTree = (TTree*)fFile->Get("data");
360 if (fTree == 0x0)
361 {
362 Error("OpenNextFile","Can not get the tree.");
363 return 1;
364 }
365
366
367 fTrackBranch = fTree->GetBranch("tracks");//get the branch with tracks
368 if (fTrackBranch == 0x0) ////check if we got the branch
369 {//if not return with error
370 Info("OpenNextFile","Can't find a branch with tracks !\n");
371 }
372 else
373 {
374 if (fTrackBuffer == 0x0)
375 fTrackBuffer = new TClonesArray("AliHBTParticle",15000);
376 fTrackBranch->SetAddress(&fTrackBuffer);
377 }
378
379 fPartBranch = fTree->GetBranch("particles");//get the branch with particles
380 if (fPartBranch == 0x0) ////check if we got the branch
381 {//if not return with error
382 Info("OpenNextFile","Can't find a branch with particles !\n");
383 }
384 else
385 {
386 if (fPartBuffer == 0x0)
387 fPartBuffer = new TClonesArray("AliHBTParticle",15000);
388 fPartBranch->SetAddress(&fPartBuffer);
389 }
390
391 Info("OpenNextFile","________________________________________________________");
392 Info("OpenNextFile","Found %d event(s) in directory %s",
18c3af34 393 (Int_t)fTree->GetEntries(),GetDirName(fCurrentDir).Data());
bed069a4 394
395 fCurrentEvent = 0;
396
91923a38 397 return 0;
398}
f701e886 399/********************************************************************/
91923a38 400
27bb2a54 401Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile, Bool_t multcheck)
91923a38 402 {
f701e886 403 //reads tracks from reader and writes runs to file
27bb2a54 404 //reader - provides data for writing in internal format
405 //name of output file
406 //multcheck - switches of checking if particle was stored with other incarnation
407 // usefull e.g. when using kine data, where all particles have 100% pid prob.and saves a lot of time
408
f701e886 409 Int_t i,j;
91923a38 410
f701e886 411 ::Info("AliHBTReaderInternal::Write","________________________________________________________");
412 ::Info("AliHBTReaderInternal::Write","________________________________________________________");
413 ::Info("AliHBTReaderInternal::Write","________________________________________________________");
9d09c308 414
91923a38 415 TFile *histoOutput = TFile::Open(outfile,"recreate");
416
417 if (!histoOutput->IsOpen())
418 {
f701e886 419 ::Error("AliHBTReaderInternal::Write","File is not opened");
91923a38 420 return 1;
421 }
422
91923a38 423 TTree *tracktree = new TTree("data","Tree with tracks");
424
425 TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
426 TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
91923a38 427
428 TClonesArray &particles = *pbuffer;
429 TClonesArray &tracks = *tbuffer;
430
431 TString name("Tracks");
91923a38 432
27bb2a54 433 Int_t nt = reader->GetNumberOfTrackEvents();
434 Int_t np = reader->GetNumberOfPartEvents();
91923a38 435
27bb2a54 436 if (AliHBTParticle::GetDebug() > 0)
437 ::Info("Write","Reader has %d track events and %d particles events.",nt,np);
438
439 Bool_t trck = (nt > 0) ? kTRUE : kFALSE;
440 Bool_t part = (np > 0) ? kTRUE : kFALSE;
91923a38 441
442 TBranch *trackbranch = 0x0, *partbranch = 0x0;
443
d4c1d7b6 444 if (trck) trackbranch = tracktree->Branch("tracks",&tbuffer,32000,0);
445 if (part) partbranch = tracktree->Branch("particles",&pbuffer,32000,0);
91923a38 446
27bb2a54 447 if ( (trck) && (part) && (np != nt))
91923a38 448 {
f701e886 449 ::Warning("AliHBTReaderInternal::Write","Number of track and particle events is different");
91923a38 450 }
451
27bb2a54 452 Int_t n;
453 if (nt >= np ) n = nt; else n = np;
454
455 if (AliHBTParticle::GetDebug() > 0)
456 ::Info("Write","Will loop over %d events",n);
91923a38 457
27bb2a54 458 for ( i =0;i< n; i++)
91923a38 459 {
f701e886 460 ::Info("AliHBTReaderInternal::Write","Event %d",i+1);
27bb2a54 461 Int_t counter = 0;
462 if (trck && (i<=nt))
463 {
91923a38 464 AliHBTEvent* trackev = reader->GetTrackEvent(i);
465 for ( j = 0; j< trackev->GetNumberOfParticles();j++)
466 {
27bb2a54 467 const AliHBTParticle& t = *(trackev->GetParticle(j));
468 if (multcheck)
469 {
470 if (FindIndex(tbuffer,t.GetUID()))
471 {
472 if (AliHBTParticle::GetDebug()>4)
473 {
474 ::Info("Write","Track with Event UID %d already stored",t.GetUID());
475 }
476 continue; //not to write the same particles with other incarnations
477 }
478 }
479 new (tracks[counter++]) AliHBTParticle(t);
91923a38 480 }
27bb2a54 481 ::Info("AliHBTReaderInternal::Write"," Tracks: %d",tracks.GetEntries());
f701e886 482 }else ::Info("AliHBTReaderInternal::Write","NO TRACKS");
91923a38 483
27bb2a54 484 counter = 0;
485 if (part && (i<=np))
91923a38 486 {
27bb2a54 487// ::Warning("AliHBTReaderInternal::Write","Find index switched off!!!");
488
91923a38 489 AliHBTEvent* partev = reader->GetParticleEvent(i);
91923a38 490 for ( j = 0; j< partev->GetNumberOfParticles();j++)
491 {
9d09c308 492 const AliHBTParticle& part= *(partev->GetParticle(j));
27bb2a54 493 if (multcheck)
494 {
495 if (FindIndex(pbuffer,part.GetUID()))
496 {
497 if (AliHBTParticle::GetDebug()>4)
498 {
499 ::Info("Write","Particle with Event UID %d already stored",part.GetUID());
500 }
501 continue; //not to write the same particles with other incarnations
502 }
503 }
504 new (particles[counter++]) AliHBTParticle(part);
91923a38 505 }
27bb2a54 506 ::Info("AliHBTReaderInternal::Write"," Particles: %d",particles.GetEntries());
f701e886 507 }else ::Info("AliHBTReaderInternal::Write","NO PARTICLES");
91923a38 508
91923a38 509 histoOutput->cd();
510 tracktree->Fill();
9d09c308 511 tracktree->AutoSave();
91923a38 512 tbuffer->Delete();
513 pbuffer->Delete();
91923a38 514 }
515
91923a38 516 histoOutput->cd();
91923a38 517 tracktree->Write(0,TObject::kOverwrite);
9d09c308 518 delete tracktree;
91923a38 519
9d09c308 520 tbuffer->SetOwner();
521 pbuffer->SetOwner();
522 delete pbuffer;
523 delete tbuffer;
524
525 histoOutput->Close();
91923a38 526 return 0;
527 }
27bb2a54 528/********************************************************************/
529
530Bool_t AliHBTReaderInternal::FindIndex(TClonesArray* arr,Int_t idx)
531{
532//Checks if in the array exists already partilce with Unique ID idx
533 if (arr == 0x0)
534 {
535 ::Error("FindIndex","Array is 0x0");
536 return kTRUE;
537 }
538 TIter next(arr);
539 AliHBTParticle* p;
540 while (( p = (AliHBTParticle*)next()))
541 {
542 if (p->GetUID() == idx) return kTRUE;
543 }
544 return kFALSE;
545}