1 #include "AliHBTReaderInternal.h"
2 //_________________________________________________________________________
3 ///////////////////////////////////////////////////////////////////////////
5 // class AliHBTReaderInternal //
7 // Multi file reader for Internal Data Format //
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 //
13 // Piotr.Skowronski@cern.ch //
14 // more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html //
16 ///////////////////////////////////////////////////////////////////////////
20 #include <TParticle.h>
25 #include "AliHBTRun.h"
26 #include "AliHBTEvent.h"
27 #include "AliHBTParticle.h"
28 #include "AliHBTParticleCut.h"
30 ClassImp(AliHBTReaderInternal)
31 /********************************************************************/
33 AliHBTReaderInternal::AliHBTReaderInternal():
44 /********************************************************************/
46 AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):
56 //filename - name of file to open
58 /********************************************************************/
60 AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
71 //dirs contains strings with directories to look data in
72 //filename - name of file to open
74 /********************************************************************/
76 AliHBTReaderInternal::AliHBTReaderInternal(const AliHBTReaderInternal& in):
78 fFileName(in.fFileName),
88 /********************************************************************/
90 AliHBTReaderInternal& AliHBTReaderInternal::operator=(const AliHBTReaderInternal& in)
93 if (this == &in) return *this;
94 Rewind();//close current session
95 AliHBTReader::operator=((const AliHBTReader&)in);
96 fFileName = in.fFileName;
99 /********************************************************************/
100 AliHBTReaderInternal::~AliHBTReaderInternal()
106 /********************************************************************/
108 void AliHBTReaderInternal::Rewind()
117 /********************************************************************/
119 Int_t AliHBTReaderInternal::ReadNext()
121 //reads data and puts put to the particles and tracks objects
122 //reurns 0 if everything is OK
125 Int_t i; //iterator and some temprary values
127 AliHBTParticle* tpart = 0x0, *ttrack = 0x0;
129 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
132 Error("ReadNext","Can not get PDG Particles Data Base");
136 if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
137 if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
139 fParticlesEvent->Reset();
140 fTracksEvent->Reset();
142 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
145 if( (i=OpenNextFile()) )
147 Error("ReadNext","Skipping directory due to problems with opening files. Errorcode %d",i);
151 if (fCurrentEvent == (Int_t)fTree->GetEntries())
162 /***************************/
163 /***************************/
164 /***************************/
167 Info("ReadNext","Event %d",fCurrentEvent);
168 fTree->GetEvent(fCurrentEvent);
171 if (fPartBranch && fTrackBranch)
173 Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());
174 Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
175 for(i = 0; i < fPartBuffer->GetEntries(); i++)
177 tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
178 ttrack = dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
180 if( tpart == 0x0 ) continue; //if returned pointer is NULL
181 if( ttrack == 0x0 ) continue; //if returned pointer is NULL
183 if (AliHBTParticle::GetDebug() > 9)
185 Info("ReadNext","Particle:");
187 Info("ReadNext","Track:");
190 if (ttrack->GetUID() != tpart->GetUID())
192 Error("ReadNext","Sth. is wrong: Track and Particle has different UID.");
193 Error("ReadNext","They probobly do not correspond to each other.");
196 for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
198 //check if we are intersted with particles of this type
199 //if not take next partilce
200 if( Rejected(ttrack->GetNthPid(s)) )
202 if (AliHBTParticle::GetDebug() > 9)
203 Info("ReadNext","Track Incarnation %d did not pass PID cut.",ttrack->GetNthPid(s));
206 TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
207 if (pdgp == 0x0)//PDG part corresponding to new incarnation
209 Error("ReadNext","Particle code unknown to PDG DB.");
213 AliHBTParticle* track = new AliHBTParticle(*ttrack);
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
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));
228 if( Rejected(track) )
230 if (AliHBTParticle::GetDebug() > 9)
231 Info("ReadNext","Track Incarnation %d did not pass cut.",ttrack->GetNthPid(s));
235 AliHBTParticle* part = new AliHBTParticle(*tpart);//particle has only 1 incarnation (real)
237 fParticlesEvent->AddParticle(part);
238 fTracksEvent->AddParticle(track);
243 Info("ReadNext"," Read: %d particles and tracks.",counter);
249 Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
250 for(i = 0; i < fPartBuffer->GetEntries(); i++)
252 tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
253 if(tpart == 0x0) continue; //if returned pointer is NULL
255 for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
257 if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
258 if( Rejected(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
259 AliHBTParticle* part = new AliHBTParticle(*tpart);
260 part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
266 fParticlesEvent->AddParticle(part);
270 Info("ReadNext"," Read: %d particles.",counter);
272 else Info("ReadNext"," Read: 0 particles.");
276 Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());
277 for(i = 0; i < fTrackBuffer->GetEntries(); i++)
279 ttrack = dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
280 if( ttrack == 0x0 ) continue; //if returned pointer is NULL
282 for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
284 if( Rejected(ttrack->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
285 //if not take next partilce
286 TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
287 if (pdgp == 0x0)//PDG part corresponding to new incarnation
289 Error("ReadNext","Particle code unknown to PDG DB.");
292 AliHBTParticle* track = new AliHBTParticle(*ttrack);
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
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));
307 if( Rejected(track) )
312 fTracksEvent->AddParticle(track);
317 Info("ReadNext"," Read: %d tracks",counter);
319 else Info("ReadNext"," Read: 0 tracks.");
322 if (fPartBuffer) fPartBuffer->Delete();
323 if (fTrackBuffer) fTrackBuffer->Delete();
330 }while(fCurrentDir < GetNumberOfDirs());
332 return 1;//no more directories to read
334 /********************************************************************/
336 Int_t AliHBTReaderInternal::OpenNextFile()
338 //open file in current directory
339 const TString& dirname = GetDirName(fCurrentDir);
342 Error("OpenNextFile","Can not get directory name");
346 TString filename = dirname +"/"+ fFileName;
347 fFile = TFile::Open(filename.Data());
350 Error("OpenNextFile","Can't open file named %s",filename.Data());
353 if (fFile->IsOpen() == kFALSE)
355 Error("OpenNextFile","Can't open filenamed %s",filename.Data());
359 fTree = (TTree*)fFile->Get("data");
362 Error("OpenNextFile","Can not get the tree.");
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");
374 if (fTrackBuffer == 0x0)
375 fTrackBuffer = new TClonesArray("AliHBTParticle",15000);
376 fTrackBranch->SetAddress(&fTrackBuffer);
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");
386 if (fPartBuffer == 0x0)
387 fPartBuffer = new TClonesArray("AliHBTParticle",15000);
388 fPartBranch->SetAddress(&fPartBuffer);
391 Info("OpenNextFile","________________________________________________________");
392 Info("OpenNextFile","Found %d event(s) in directory %s",
393 (Int_t)fTree->GetEntries(),GetDirName(fCurrentDir).Data());
399 /********************************************************************/
401 Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile, Bool_t multcheck)
403 //reads tracks from reader and writes runs to file
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
411 ::Info("AliHBTReaderInternal::Write","________________________________________________________");
412 ::Info("AliHBTReaderInternal::Write","________________________________________________________");
413 ::Info("AliHBTReaderInternal::Write","________________________________________________________");
415 TFile *histoOutput = TFile::Open(outfile,"recreate");
417 if (!histoOutput->IsOpen())
419 ::Error("AliHBTReaderInternal::Write","File is not opened");
423 TTree *tracktree = new TTree("data","Tree with tracks");
425 TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
426 TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
428 TClonesArray &particles = *pbuffer;
429 TClonesArray &tracks = *tbuffer;
431 TString name("Tracks");
433 Int_t nt = reader->GetNumberOfTrackEvents();
434 Int_t np = reader->GetNumberOfPartEvents();
436 if (AliHBTParticle::GetDebug() > 0)
437 ::Info("Write","Reader has %d track events and %d particles events.",nt,np);
439 Bool_t trck = (nt > 0) ? kTRUE : kFALSE;
440 Bool_t part = (np > 0) ? kTRUE : kFALSE;
442 TBranch *trackbranch = 0x0, *partbranch = 0x0;
444 if (trck) trackbranch = tracktree->Branch("tracks",&tbuffer,32000,0);
445 if (part) partbranch = tracktree->Branch("particles",&pbuffer,32000,0);
447 if ( (trck) && (part) && (np != nt))
449 ::Warning("AliHBTReaderInternal::Write","Number of track and particle events is different");
453 if (nt >= np ) n = nt; else n = np;
455 if (AliHBTParticle::GetDebug() > 0)
456 ::Info("Write","Will loop over %d events",n);
458 for ( i =0;i< n; i++)
460 ::Info("AliHBTReaderInternal::Write","Event %d",i+1);
464 AliHBTEvent* trackev = reader->GetTrackEvent(i);
465 for ( j = 0; j< trackev->GetNumberOfParticles();j++)
467 const AliHBTParticle& t = *(trackev->GetParticle(j));
470 if (FindIndex(tbuffer,t.GetUID()))
472 if (AliHBTParticle::GetDebug()>4)
474 ::Info("Write","Track with Event UID %d already stored",t.GetUID());
476 continue; //not to write the same particles with other incarnations
479 new (tracks[counter++]) AliHBTParticle(t);
481 ::Info("AliHBTReaderInternal::Write"," Tracks: %d",tracks.GetEntries());
482 }else ::Info("AliHBTReaderInternal::Write","NO TRACKS");
487 // ::Warning("AliHBTReaderInternal::Write","Find index switched off!!!");
489 AliHBTEvent* partev = reader->GetParticleEvent(i);
490 for ( j = 0; j< partev->GetNumberOfParticles();j++)
492 const AliHBTParticle& part= *(partev->GetParticle(j));
495 if (FindIndex(pbuffer,part.GetUID()))
497 if (AliHBTParticle::GetDebug()>4)
499 ::Info("Write","Particle with Event UID %d already stored",part.GetUID());
501 continue; //not to write the same particles with other incarnations
504 new (particles[counter++]) AliHBTParticle(part);
506 ::Info("AliHBTReaderInternal::Write"," Particles: %d",particles.GetEntries());
507 }else ::Info("AliHBTReaderInternal::Write","NO PARTICLES");
511 tracktree->AutoSave();
517 tracktree->Write(0,TObject::kOverwrite);
525 histoOutput->Close();
528 /********************************************************************/
530 Bool_t AliHBTReaderInternal::FindIndex(TClonesArray* arr,Int_t idx)
532 //Checks if in the array exists already partilce with Unique ID idx
535 ::Error("FindIndex","Array is 0x0");
540 while (( p = (AliHBTParticle*)next()))
542 if (p->GetUID() == idx) return kTRUE;