1 #include "AliHBTReaderInternal.h"
10 #include "AliHBTRun.h"
11 #include "AliHBTEvent.h"
12 #include "AliHBTParticle.h"
13 #include "AliHBTParticleCut.h"
15 AliHBTReaderInternal test;
17 ClassImp(AliHBTReaderInternal)
18 /********************************************************************/
20 AliHBTReaderInternal::AliHBTReaderInternal():
29 /********************************************************************/
31 AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):
39 //filename - name of file to open
41 /********************************************************************/
43 AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
52 //dirs contains strings with directories to look data in
53 //filename - name of file to open
55 /********************************************************************/
57 AliHBTReaderInternal::~AliHBTReaderInternal()
63 /********************************************************************/
65 void AliHBTReaderInternal::Rewind()
74 /********************************************************************/
76 Int_t AliHBTReaderInternal::ReadNext()
78 //reads data and puts put to the particles and tracks objects
79 //reurns 0 if everything is OK
82 Int_t i; //iterator and some temprary values
84 AliHBTParticle* tpart = 0x0, *ttrack = 0x0;
86 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
89 Error("ReadNext","Can not get PDG Particles Data Base");
93 if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
94 if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
96 fParticlesEvent->Reset();
97 fTracksEvent->Reset();
99 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
102 if( (i=OpenNextFile()) )
104 Error("ReadNext","Skipping directory due to problems with opening files. Errorcode %d",i);
108 if (fCurrentEvent == (Int_t)fTree->GetEntries())
119 /***************************/
120 /***************************/
121 /***************************/
124 Info("ReadNext","Event %d",fCurrentEvent);
125 fTree->GetEvent(fCurrentEvent);
128 if (fPartBranch && fTrackBranch)
130 Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());
131 Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
132 for(i = 0; i < fPartBuffer->GetEntries(); i++)
134 tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
135 ttrack = dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
137 if( tpart == 0x0 ) continue; //if returned pointer is NULL
138 if( ttrack == 0x0 ) continue; //if returned pointer is NULL
140 if (AliHBTParticle::GetDebug() > 9)
142 Info("ReadNext","Particle:");
144 Info("ReadNext","Track:");
147 if (ttrack->GetUID() != tpart->GetUID())
149 Error("ReadNext","Sth. is wrong: Track and Particle has different UID.");
150 Error("ReadNext","They probobly do not correspond to each other.");
153 for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
155 //check if we are intersted with particles of this type
156 //if not take next partilce
157 if( Pass(ttrack->GetNthPid(s)) )
159 if (AliHBTParticle::GetDebug() > 9)
160 Info("ReadNext","Track Incarnation %d did not pass PID cut.",ttrack->GetNthPid(s));
163 TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
164 if (pdgp == 0x0)//PDG part corresponding to new incarnation
166 Error("ReadNext","Particle code unknown to PDG DB.");
170 AliHBTParticle* track = new AliHBTParticle(*ttrack);
172 //apart of setting PDG code of an incarnation
173 //it is necessary tu recalculate energy on the basis of
174 //new PDG code (mass) hypothesis
175 Double_t mass = pdgp->Mass();//mass of new incarnation
176 Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() +
177 ttrack->Py()*ttrack->Py() +
178 ttrack->Pz()*ttrack->Pz() +
179 mass*mass);//total energy of the new incarnation
180 //update Energy and Calc Mass
181 track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
182 track->SetCalcMass(mass);
183 track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
187 if (AliHBTParticle::GetDebug() > 9)
188 Info("ReadNext","Track Incarnation %d did not pass cut.",ttrack->GetNthPid(s));
192 AliHBTParticle* part = new AliHBTParticle(*tpart);//particle has only 1 incarnation (real)
194 fParticlesEvent->AddParticle(part);
195 fTracksEvent->AddParticle(track);
200 Info("ReadNext"," Read: %d particles and tracks.",counter);
206 Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
207 for(i = 0; i < fPartBuffer->GetEntries(); i++)
209 tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
210 if(tpart == 0x0) continue; //if returned pointer is NULL
212 for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
214 if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
215 if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
216 AliHBTParticle* part = new AliHBTParticle(*tpart);
217 part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
223 fParticlesEvent->AddParticle(part);
227 Info("ReadNext"," Read: %d particles.",counter);
229 else Info("ReadNext"," Read: 0 particles.");
233 Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());
234 for(i = 0; i < fTrackBuffer->GetEntries(); i++)
236 ttrack = dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
237 if( ttrack == 0x0 ) continue; //if returned pointer is NULL
239 for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
241 if( Pass(ttrack->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
242 //if not take next partilce
243 TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
244 if (pdgp == 0x0)//PDG part corresponding to new incarnation
246 Error("ReadNext","Particle code unknown to PDG DB.");
249 AliHBTParticle* track = new AliHBTParticle(*ttrack);
251 //apart of setting PDG code of an incarnation
252 //it is necessary tu recalculate energy on the basis of
253 //new PDG code (mass) hypothesis
254 Double_t mass = pdgp->Mass();//mass of new incarnation
255 Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() +
256 ttrack->Py()*ttrack->Py() +
257 ttrack->Pz()*ttrack->Pz() +
258 mass*mass);//total energy of the new incarnation
259 //update Energy and Calc Mass
260 track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
261 track->SetCalcMass(mass);
262 track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
269 fTracksEvent->AddParticle(track);
274 Info("ReadNext"," Read: %d tracks",counter);
276 else Info("ReadNext"," Read: 0 tracks.");
279 if (fPartBuffer) fPartBuffer->Delete();
280 if (fTrackBuffer) fTrackBuffer->Delete();
287 }while(fCurrentDir < GetNumberOfDirs());
289 return 1;//no more directories to read
291 /********************************************************************/
293 Int_t AliHBTReaderInternal::OpenNextFile()
295 //open file in current directory
296 const TString& dirname = GetDirName(fCurrentDir);
299 Error("OpenNextFile","Can not get directory name");
303 TString filename = dirname +"/"+ fFileName;
304 fFile = TFile::Open(filename.Data());
307 Error("OpenNextFile","Can't open file named %s",filename.Data());
310 if (fFile->IsOpen() == kFALSE)
312 Error("OpenNextFile","Can't open filenamed %s",filename.Data());
316 fTree = (TTree*)fFile->Get("data");
319 Error("OpenNextFile","Can not get the tree.");
324 fTrackBranch = fTree->GetBranch("tracks");//get the branch with tracks
325 if (fTrackBranch == 0x0) ////check if we got the branch
326 {//if not return with error
327 Info("OpenNextFile","Can't find a branch with tracks !\n");
331 if (fTrackBuffer == 0x0)
332 fTrackBuffer = new TClonesArray("AliHBTParticle",15000);
333 fTrackBranch->SetAddress(&fTrackBuffer);
336 fPartBranch = fTree->GetBranch("particles");//get the branch with particles
337 if (fPartBranch == 0x0) ////check if we got the branch
338 {//if not return with error
339 Info("OpenNextFile","Can't find a branch with particles !\n");
343 if (fPartBuffer == 0x0)
344 fPartBuffer = new TClonesArray("AliHBTParticle",15000);
345 fPartBranch->SetAddress(&fPartBuffer);
348 Info("OpenNextFile","________________________________________________________");
349 Info("OpenNextFile","Found %d event(s) in directory %s",
350 (Int_t)fTree->GetEntries(),GetDirName(fCurrentDir).Data());
356 /********************************************************************/
358 Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile, Bool_t multcheck)
360 //reads tracks from reader and writes runs to file
361 //reader - provides data for writing in internal format
362 //name of output file
363 //multcheck - switches of checking if particle was stored with other incarnation
364 // usefull e.g. when using kine data, where all particles have 100% pid prob.and saves a lot of time
368 ::Info("AliHBTReaderInternal::Write","________________________________________________________");
369 ::Info("AliHBTReaderInternal::Write","________________________________________________________");
370 ::Info("AliHBTReaderInternal::Write","________________________________________________________");
372 TFile *histoOutput = TFile::Open(outfile,"recreate");
374 if (!histoOutput->IsOpen())
376 ::Error("AliHBTReaderInternal::Write","File is not opened");
380 TTree *tracktree = new TTree("data","Tree with tracks");
382 TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
383 TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
385 TClonesArray &particles = *pbuffer;
386 TClonesArray &tracks = *tbuffer;
388 TString name("Tracks");
390 Int_t nt = reader->GetNumberOfTrackEvents();
391 Int_t np = reader->GetNumberOfPartEvents();
393 if (AliHBTParticle::GetDebug() > 0)
394 ::Info("Write","Reader has %d track events and %d particles events.",nt,np);
396 Bool_t trck = (nt > 0) ? kTRUE : kFALSE;
397 Bool_t part = (np > 0) ? kTRUE : kFALSE;
399 TBranch *trackbranch = 0x0, *partbranch = 0x0;
401 if (trck) trackbranch = tracktree->Branch("tracks",&tbuffer,32000,0);
402 if (part) partbranch = tracktree->Branch("particles",&pbuffer,32000,0);
404 if ( (trck) && (part) && (np != nt))
406 ::Warning("AliHBTReaderInternal::Write","Number of track and particle events is different");
410 if (nt >= np ) n = nt; else n = np;
412 if (AliHBTParticle::GetDebug() > 0)
413 ::Info("Write","Will loop over %d events",n);
415 for ( i =0;i< n; i++)
417 ::Info("AliHBTReaderInternal::Write","Event %d",i+1);
421 AliHBTEvent* trackev = reader->GetTrackEvent(i);
422 for ( j = 0; j< trackev->GetNumberOfParticles();j++)
424 const AliHBTParticle& t = *(trackev->GetParticle(j));
427 if (FindIndex(tbuffer,t.GetUID()))
429 if (AliHBTParticle::GetDebug()>4)
431 ::Info("Write","Track with Event UID %d already stored",t.GetUID());
433 continue; //not to write the same particles with other incarnations
436 new (tracks[counter++]) AliHBTParticle(t);
438 ::Info("AliHBTReaderInternal::Write"," Tracks: %d",tracks.GetEntries());
439 }else ::Info("AliHBTReaderInternal::Write","NO TRACKS");
444 // ::Warning("AliHBTReaderInternal::Write","Find index switched off!!!");
446 AliHBTEvent* partev = reader->GetParticleEvent(i);
447 for ( j = 0; j< partev->GetNumberOfParticles();j++)
449 const AliHBTParticle& part= *(partev->GetParticle(j));
452 if (FindIndex(pbuffer,part.GetUID()))
454 if (AliHBTParticle::GetDebug()>4)
456 ::Info("Write","Particle with Event UID %d already stored",part.GetUID());
458 continue; //not to write the same particles with other incarnations
461 new (particles[counter++]) AliHBTParticle(part);
463 ::Info("AliHBTReaderInternal::Write"," Particles: %d",particles.GetEntries());
464 }else ::Info("AliHBTReaderInternal::Write","NO PARTICLES");
468 tracktree->AutoSave();
474 tracktree->Write(0,TObject::kOverwrite);
482 histoOutput->Close();
485 /********************************************************************/
487 Bool_t AliHBTReaderInternal::FindIndex(TClonesArray* arr,Int_t idx)
489 //Checks if in the array exists already partilce with Unique ID idx
492 ::Error("FindIndex","Array is 0x0");
497 while (( p = (AliHBTParticle*)next()))
499 if (p->GetUID() == idx) return kTRUE;