#include "AliHBTPairCut.h"
#include "AliHBTFunction.h"
#include "AliHBTMonitorFunction.h"
+#include "AliHBTEventBuffer.h"
#include <TBenchmark.h>
#include <TList.h>
Bool_t nonid = IsNonIdentAnalysis();
+ Init();
+
if(oT && oP)
{
- if (fReader->GetNumberOfPartEvents() <1)
- {
- Error("Process","There is no Particles. Maybe change the option?");
- return;
- }
- if (fReader->GetNumberOfTrackEvents() <1)
- {
- Error("Process","There is no Tracks. Maybe change the option?");
- return;
- }
-
- if ( RunCoherencyCheck() )
- {
- Error("Process",
- "Coherency check not passed. Maybe change the option?\n");
- return;
- }
- Init();
if (nonid) ProcessTracksAndParticlesNonIdentAnal();
else ProcessTracksAndParticles();
return;
if(oT)
{
- if (fReader->GetNumberOfTrackEvents() <1)
- {
- Error("Process","There is no data to analyze.");
- return;
- }
- Init();
if (nonid) ProcessTracksNonIdentAnal();
else ProcessTracks();
return;
if(oP)
{
- if (fReader->GetNumberOfPartEvents() <1)
- {
- Error("Process","There is no data to analyze.");
- return;
- }
- Init();
if (nonid) ProcessParticlesNonIdentAnal();
else ProcessParticles();
return;
AliHBTParticle * track1, * track2;
AliHBTEvent * trackEvent, *partEvent;
+ AliHBTEvent * trackEvent1 = 0x0,*partEvent1 = 0x0;
AliHBTEvent * trackEvent2,*partEvent2;
// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
- Int_t nev = fReader->GetNumberOfTrackEvents();
-
- /***************************************/
- /****** Looping same events ********/
- /****** filling numerators ********/
- /***************************************/
+// Int_t nev = fReader->GetNumberOfTrackEvents();
AliHBTPair * trackpair = new AliHBTPair();
AliHBTPair * partpair = new AliHBTPair();
AliHBTPair * tmptrackpair;//temprary pointers to pairs
AliHBTPair * tmppartpair;
+ AliHBTEventBuffer partbuffer(fBufferSize);
+ AliHBTEventBuffer trackbuffer(fBufferSize);
+
+
register UInt_t ii;
- for (Int_t i = 0;i<nev;i++)
+ Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
+
+ fReader->Rewind();
+ Int_t i = -1;
+ while (fReader->Next() == kFALSE)
{
- partEvent= fReader->GetParticleEvent(i);
- trackEvent = fReader->GetTrackEvent(i);
+ i++;
+ partEvent= fReader->GetParticleEvent();
+ trackEvent = fReader->GetTrackEvent();
- if (!partEvent) continue;
-
- //N = 0;
+ if ( !partEvent || !trackEvent )
+ {
+ Error("ProcessTracksAndParticles","Can not get event");
+ return;
+ }
+
+ if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
+ {
+ Fatal("ProcessTracksAndParticles",
+ "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
+ i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
+ return;
+ }
- for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+ if(partEvent1 == 0x0)
{
+ partEvent1 = new AliHBTEvent();
+ partEvent1->SetOwner(kTRUE);
- if ( (j%fDisplayMixingInfo) == 0)
+ trackEvent1 = new AliHBTEvent();
+ trackEvent1->SetOwner(kTRUE);
+ }
+ else
+ {
+ partEvent1->Reset();
+ trackEvent1->Reset();
+ }
+
+ for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+ {
+ /***************************************/
+ /****** Looping same events ********/
+ /****** filling numerators ********/
+ /***************************************/
+ if ( (j%fDisplayMixingInfo) == 0)
Info("ProcessTracksAndParticles",
"Mixing particle %d from event %d with particles from event %d",j,i,i);
part1= partEvent->GetParticle(j);
- track1= trackEvent->GetParticle(j);
-
- if (fPairCut->GetFirstPartCut()->Pass(part1)) continue;
+ track1= trackEvent->GetParticle(j);
+
+ if( part1->GetPdgCode() != track1->GetPdgCode() )
+ {
+ Fatal("ProcessTracksAndParticles",
+ "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
+ i,j, part1->GetPdgCode(),track1->GetPdgCode() );
+ return;
+ }
+
+ Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1);
+
+ if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) )
+ {
+ //accepted by any cut
+ // we have to copy because reader keeps only one event
+
+ partEvent1->AddParticle(new AliHBTParticle(*part1));
+ trackEvent1->AddParticle(new AliHBTParticle(*track1));
+ }
+ if (firstcut) continue;
+
for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
fParticleMonitorFunctions[ii]->Process(part1);
for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
- if ( (fNParticleFunctions == 0) && (fNTrackFunctions ==0) && (fNParticleAndTrackFunctions == 0))
- continue;
+ if (nocorrfctns) continue;
for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
{
track2= trackEvent->GetParticle(k);
trackpair->SetParticles(track1,track2);
- if(fPairCut->Pass(partpair) ) //check pair cut
- { //do not meets crietria of the pair cut, try with swaped pairs
+ if(fPairCut->Pass(partpair) ) //check pair cut
+ { //do not meets crietria of the pair cut, try with swapped pairs
if( fPairCut->Pass(partpair->GetSwapedPair()) )
continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
else
- { //swaped pair meets all the criteria
+ { //swaped pair meets all the criteria
tmppartpair = partpair->GetSwapedPair();
tmptrackpair = trackpair->GetSwapedPair();
}
for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
- }
- }
- }
-
-
- /***** Filling denominators *********/
- /***************************************/
- if (fBufferSize != 0)
- for (Int_t i = 0;i<nev-1;i++) //In each event (but last) ....
- {
-
- if ((fNParticleFunctions == 0) && (fNTrackFunctions ==0) && (fNParticleAndTrackFunctions == 0))
- continue;
-
- partEvent= fReader->GetParticleEvent(i);
- if (!partEvent) continue;
-
- trackEvent = fReader->GetTrackEvent(i);
-
- for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
- {
- part1= partEvent->GetParticle(j);
- track1= trackEvent->GetParticle(j);
-
- if (fPairCut->GetFirstPartCut()->Pass(part1)) continue;
-
- Int_t maxeventnumber;
-
- if ( ((i+fBufferSize) >= nev) ||( fBufferSize < 0) ) //if buffer size is negative
- //or current event+buffersize is greater
- //than max nuber of events
- {
- maxeventnumber = nev; //set the max event number
- }
- else
- {
- maxeventnumber = i+fBufferSize; //set the current event number + fBufferSize
- }
-
- for (Int_t k = i+1; k<maxeventnumber;k++) // ... Loop over next event
- {
-
- partEvent2= fReader->GetParticleEvent(k);
- if (!partEvent2) continue;
-
- trackEvent2 = fReader->GetTrackEvent(k);
-
- if ( (j%fDisplayMixingInfo) == 0)
- Info("ProcessTracksAndParticles",
- "Mixing particle %d from event %d with particles from event %d",j,i,k);
+ //end of 2nd loop over particles from the same event
+ }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+
+ /***************************************/
+ /***** Filling denominators *********/
+ /***************************************/
+ if (fBufferSize == 0) continue;
+
+ partbuffer.ResetIter();
+ trackbuffer.ResetIter();
+ Int_t m = 0;
+ while (( partEvent2 = partbuffer.Next() ))
+ {
+ trackEvent2 = trackbuffer.Next();
- for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
+ m++;
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessTracksAndParticles",
+ "Mixing particle %d from event %d with particles from event %d",j,i,i-m);
+
+ for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
{
part2= partEvent2->GetParticle(l);
partpair->SetParticles(part1,part2);
for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
}//for(Int_t l = 0; l<N2;l++) // ... on all particles
- }//for (Int_t k = i+1; k<maxeventnumber;k++) // ... Loop over next event
- }
- }
- /***************************************/
+
+ }
+ //end of loop over particles from first event
+ }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+ partEvent1 = partbuffer.Push(partEvent1);
+ trackEvent1 = trackbuffer.Push(trackEvent1);
+ //end of loop over events
+ }//while (fReader->Next() == kFALSE)
}
/*************************************************************************************/
//the loops are splited
AliHBTParticle * track1, * track2;
AliHBTEvent * trackEvent;
+ AliHBTEvent * trackEvent1 = 0x0;
AliHBTEvent * trackEvent2;
- UInt_t ii;
- Int_t nev = fReader->GetNumberOfTrackEvents();
+ register UInt_t ii;
AliHBTPair * trackpair = new AliHBTPair();
AliHBTPair * tmptrackpair; //temporary pointer
- /***************************************/
- /****** Looping same events ********/
- /****** filling numerators ********/
- /***************************************/
- for (Int_t i = 0;i<nev;i++)
+ AliHBTEventBuffer trackbuffer(fBufferSize);
+
+ fReader->Rewind();
+ Int_t i = -1;
+ while (fReader->Next() == kFALSE)
{
- trackEvent = fReader->GetTrackEvent(i);
+ i++;
+ trackEvent = fReader->GetTrackEvent();
if (!trackEvent) continue;
+ if(trackEvent1 == 0x0)
+ {
+ trackEvent1 = new AliHBTEvent();
+ trackEvent1->SetOwner(kTRUE);
+ }
+ else
+ {
+ trackEvent1->Reset();
+ }
+
for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
{
+ /***************************************/
+ /****** Looping same events ********/
+ /****** filling numerators ********/
+ /***************************************/
if ( (j%fDisplayMixingInfo) == 0)
Info("ProcessTracks",
"Mixing particle %d from event %d with particles from event %d",j,i,i);
- track1= trackEvent->GetParticle(j);
- if (fPairCut->GetFirstPartCut()->Pass(track1)) continue;
+ track1= trackEvent->GetParticle(j);
+ Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(track1);
+
+ if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(track1) == kFALSE) )
+ {
+ //accepted by any cut
+ // we have to copy because reader keeps only one event
+ trackEvent1->AddParticle(new AliHBTParticle(*track1));
+ }
+
+ if (firstcut) continue;
for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
fTrackMonitorFunctions[ii]->Process(track1);
- if ( fNTrackFunctions ==0 )
- continue;
+ if ( fNTrackFunctions ==0 ) continue;
for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
{
for(ii = 0;ii<fNTrackFunctions;ii++)
fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
}
- }
- }
-
- /***************************************/
- /***** Filling diff histogram *********/
- /***************************************/
- if (fBufferSize != 0)
- for (Int_t i = 0;i<nev-1;i++) //In each event (but last) ....
- {
- if ( fNTrackFunctions ==0 )
- continue;
-
- trackEvent = fReader->GetTrackEvent(i);
- if (!trackEvent) continue;
-
- for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
- {
- track1= trackEvent->GetParticle(j);
- if (fPairCut->GetFirstPartCut()->Pass(track1)) continue;
-
- Int_t maxeventnumber;
-
- if ( ((i+fBufferSize) >= nev) ||( fBufferSize < 0) ) //if buffer size is negative
- //or current event+buffersize is greater
- //than max nuber of events
- {
- maxeventnumber = nev; //set the max event number
- }
- else
- {
- maxeventnumber = i+fBufferSize; //set the current event number + fBufferSize
- }
-
- for (Int_t k = i+1; k<maxeventnumber;k++) // ... Loop over next event
+ /***************************************/
+ /***** Filling denominators *********/
+ /***************************************/
+
+ if (fBufferSize == 0) continue;
+
+ trackbuffer.ResetIter();
+
+ while (( trackEvent2 = trackbuffer.Next() ))
{
- trackEvent2 = fReader->GetTrackEvent(k);
- if (!trackEvent2) continue;
-
- if ( (j%fDisplayMixingInfo) == 0)
- Info("ProcessTracks",
- "Mixing particle %d from event %d with particles from event %d",j,i,k);
-
for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all particles
- {
- track2= trackEvent2->GetParticle(l);
- trackpair->SetParticles(track1,track2);
-
- if(fPairCut->Pass(trackpair)) //check pair cut
- { //do not meets crietria of the
- if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
- else tmptrackpair = trackpair->GetSwapedPair();
- }
- else
- {
+ {
+
+ track2= trackEvent2->GetParticle(l);
+ trackpair->SetParticles(track1,track2);
+
+ if( fPairCut->Pass(trackpair) ) //check pair cut
+ { //do not meets crietria of the
+ if( fPairCut->Pass(trackpair->GetSwapedPair()) )
+ continue;
+ else
+ {
+ tmptrackpair = trackpair->GetSwapedPair();
+ }
+ }
+ else
+ {//meets criteria of the pair cut
tmptrackpair = trackpair;
- }
- for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
-
+ }
+
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
+
}//for(Int_t l = 0; l<N2;l++) // ... on all particles
- }//for (Int_t k = i+1; k<maxeventnumber;k++) // ... Loop over next event
+
+ }
}
- }
- /***************************************/
+ trackEvent1 = trackbuffer.Push(trackEvent1);
+ }//while (fReader->Next() == kFALSE)
}
/*************************************************************************************/
//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
//the loops are splited
AliHBTParticle * part1, * part2;
-
AliHBTEvent * partEvent;
+ AliHBTEvent * partEvent1 = 0x0;
AliHBTEvent * partEvent2;
+ register UInt_t ii;
+
AliHBTPair * partpair = new AliHBTPair();
- AliHBTPair * tmppartpair; //temporary pointer to the pair
+ AliHBTPair * tmppartpair; //temporary pointer
- Int_t nev = fReader->GetNumberOfPartEvents();
+ AliHBTEventBuffer partbuffer(fBufferSize);
- /***************************************/
- /****** Looping same events ********/
- /****** filling numerators ********/
- /***************************************/
- for (Int_t i = 0;i<nev;i++)
+ fReader->Rewind();
+ Int_t i = -1;
+ while (fReader->Next() == kFALSE)
{
- partEvent= fReader->GetParticleEvent(i);
+ i++;
+ partEvent = fReader->GetParticleEvent();
if (!partEvent) continue;
- //N = 0;
+
+ if(partEvent1 == 0x0)
+ {
+ partEvent1 = new AliHBTEvent();
+ partEvent1->SetOwner(kTRUE);
+ }
+ else
+ {
+ partEvent1->Reset();
+ }
for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
{
+ /***************************************/
+ /****** Looping same events ********/
+ /****** filling numerators ********/
+ /***************************************/
if ( (j%fDisplayMixingInfo) == 0)
- Info("ProcessParticles",
- "Mixing particle %d from event %d with particles from event %d",j,i,i);
-
+ Info("ProcessParts",
+ "Mixing particle %d from event %d with particles from event %d",j,i,i);
+
part1= partEvent->GetParticle(j);
- if (fPairCut->GetFirstPartCut()->Pass(part1)) continue;
-
- UInt_t zz;
- for(zz = 0; zz<fNParticleMonitorFunctions; zz++)
- fParticleMonitorFunctions[zz]->Process(part1);
-
- if ( fNParticleFunctions ==0 )
- continue;
-
- for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+ Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1);
+
+ if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) )
{
- part2= partEvent->GetParticle(k);
- if (part1->GetUID() == part2->GetUID()) continue; //this is the same particle but different incarnation (PID)
- partpair->SetParticles(part1,part2);
-
- if( fPairCut->Pass(partpair) ) //check pair cut
- { //do not meets crietria of the pair cut, try with swaped pairs
- if( fPairCut->Pass(partpair->GetSwapedPair() ) )
- continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
- else
- { //swaped pair meets all the criteria
- tmppartpair = partpair->GetSwapedPair();
- }
- }
- else
- {
- tmppartpair = partpair;
- }
-
- UInt_t ii;
- for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
- }
- }
- }
+ //accepted by any cut
+ // we have to copy because reader keeps only one event
+ partEvent1->AddParticle(new AliHBTParticle(*part1));
+ }
- /***************************************/
- /***** Filling diff histogram *********/
- /***************************************/
- if (fBufferSize != 0) //less then 0 mix everything, == 0 do not mix denominator
- for (Int_t i = 0;i<nev-1;i++) //In each event (but last)....
- {
- if ( fNParticleFunctions ==0 )
- continue;
+ if (firstcut) continue;
- partEvent= fReader->GetParticleEvent(i);
- if (!partEvent) continue;
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+ fParticleMonitorFunctions[ii]->Process(part1);
- for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
- {
- part1= partEvent->GetParticle(j);
- if (fPairCut->GetFirstPartCut()->Pass(part1)) continue;
- Int_t maxeventnumber;
+ if ( fNParticleFunctions == 0 ) continue;
+
+ for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+ {
+ part2= partEvent->GetParticle(k);
+ if (part1->GetUID() == part2->GetUID()) continue;
- if ( ((i+fBufferSize) >= nev) ||( fBufferSize < 0) ) //if buffer size is negative
- //or current event+buffersize is greater
- //than max nuber of events
- {
- maxeventnumber = nev; //set the max event number
+ partpair->SetParticles(part1,part2);
+ if(fPairCut->Pass(partpair)) //check pair cut
+ { //do not meets crietria of the
+ if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
+ else tmppartpair = partpair->GetSwapedPair();
}
- else
+ else
{
- maxeventnumber = i+fBufferSize; //set the current event number + fBufferSize
+ tmppartpair = partpair;
}
-
- for (Int_t k = i+1; k<maxeventnumber;k++) // ... Loop over next event
- {
-
- partEvent2= fReader->GetParticleEvent(k);
- if (!partEvent2) continue;
-
- if ( (j%fDisplayMixingInfo) == 0)
- Info("ProcessParticles",
- "Mixing particle %d from event %d with particles from event %d",j,i,k);
-
- for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
+ }
+ /***************************************/
+ /***** Filling denominators *********/
+ /***************************************/
+
+ if (fBufferSize == 0) continue;
+
+ partbuffer.ResetIter();
+
+ while (( partEvent2 = partbuffer.Next() ))
+ {
+ for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
{
+
part2= partEvent2->GetParticle(l);
partpair->SetParticles(part1,part2);
-
- if(fPairCut->Pass(partpair)) //check pair cut
+
+ if( fPairCut->Pass(partpair) ) //check pair cut
{ //do not meets crietria of the
- if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
- else tmppartpair = partpair->GetSwapedPair();
+ if( fPairCut->Pass(partpair->GetSwapedPair()) )
+ continue;
+ else
+ {
+ tmppartpair = partpair->GetSwapedPair();
+ }
}
else
- {
- tmppartpair = partpair;
- }
- UInt_t ii;
+ {//meets criteria of the pair cut
+ tmppartpair = partpair;
+ }
+
for(ii = 0;ii<fNParticleFunctions;ii++)
fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
-
- }//for(Int_t l = 0; l<N2;l++) // ... on all particles
- }//for (Int_t k = i+1; k<maxeventnumber;k++) // ... Loop over next event
+
+ }//for(Int_t l = 0; l<N2;l++) // ... on all particles
+ }
}
- }
- /***************************************/
-
+ partEvent1 = partbuffer.Push(partEvent1);
+ }//while (fReader->Next() == kFALSE)
}
/*************************************************************************************/
AliHBTEvent * trackEvent2=0x0,*partEvent2=0x0;
AliHBTEvent * trackEvent3=0x0,*partEvent3=0x0;
- AliHBTEvent * rawtrackEvent, *rawpartEvent;
+ AliHBTEvent * rawtrackEvent, *rawpartEvent;//this we get from Reader
- Int_t nev = fReader->GetNumberOfTrackEvents();
-
AliHBTPair * trackpair = new AliHBTPair();
AliHBTPair * partpair = new AliHBTPair();
- TList tbuffer;
- TList pbuffer;
- Int_t ninbuffer = 0;
- UInt_t ii;
+ AliHBTEventBuffer partbuffer(fBufferSize);
+ AliHBTEventBuffer trackbuffer(fBufferSize);
+
+ register UInt_t ii;
trackEvent1 = new AliHBTEvent();
partEvent1 = new AliHBTEvent();
trackEvent1->SetOwner(kFALSE);
partEvent1->SetOwner(kFALSE);;
+ Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
+
+ fReader->Rewind();
+
Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
Info("ProcessTracksAndParticlesNonIdentAnal","***** NON IDENT MODE ****************");
Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
- for (Int_t i = 0;i<nev;i++)
+ for (Int_t i = 0;;i++)//infinite loop
{
- rawpartEvent = fReader->GetParticleEvent(i);
- rawtrackEvent = fReader->GetTrackEvent(i);
+ if (fReader->Next()) break; //end when no more events available
+
+ rawpartEvent = fReader->GetParticleEvent();
+ rawtrackEvent = fReader->GetTrackEvent();
if ( (rawpartEvent == 0x0) || (rawtrackEvent == 0x0) ) continue;//in case of any error
-
+
+ if ( rawpartEvent->GetNumberOfParticles() != rawtrackEvent->GetNumberOfParticles() )
+ {
+ Fatal("ProcessTracksAndParticlesNonIdentAnal",
+ "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
+ i,rawpartEvent->GetNumberOfParticles() , rawtrackEvent->GetNumberOfParticles());
+ return;
+ }
+
/********************************/
/* Filtering out */
/********************************/
- if ( (fBufferSize != 0) || ( (partEvent2==0x0)||(trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
- if ((ninbuffer > fBufferSize) && (fBufferSize > 0))
- {//if we have in buffer desired number of events, use the last. If fBufferSize<0 mix all events for background
- partEvent2 = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
- trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
- ninbuffer--;
- }
- else
- {
+ if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
+ {
partEvent2 = new AliHBTEvent();
trackEvent2 = new AliHBTEvent();
- partEvent2->SetOwner(kFALSE);
- trackEvent2->SetOwner(kFALSE);
- }
+ partEvent2->SetOwner(kTRUE);
+ trackEvent2->SetOwner(kTRUE);
+ }
+
FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
part1= partEvent1->GetParticle(j);
track1= trackEvent1->GetParticle(j);
-
+
+ if( part1->GetPdgCode() != track1->GetPdgCode() )
+ {
+ Fatal("ProcessTracksAndParticlesNonIdentAnal",
+ "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
+ i,j, part1->GetPdgCode(),track1->GetPdgCode() );
+ return;
+ }
+
for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
fParticleMonitorFunctions[ii]->Process(part1);
for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
+ if (nocorrfctns) continue;
+
/***************************************/
/****** filling numerators ********/
/****** (particles from event2) ********/
/***************************************/
- for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
+
+ for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
{
part2= partEvent2->GetParticle(k);
if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
/***************************************/
/***** Filling denominators *********/
/***************************************/
- TIter piter(&pbuffer);
- TIter titer(&tbuffer);
+ partbuffer.ResetIter();
+ trackbuffer.ResetIter();
+
Int_t nmonitor = 0;
- while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
+ while ( (partEvent3 = partbuffer.Next() ) != 0x0)
{
- trackEvent3 = (AliHBTEvent*)titer.Next();
+ trackEvent3 = trackbuffer.Next();
+
if ( (j%fDisplayMixingInfo) == 0)
Info("ProcessTracksAndParticlesNonIdentAnal",
"Mixing particle %d from event %d with particles from event %d",j,i,i-(++nmonitor));
}//while event2
}//for over particles in event1
- if ( fBufferSize == 0) continue;//do not mix diff histograms-> do not add to buffer list
- pbuffer.AddFirst(partEvent2);
- tbuffer.AddFirst(trackEvent2);
- partEvent2 = 0x0;
- trackEvent2 = 0x0;
- ninbuffer++;
+ partEvent2 = partbuffer.Push(partEvent2);
+ trackEvent2 = trackbuffer.Push(trackEvent2);
}//end of loop over events (1)
- pbuffer.SetOwner(); //to clean stored events by the way of buffer destruction
- tbuffer.SetOwner();
+
+ partbuffer.SetOwner(kTRUE);
+ trackbuffer.SetOwner(kTRUE);
+
delete partEvent1;
delete trackEvent1;
+ delete partEvent2;
+ delete trackEvent2;
delete partpair;
delete trackpair;
- if ( fBufferSize == 0)
- {//in that case we did not add these events to list
- delete partEvent2;
- delete trackEvent2;
- }
}
/*************************************************************************************/
AliHBTEvent * trackEvent2=0x0;
AliHBTEvent * trackEvent3=0x0;
-
AliHBTEvent * rawtrackEvent;
- Int_t nev = fReader->GetNumberOfTrackEvents();
-
AliHBTPair * trackpair = new AliHBTPair();
+ AliHBTEventBuffer trackbuffer(fBufferSize);
- TList tbuffer;
- Int_t ninbuffer = 0;
register UInt_t ii;
trackEvent1 = new AliHBTEvent();
trackEvent1->SetOwner(kFALSE);
+ fReader->Rewind();
+
Info("ProcessTracksNonIdentAnal","**************************************");
Info("ProcessTracksNonIdentAnal","***** NON IDENT MODE ****************");
Info("ProcessTracksNonIdentAnal","**************************************");
- for (Int_t i = 0;i<nev;i++)
+
+ for (Int_t i = 0;;i++)//infinite loop
{
- rawtrackEvent = fReader->GetTrackEvent(i);
+ if (fReader->Next()) break; //end when no more events available
+ rawtrackEvent = fReader->GetTrackEvent();
+
if (rawtrackEvent == 0x0) continue;//in case of any error
/********************************/
/* Filtering out */
/********************************/
- if ( (fBufferSize != 0) || ( trackEvent2==0x0) )//in case fBufferSize == 0 and pointers are created do not eneter
- if ((ninbuffer > fBufferSize) && (fBufferSize > 0))
- {//if we have in buffer desired number of events, use the last. If fBufferSize<0 mix all events for background
- trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
- ninbuffer--;
- }
- else
- {
+ if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
+ {
trackEvent2 = new AliHBTEvent();
- trackEvent2->SetOwner(kFALSE);
- }
+ trackEvent2->SetOwner(kTRUE);
+ }
+
FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
fTrackMonitorFunctions[ii]->Process(track1);
-
+
+ if (fNTrackFunctions == 0x0) continue;
+
/***************************************/
/****** filling numerators ********/
/****** (particles from event2) ********/
/***************************************/
/***** Filling denominators *********/
/***************************************/
- TIter titer(&tbuffer);
+ trackbuffer.ResetIter();
Int_t nmonitor = 0;
- while ( (trackEvent3 = (AliHBTEvent*)titer.Next()) != 0x0)
+ while ( (trackEvent3 = trackbuffer.Next() ) != 0x0)
{
if ( (j%fDisplayMixingInfo) == 0)
}//while event2
}//for over particles in event1
- if ( fBufferSize == 0) continue;//do not mix diff histograms
- tbuffer.AddFirst(trackEvent2);
- trackEvent2 = 0x0;
- ninbuffer++;
+ trackEvent2 = trackbuffer.Push(trackEvent2);
}//end of loop over events (1)
+ trackbuffer.SetOwner(kTRUE);
delete trackpair;
delete trackEvent1;
- if (fBufferSize == 0) delete trackEvent2;
- tbuffer.SetOwner();
+ delete trackEvent2;
}
/*************************************************************************************/
AliHBTEvent * rawpartEvent = 0x0;
- Int_t nev = fReader->GetNumberOfPartEvents();
-
AliHBTPair * partpair = new AliHBTPair();
+ AliHBTEventBuffer partbuffer(fBufferSize);
- TList pbuffer;
- Int_t ninbuffer = 0;
+ register UInt_t ii;
partEvent1 = new AliHBTEvent();
- partEvent1->SetOwner(kFALSE);;
+ partEvent1->SetOwner(kFALSE);
+
+ fReader->Rewind();
Info("ProcessParticlesNonIdentAnal","**************************************");
Info("ProcessParticlesNonIdentAnal","***** NON IDENT MODE ****************");
Info("ProcessParticlesNonIdentAnal","**************************************");
- for (Int_t i = 0;i<nev;i++)
+ for (Int_t i = 0;;i++)//infinite loop
{
- rawpartEvent = fReader->GetParticleEvent(i);
+ if (fReader->Next()) break; //end when no more events available
+
+ rawpartEvent = fReader->GetParticleEvent();
if ( rawpartEvent == 0x0 ) continue;//in case of any error
/********************************/
/* Filtering out */
/********************************/
- if ( (fBufferSize != 0) || ( partEvent2==0x0) )//in case fBufferSize == 0 and pointers are created do not eneter
- if ((ninbuffer > fBufferSize) && (fBufferSize > 0))
- {//if we have in buffer desired number of events, use the last. If fBufferSize<0 mix all events for background
- partEvent2 = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
- ninbuffer--;
- }
- else
- {
+ if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
+ {
partEvent2 = new AliHBTEvent();
- partEvent2->SetOwner(kFALSE);
- }
+ partEvent2->SetOwner(kTRUE);
+ }
+
FilterOut(partEvent1, partEvent2, rawpartEvent);
for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
part1= partEvent1->GetParticle(j);
- UInt_t zz;
- for(zz = 0; zz<fNParticleMonitorFunctions; zz++)
- fParticleMonitorFunctions[zz]->Process(part1);
-
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+ fParticleMonitorFunctions[ii]->Process(part1);
+
+ if (fNParticleFunctions == 0) continue;
+
/***************************************/
/****** filling numerators ********/
/****** (particles from event2) ********/
}
else
{//meets criteria of the pair cut
- UInt_t ii;
for(ii = 0;ii<fNParticleFunctions;ii++)
- {
fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
- }
}
}
if ( fBufferSize == 0) continue;//do not mix diff histograms
/***************************************/
/***** Filling denominators *********/
/***************************************/
- TIter piter(&pbuffer);
+ partbuffer.ResetIter();
Int_t nmonitor = 0;
- while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
+ while ( (partEvent3 = partbuffer.Next() ) != 0x0)
{
if ( (j%fDisplayMixingInfo) == 0)
Info("ProcessParticlesNonIdentAnal",
}
else
{//meets criteria of the pair cut
- UInt_t ii;
for(ii = 0;ii<fNParticleFunctions;ii++)
{
fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
}// for particles event2
}//while event2
}//for over particles in event1
- if ( fBufferSize == 0) continue;//do not mix diff histograms
- pbuffer.AddFirst(partEvent2);
- partEvent2 = 0x0;
- ninbuffer++;
-
+ partEvent2 = partbuffer.Push(partEvent2);
}//end of loop over events (1)
+
+ partbuffer.SetOwner(kTRUE);
delete partpair;
delete partEvent1;
- if ( fBufferSize == 0) delete partEvent2;
- pbuffer.SetOwner();//to delete stered events.
+ delete partEvent2;
}
/*************************************************************************************/
if (in2)
{
- outpart2->AddParticle(part);
- outtrack2->AddParticle(track);
+ outpart2->AddParticle(new AliHBTParticle(*part));
+ outtrack2->AddParticle(new AliHBTParticle(*track));
continue;
}
}
}
}
/**************************************************************************/
+
+AliHBTEvent::AliHBTEvent(const AliHBTEvent& source):
+ TObject(source),
+ fSize(source.fSize),
+ fParticles(new AliHBTParticle* [fSize]),
+ fNParticles(source.fNParticles),
+ fOwner(source.fNParticles),
+ fRandomized(source.fRandomized)
+{
+//copy constructor
+ for(Int_t i =0; i<fNParticles; i++)
+ {
+ fParticles[i] = new AliHBTParticle( *(source.fParticles[i]) );
+ }
+
+}
+/**************************************************************************/
+AliHBTEvent& AliHBTEvent::operator=(const AliHBTEvent source)
+{
+ // assigment operator
+ Reset();
+ if (fParticles) delete [] fParticles;
+ fSize = source.fSize;
+ fParticles = new AliHBTParticle* [fSize];
+ fNParticles = source.fNParticles;
+ fOwner = source.fNParticles;
+ fRandomized = source.fRandomized;
+
+ for(Int_t i =0; i<fNParticles; i++)
+ {
+ fParticles[i] = new AliHBTParticle( *(source.fParticles[i]) );
+ }
+ return *this;
+}
+/**************************************************************************/
AliHBTEvent::~AliHBTEvent()
{
//destructor
if(fParticles && fOwner)
{
for(Int_t i =0; i<fNParticles; i++)
- delete fParticles[i];
+ {
+ for (Int_t j = i+1; j<fNParticles; j++)
+ if (fParticles[j] == fParticles[i]) fParticles[j] = 0x0;
+ delete fParticles[i];
+ }
}
fNParticles = 0;
}
{
public:
AliHBTEvent();
+ AliHBTEvent(const AliHBTEvent& source);
virtual ~AliHBTEvent();
+
+ AliHBTEvent& operator=(const AliHBTEvent source);
+
static const UInt_t fgkInitEventSize; //initial number of the array
//if expanded, this size is used also
AliHBTParticle* GetParticle(Int_t n); //gets particle
--- /dev/null
+#include "AliHBTEventBuffer.h"
+
+ClassImp(AliHBTEventBuffer)
+
+//______________________________________________________
+////////////////////////////////////////////////////////
+//
+// class AliHBTEventBuffer
+//
+// FIFO type event buffer
+
+AliHBTEventBuffer::AliHBTEventBuffer():
+ fSize(-1),fEvents(),fIter(&fEvents)
+{
+ //ctor
+}
+/***********************************************************/
+AliHBTEventBuffer::AliHBTEventBuffer(Int_t size):
+ fSize(size),fEvents(),fIter(&fEvents)
+{
+ //ctor
+}
+
+AliHBTEvent* AliHBTEventBuffer::Push(AliHBTEvent* event)
+{
+ //adds a new event, and returns old of do not fit in size
+ if (fSize == 0) return event;
+
+ AliHBTEvent* ret = 0x0;
+
+ if (fSize == fEvents.GetSize())
+ ret = dynamic_cast<AliHBTEvent*>(fEvents.Remove(fEvents.Last()));
+ if (event) fEvents.AddFirst(event);
+ return ret;
+}
+
--- /dev/null
+#ifndef ALIHBTEVENTBUFFER_H
+#define ALIHBTEVENTBUFFER_H
+
+#include <TObject.h>
+#include <TList.h>
+#include "AliHBTEvent.h"
+
+class AliHBTEventBuffer: public TObject
+{
+ public:
+ AliHBTEventBuffer();
+ AliHBTEventBuffer(Int_t size);
+ virtual ~AliHBTEventBuffer(){}
+
+ AliHBTEvent* Push(AliHBTEvent* event);//adds a new event, and returns old of do not fit in size
+ AliHBTEvent* RemoveLast(){return dynamic_cast<AliHBTEvent*>(fEvents.Remove(fEvents.Last()));}
+ void ResetIter(){fIter.Reset();}
+ AliHBTEvent* Next(){return dynamic_cast<AliHBTEvent*>( fIter.Next() );}
+ void SetSize(Int_t size){fSize = size;}
+ Int_t GetSize() const {return fSize;}
+ void SetOwner(Bool_t flag) {fEvents.SetOwner(flag);}
+ protected:
+ private:
+ Int_t fSize;//size of buffer; if 0 infinite size
+ TList fEvents;//list with arrays
+ TIter fIter;//iterator
+ ClassDef(AliHBTEventBuffer,1)
+};
+
+
+#endif
//denominator histogram
fNumerator = new TH3D(numstr.Data(),numstr.Data(),
- nxbins,xmin,xmin,nybins,ymin,ymin,nzbins,zmin,zmin);
+ nxbins,xmin,xmax,nybins,ymin,ymax,nzbins,zmin,zmax);
fDenominator = new TH3D(denstr.Data(),denstr.Data(),
- nxbins,xmin,xmin,nybins,ymin,ymin,nzbins,zmin,zmin);
+ nxbins,xmin,xmax,nybins,ymin,ymax,nzbins,zmin,zmax);
fNumerator->Sumw2();
fDenominator->Sumw2();
}
/**********************************************************/
-AliHBTPairCut::AliHBTPairCut(const AliHBTPairCut& in)
+AliHBTPairCut::AliHBTPairCut(const AliHBTPairCut& in):
+ TNamed(in)
{
//copy constructor
fCuts = new AliHbtBasePairCut*[fgkMaxCuts];
}
//______________________________________________________________________________
AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
+ TObject(in),
fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
fCalcMass(in.GetCalcMass()),
}
//______________________________________________________________________________
+AliHBTParticle::~AliHBTParticle()
+{
+//dtor
+ delete [] fPids;
+ delete [] fPidProb;
+}
+//______________________________________________________________________________
+
+AliHBTParticle& AliHBTParticle::operator=(const AliHBTParticle& in)
+{
+//assigment operator
+
+ fNPids = in.fNPids;
+ delete [] fPids;
+ delete [] fPidProb;
+ Int_t* fPids = new Int_t[fNPids];
+ Float_t* fPidProb = new Float_t[fNPids];
+ for (Int_t i = 0; i < fNPids;i++)
+ {
+ fPids[i] = in.fPids[i];
+ fPidProb[i] = in.fPidProb[i];
+ }
+
+ fPdgIdx = in.fPdgIdx;
+ fIdxInEvent = in.fIdxInEvent;
+ fCalcMass = in.GetCalcMass();
+ fPx = in.Px();
+ fPy = in.Py();
+ fPz = in.Pz();
+ fE = in.Energy();
+ fVx = in.Vx();
+ fVy = in.Vy();
+ fVz = in.Vz();
+ fVt = in.T();
+
+ return *this;
+}
+//______________________________________________________________________________
+
void AliHBTParticle::SetPdgCode(Int_t pdg,Float_t prob)
{
SetPIDprobability(pdg,prob);
AliHBTParticle(const TParticle& p,Int_t idx);
- virtual ~AliHBTParticle(){};
-
+ virtual ~AliHBTParticle();
+
+ AliHBTParticle& operator=(const AliHBTParticle& in);
+
void SetPIDprobability(Int_t pdg, Float_t prob = 1.0);
Float_t GetPIDprobability(Int_t pdg);
}
/******************************************************************/
-AliHBTParticleCut::AliHBTParticleCut(const AliHBTParticleCut& in)
+AliHBTParticleCut::AliHBTParticleCut(const AliHBTParticleCut& in):
+ TObject()
{
fCuts = new AliHbtBaseCut* [fkgMaxCuts];//last property in the property
//property enum => defines number of properties
}
/******************************************************************/
-Bool_t AliHBTLogicalOperCut::AliHBTDummyBaseCut::Pass(AliHBTParticle*p)
+Bool_t AliHBTLogicalOperCut::AliHBTDummyBaseCut::Pass(AliHBTParticle* /*part*/)
{
Warning("Pass","You are using dummy base cut! Probobly some logical cut is not set up properly");
return kFALSE;//accept
AliHBTLogicalOperCut(AliHbtBaseCut* first, AliHbtBaseCut* second);
virtual ~AliHBTLogicalOperCut();
protected:
- Double_t GetValue(AliHBTParticle * p){MayNotUse("GetValue");return 0.0;}
+ Double_t GetValue(AliHBTParticle * /*part*/){MayNotUse("GetValue");return 0.0;}
AliHbtBaseCut* fFirst; //second cut
AliHbtBaseCut* fSecond; //first cut
private:
class AliHBTDummyBaseCut: public AliHbtBaseCut
{
- Double_t GetValue(AliHBTParticle * p){return 0.0;}
- Bool_t Pass(AliHBTParticle*p);
+ Double_t GetValue(AliHBTParticle * /*part*/){return 0.0;}
+ Bool_t Pass(AliHBTParticle* /*part*/);
};
ClassDef(AliHBTLogicalOperCut,1)
}
/*********************************************************************/
-void AliHBTRndmGaussBall::Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*p)
+void AliHBTRndmGaussBall::Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*/*particle*/)
{
//randomizez gauss for each coordinate separately
x = gRandom->Gaus(0.0,fRx);
// //
///////////////////////////////////////////////////////////////////////
-void AliHBTRndmCyllSurf::Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*p)
+void AliHBTRndmCyllSurf::Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle* particle)
{
Double_t r = fR + gRandom->Gaus(0.0, 1.0);
- Double_t sf = r/p->Pt();//scaling factor for position transformation ->
+ Double_t sf = r/particle->Pt();//scaling factor for position transformation ->
//we move direction of string momentum but legth defined by r
- x = sf*p->Px();
- y = sf*p->Py();
+ x = sf*particle->Px();
+ y = sf*particle->Py();
z = gRandom->Uniform(-fL,fL);
}
AliHBTRndmGaussBall(Float_t r);
AliHBTRndmGaussBall(Float_t rx, Float_t ry, Float_t rz);
~AliHBTRndmGaussBall(){}
- void Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*p);
+ void Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*/*particle*/);
private:
Float_t fRx;
Float_t fRy;
AliHBTRndmCyllSurf(Float_t r, Float_t l):fR(r),fL(l){}
~AliHBTRndmCyllSurf(){}
- void Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle*p);
+ void Randomize(Double_t& x,Double_t& y,Double_t&z, AliHBTParticle* particle);
private:
Float_t fR;
Float_t fL;
#include <TObjString.h>
#include <TObjArray.h>
#include <TClass.h>
-#include <Riostream.h>
#include "AliHBTParticleCut.h"
-
+#include "AliHBTEvent.h"
+#include "AliHBTRun.h"
ClassImp(AliHBTReader)
//pure virtual
-
+
/*************************************************************************************/
-AliHBTReader::AliHBTReader()
+AliHBTReader::AliHBTReader():
+ fCuts(new TObjArray()),
+ fDirs(0x0),
+ fCurrentEvent(0),
+ fCurrentDir(0),
+ fNEventsRead(0),
+ fTracksEvent(0x0),
+ fParticlesEvent(0x0),
+ fParticles(0x0),
+ fTracks(0x0),
+ fIsRead(kFALSE),
+ fBufferEvents(kFALSE)
{
//constructor
- fCuts = new TObjArray();
- fDirs = 0x0;
}
+/*************************************************************************************/
+AliHBTReader::AliHBTReader(TObjArray* dirs):
+ fCuts(new TObjArray()),
+ fDirs(dirs),
+ fCurrentEvent(0),
+ fCurrentDir(0),
+ fNEventsRead(0),
+ fTracksEvent(0x0),
+ fParticlesEvent(0x0),
+ fParticles(0x0),
+ fTracks(0x0),
+ fIsRead(kFALSE),
+ fBufferEvents(kFALSE)
+{
+//ctor with array of directories to read as parameter
+}
/*************************************************************************************/
-AliHBTReader::AliHBTReader(TObjArray* dirs)
- {
- fCuts = new TObjArray();
- fDirs = dirs;
- }
AliHBTReader::~AliHBTReader()
{
fCuts->SetOwner();
delete fCuts;
}
+ delete fParticlesEvent;
+ delete fTracksEvent;
}
+/*************************************************************************************/
+Int_t AliHBTReader::Next()
+{
+//moves to next event
+ if ( ReadNext() == kTRUE)
+ return kTRUE;
+
+ if (fBufferEvents)
+ {
+ if ( ReadsTracks() && fTracksEvent)
+ fTracks->SetEvent(fNEventsRead-1,fTracksEvent);
+ if ( ReadsParticles() && fParticlesEvent)
+ fParticles->SetEvent(fNEventsRead-1,fParticlesEvent);
+ }
+ return kFALSE;
+}
/*************************************************************************************/
void AliHBTReader::AddParticleCut(AliHBTParticleCut* cut)
AliHBTParticleCut *c = (AliHBTParticleCut*)cut->Clone();
fCuts->Add(c);
}
+/********************************************************************/
+
+AliHBTEvent* AliHBTReader::GetParticleEvent(Int_t n)
+ {
+ //returns Nth event with simulated particles
+ if (ReadsParticles() == kFALSE)
+ {
+ Error("GetParticleEvent","This reader is not able to provide simulated particles.");
+ return 0;
+ }
+
+ if (!fIsRead)
+ {
+ if (ReadsParticles() && (fParticles == 0x0)) fParticles = new AliHBTRun();
+ if (ReadsTracks() && (fTracks == 0x0)) fTracks = new AliHBTRun();
+
+ if (Read(fParticles,fTracks))
+ {
+ Error("GetParticleEvent","Error in reading");
+ return 0x0;
+ }
+ else fIsRead = kTRUE;
+ }
+ return fParticles->GetEvent(n);
+ }
+/********************************************************************/
+
+AliHBTEvent* AliHBTReader::GetTrackEvent(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+ if (ReadsTracks() == kFALSE)
+ {
+ Error("GetTrackEvent","This reader is not able to provide recosntructed tracks.");
+ return 0;
+ }
+ if (!fIsRead)
+ {
+ if (ReadsParticles() && (fParticles == 0x0)) fParticles = new AliHBTRun();
+ if (ReadsTracks() && (fTracks == 0x0)) fTracks = new AliHBTRun();
+
+ if(Read(fParticles,fTracks))
+ {
+ Error("GetTrackEvent","Error in reading");
+ return 0x0;
+ }
+ else fIsRead = kTRUE;
+ }
+ return fTracks->GetEvent(n);
+ }
+/********************************************************************/
+Int_t AliHBTReader::GetNumberOfPartEvents()
+ {
+ //returns number of events of particles
+ if (ReadsParticles() == kFALSE)
+ {
+ Error("GetNumberOfPartEvents","This reader is not able to provide simulated particles.");
+ return 0;
+ }
+
+ if (!fIsRead)
+ {
+ if (ReadsParticles() && (fParticles == 0x0)) fParticles = new AliHBTRun();
+ if (ReadsTracks() && (fTracks == 0x0)) fTracks = new AliHBTRun();
+
+ if (Read(fParticles,fTracks))
+ {
+ Error("GetNumberOfPartEvents","Error in reading");
+ return 0;
+ }
+ else fIsRead = kTRUE;
+ }
+ return fParticles->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+Int_t AliHBTReader::GetNumberOfTrackEvents()
+ {
+ //returns number of events of tracks
+ if (ReadsTracks() == kFALSE)
+ {
+ Error("GetNumberOfTrackEvents","This reader is not able to provide recosntructed tracks.");
+ return 0;
+ }
+ if (!fIsRead)
+ {
+ if (ReadsParticles() && (fParticles == 0x0)) fParticles = new AliHBTRun();
+ if (ReadsTracks() && (fTracks == 0x0)) fTracks = new AliHBTRun();
+
+ if(Read(fParticles,fTracks))
+ {
+ Error("GetNumberOfTrackEvents","Error in reading");
+ return 0;
+ }
+ else fIsRead = kTRUE;
+ }
+ return fTracks->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+Int_t AliHBTReader::Read(AliHBTRun* particles, AliHBTRun *tracks)
+{
+ //reads data and puts put to the particles and tracks objects
+ //reurns 0 if everything is OK
+ //
+ Info("Read","");
+
+ if (!particles) //check if an object is instatiated
+ {
+ Error("Read"," particles object must be instatiated before passing it to the reader");
+ return 1;
+ }
+ if (!tracks) //check if an object is instatiated
+ {
+ Error("Read"," tracks object must be instatiated before passing it to the reader");
+ return 1;
+ }
+ particles->Reset();//clear runs == delete all old events
+ tracks->Reset();
+
+ Rewind();
+
+ Int_t i = 0;
+ while(Next() == kFALSE)
+ {
+ if (ReadsTracks()) tracks->SetEvent(i,fTracksEvent);
+ if (ReadsParticles()) particles->SetEvent(i,fParticlesEvent);
+ i++;
+ }
+ return 0;
+}
/*************************************************************************************/
Bool_t AliHBTReader::Pass(AliHBTParticle* p)
/*************************************************************************************/
TString& AliHBTReader::GetDirName(Int_t entry)
- {
- TString* retval;//return value
- if (fDirs == 0x0)
- {
- retval = new TString(".");
- return *retval;
- }
-
- if ( (entry>fDirs->GetEntries()) || (entry<0))//if out of bounds return empty string
- { //note that entry==0 is accepted even if array is empty (size=0)
- Error("GetDirName","Name out of bounds");
- retval = new TString();
- return *retval;
- }
-
- if (fDirs->GetEntries() == 0)
- {
- retval = new TString(".");
- return *retval;
- }
-
- TClass *objclass = fDirs->At(entry)->IsA();
- TClass *stringclass = TObjString::Class();
-
- TObjString *dir = (TObjString*)objclass->DynamicCast(stringclass,fDirs->At(entry));
-
- if(dir == 0x0)
- {
- Error("GetDirName","Object in TObjArray is not a TObjString or its descendant");
- retval = new TString();
- return *retval;
- }
- if (gDebug > 0) cout<<"Returned ok "<<dir->String().Data()<<endl;
- return dir->String();
- }
+{
+//returns directory name of next one to read
+ TString* retval;//return value
+ if (fDirs == 0x0)
+ {
+ retval = new TString(".");
+ return *retval;
+ }
+
+ if ( (entry>fDirs->GetEntries()) || (entry<0))//if out of bounds return empty string
+ { //note that entry==0 is accepted even if array is empty (size=0)
+ Error("GetDirName","Name out of bounds");
+ retval = new TString();
+ return *retval;
+ }
+
+ if (fDirs->GetEntries() == 0)
+ {
+ retval = new TString(".");
+ return *retval;
+ }
+
+ TClass *objclass = fDirs->At(entry)->IsA();
+ TClass *stringclass = TObjString::Class();
+
+ TObjString *dir = (TObjString*)objclass->DynamicCast(stringclass,fDirs->At(entry));
+
+ if(dir == 0x0)
+ {
+ Error("GetDirName","Object in TObjArray is not a TObjString or its descendant");
+ retval = new TString();
+ return *retval;
+ }
+ if (gDebug > 0) Info("GetDirName","Returned ok %s",dir->String().Data());
+ return dir->String();
+}
#define ALIHBTREADER_H
#include <TNamed.h>
+#include <TObjArray.h>
//Reader Base class (reads particles and tracks and
//puts it to the AliHBTRun objects
class AliHBTRun;
class AliHBTEvent;
class AliHBTParticleCut;
-class TObjArray;
class AliHBTParticle;
class TString;
AliHBTReader();
AliHBTReader(TObjArray*);
virtual ~AliHBTReader();
- //in the future this class is will read global tracking
- virtual Int_t Read(AliHBTRun* particles, AliHBTRun *tracks) = 0;
- virtual AliHBTEvent* GetParticleEvent(Int_t) = 0;
- virtual AliHBTEvent* GetTrackEvent(Int_t) = 0;
- virtual Int_t GetNumberOfPartEvents() = 0;
- virtual Int_t GetNumberOfTrackEvents() = 0;
+ virtual Int_t Next();
+ virtual void Rewind() = 0;
- void AddParticleCut(AliHBTParticleCut* cut);
+ virtual Bool_t ReadsTracks() const = 0;
+ virtual Bool_t ReadsParticles() const = 0;
- void SetDirs(TObjArray* dirs){fDirs = dirs;} //sets array directories names
-
- virtual AliHBTEvent* GetNextParticleEvent(){return 0x0;}
- virtual AliHBTEvent* GetNextTrackEvent(){return 0x0;}
-
+ void AddParticleCut(AliHBTParticleCut* cut);
+
+ virtual Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);
+
+ virtual AliHBTEvent* GetParticleEvent() const {return fParticlesEvent;}
+ virtual AliHBTEvent* GetTrackEvent() const {return fTracksEvent;}
+
+ virtual AliHBTEvent* GetParticleEvent(Int_t);
+ virtual AliHBTEvent* GetTrackEvent(Int_t);
+
+ virtual Int_t GetNumberOfPartEvents();
+ virtual Int_t GetNumberOfTrackEvents();
+
+ void SetDirs(TObjArray* dirs){fDirs = dirs;} //sets array directories names
+ void SetEventBuffering(Bool_t flag){fBufferEvents = flag;}
+
+ virtual Int_t GetNumberOfDirs() const {return (fDirs)?fDirs->GetEntries():0;}
protected:
- TObjArray *fCuts;//array with particle cuts
- TObjArray *fDirs;//arry with directories to read data from
-
+ TObjArray* fCuts;//array with particle cuts
+ TObjArray* fDirs;//arry with directories to read data from
+
+ Int_t fCurrentEvent;//! number of current event in current directory
+ Int_t fCurrentDir;//! number of current directory
+
+ Int_t fNEventsRead;//!total
+
+ AliHBTEvent* fTracksEvent; //! tracks read from current event
+ AliHBTEvent* fParticlesEvent; //! particles read from current event
+
+ AliHBTRun* fParticles; //!simulated particles
+ AliHBTRun* fTracks; //!reconstructed tracks (particles)
+
+ Bool_t fIsRead;//!flag indicating if the data are already read
+ Bool_t fBufferEvents;//flag indicating if the data should be bufferred
+
+ virtual Int_t ReadNext() = 0; //this methods reads next event and put result in fTracksEvent and/or fParticlesEvent
Bool_t Pass(AliHBTParticle*);
Bool_t Pass(Int_t pid);
#include <TFile.h>
#include <TParticle.h>
+#include <AliRunLoader.h>
+#include <AliStack.h>
#include <AliESDtrack.h>
#include <AliESD.h>
#include "AliHBTParticle.h"
#include "AliHBTParticleCut.h"
+AliHBTReaderESD r;
+
ClassImp(AliHBTReaderESD)
-AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename):
- fParticles(new AliHBTRun()),
- fTracks(new AliHBTRun()),
+AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
fESDFileName(esdfilename),
- fIsRead(kFALSE)
+ fGAlFileName(galfilename),
+ fFile(0x0),
+ fRunLoader(0x0),
+ fReadParticles(kFALSE)
{
//cosntructor
if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
}
/********************************************************************/
-AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename):
+AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
AliHBTReader(dirs),
- fParticles(new AliHBTRun()),
- fTracks(new AliHBTRun()),
fESDFileName(esdfilename),
- fIsRead(kFALSE)
+ fGAlFileName(galfilename),
+ fFile(0x0),
+ fRunLoader(0x0),
+ fReadParticles(kFALSE)
{
//cosntructor
if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
AliHBTReaderESD::~AliHBTReaderESD()
{
//desctructor
- delete fParticles;
- delete fTracks;
+ delete fFile;
+ delete fRunLoader;
}
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderESD::GetParticleEvent(Int_t n)
- {
- //returns Nth event with simulated particles
- if (!fIsRead)
- if(Read(fParticles,fTracks))
- {
- Error("GetParticleEvent","Error in reading");
- return 0x0;
- }
- return fParticles->GetEvent(n);
- }
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderESD::GetTrackEvent(Int_t n)
- {
- //returns Nth event with reconstructed tracks
- if (!fIsRead)
- if(Read(fParticles,fTracks))
- {
- Error("GetTrackEvent","Error in reading");
- return 0x0;
- }
- return fTracks->GetEvent(n);
- }
-/********************************************************************/
-
-Int_t AliHBTReaderESD::GetNumberOfPartEvents()
- {
- //returns number of events of particles
- if (!fIsRead)
- if ( Read(fParticles,fTracks))
- {
- Error("GetNumberOfPartEvents","Error in reading");
- return 0;
- }
- return fParticles->GetNumberOfEvents();
- }
-
-/********************************************************************/
-Int_t AliHBTReaderESD::GetNumberOfTrackEvents()
- {
- //returns number of events of tracks
- if (!fIsRead)
- if(Read(fParticles,fTracks))
- {
- Error("GetNumberOfTrackEvents","Error in reading");
- return 0;
- }
- return fTracks->GetNumberOfEvents();
- }
-/********************************************************************/
-
-
-Int_t AliHBTReaderESD::Read(AliHBTRun* particles, AliHBTRun *tracks)
+/**********************************************************/
+Int_t AliHBTReaderESD::ReadNext()
{
- //reads data and puts put to the particles and tracks objects
- //reurns 0 if everything is OK
- //
- Info("Read","");
- Int_t totalNevents = 0;
- Int_t currentdir = 0; //number of current directory name is fDirs array
+//reads next event from fFile
+
+ //****** Tentative particle type "concentrations"
+ static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
Double_t w[kNSpecies];
Double_t mom[3];//momentum
Double_t pos[3];//position
- //****** Tentative particle type "concentrations"
- const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
+ if (AliHBTParticle::fgDebug)
+ Info("ReadNext","Entered");
+
TDatabasePDG* pdgdb = TDatabasePDG::Instance();
if (pdgdb == 0x0)
{
- Error("Read","Can not get PDG Database Instance.");
+ Error("Next","Can not get PDG Database Instance.");
return 1;
}
- if (!particles) //check if an object is instatiated
- {
- Error("Read"," particles object must instatiated before passing it to the reader");
- return 1;
- }
- if (!tracks) //check if an object is instatiated
- {
- Error("Read"," tracks object must instatiated before passing it to the reader");
- return 1;
- }
- particles->Reset();//clear runs == delete all old events
- tracks->Reset();
-
- Int_t Ndirs;
- if (fDirs) //if array with directories is supplied by user
- {
- Ndirs = fDirs->GetEntries(); //get the number if directories
- }
- else
- {
- Ndirs = 0; //if the array is not supplied read only from current directory
- }
-
+
+ if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
+ if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
+
+ fParticlesEvent->Reset();
+ fTracksEvent->Reset();
+
do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
{
- TFile* file = OpenFile(currentdir);
- if (file == 0x0)
- {
- Error("Read","Cannot get File for dir no. %d",currentdir);
- currentdir++;
- continue;
- }
-
- for(Int_t currentEvent = 0; ;currentEvent++)//loop over all events
+ if (fFile == 0x0)
{
- TString esdname;
- esdname+=currentEvent;
- AliESD* esd = dynamic_cast<AliESD*>(file->Get(esdname));
- if (esd == 0x0)
+ fFile = OpenFile(fCurrentDir);//rl is opened here
+ if (fFile == 0x0)
{
- if (AliHBTParticle::fgDebug > 2 )
- {
- Info("Read","Can not find ESD object named %s.",esdname.Data());
- }
- currentdir++;
- break;//we have to assume there is no more ESD objects in the file
- }
-
- Info("Read","Reading Event %d",currentEvent);
+ Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
+ fCurrentDir++;
+ continue;
+ }
+ fCurrentEvent = 0;
+ }
+
+ //try to read
+ TString esdname;
+ esdname+=fCurrentEvent;
+ AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
+ if (esd == 0x0)
+ {
+ if (AliHBTParticle::fgDebug > 2 )
+ {
+ Info("ReadNext","Can not find ESD object named %s.",esdname.Data());
+ }
+ fCurrentDir++;
+ delete fFile;//we have to assume there is no more ESD objects in the fFile
+ delete fRunLoader;
+ fFile = 0x0;
+ continue;
+ }
+ AliStack* stack = 0x0;
+ if (fReadParticles && fRunLoader)
+ {
+ fRunLoader->GetEvent(fCurrentEvent);
+ stack = fRunLoader->Stack();
+ }
+ Info("ReadNext","Reading Event %d",fCurrentEvent);
- Int_t ntr = esd->GetNumberOfTracks();
- Info("Read","Found %d tracks.",ntr);
- for (Int_t i = 0;i<ntr; i++)
- {
- AliESDtrack *esdtrack = esd->GetTrack(i);
- if (esdtrack == 0x0)
- {
- Error("Read","Can not get track %d", i);
- continue;
- }
- if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE)
- {
- if (AliHBTParticle::fgDebug > 2)
- Info("Read","Particle skipped: PID BIT is not set.");
- continue;
- }
-
- esdtrack->GetESDpid(pidtable);
- esdtrack->GetPxPyPz(mom);
- esdtrack->GetXYZ(pos);
- Double_t rc=0.;
+ Int_t ntr = esd->GetNumberOfTracks();
+ Info("ReadNext","Found %d tracks.",ntr);
+ for (Int_t i = 0;i<ntr; i++)
+ {
+ AliESDtrack *esdtrack = esd->GetTrack(i);
+ if (esdtrack == 0x0)
+ {
+ Error("Next","Can not get track %d", i);
+ continue;
+ }
+ if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE)
+ {
+ if (AliHBTParticle::fgDebug > 2)
+ Info("ReadNext","Particle skipped: PID BIT is not set.");
+ continue;
+ }
- //Here we apply Bayes' formula
- for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
- if (rc==0.0)
- {
- if (AliHBTParticle::fgDebug > 2)
- Info("Read","Particle rejected since total bayessian PID probab. is zero.");
- continue;
- }
+ esdtrack->GetESDpid(pidtable);
+ esdtrack->GetPxPyPz(mom);
+ esdtrack->GetXYZ(pos);
+
+ //Particle from kinematics
+ AliHBTParticle* particle = 0;
+ Bool_t keeppart = kFALSE;
+ if ( fReadParticles && stack )
+ {
+ if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track
+ TParticle *p = stack->Particle(esdtrack->GetLabel());
+ if (p==0x0)
+ {
+ Error("ReadNext","Can not find track with such label.");
+ continue;
+ }
+ particle = new AliHBTParticle(*p,i);
+ }
+
+ //Here we apply Bayes' formula
+ Double_t rc=0.;
+ for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
+ if (rc==0.0)
+ {
+ if (AliHBTParticle::fgDebug > 2)
+ Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
+ continue;
+ }
-
for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
+ for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
- if (AliHBTParticle::fgDebug > 4)
- {
- Info("Read","###########################################################################");
- Info("Read","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
- Info("Read","Position: %f %f %f",pos[0],pos[1],pos[2]);
- Info("Read","Radius: %f ,rc: %f",TMath::Hypot(pos[0],pos[1]),rc);
- TString msg("Pid list got from track:");
- for (Int_t s = 0;s<kNSpecies;s++)
- {
- msg+="\n ";
- msg+=s;
- msg+="(";
- msg+=GetSpeciesPdgCode((ESpecies)s);
- msg+="): ";
- msg+=w[s];
- msg+=" (";
- msg+=pidtable[s];
- msg+=")";
- }
- Info("Read","%s",msg.Data());
- }//if (AliHBTParticle::fgDebug>4)
-
- for (Int_t s = 0; s<kNSpecies; s++)
- {
- if (w[s] == 0.0) continue;
+ if (AliHBTParticle::fgDebug > 4)
+ {
+ Info("ReadNext","###########################################################################");
+ Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
+ Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
+ TString msg("Pid list got from track:");
+ for (Int_t s = 0;s<kNSpecies;s++)
+ {
+ msg+="\n ";
+ msg+=s;
+ msg+="(";
+ msg+=GetSpeciesPdgCode((ESpecies)s);
+ msg+="): ";
+ msg+=w[s];
+ msg+=" (";
+ msg+=pidtable[s];
+ msg+=")";
+ }
+ Info("ReadNext","%s",msg.Data());
+ }//if (AliHBTParticle::fgDebug>4)
- Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s);
- if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type
+ for (Int_t s = 0; s<kNSpecies; s++)
+ {
+ if (w[s] == 0.0) continue;
- Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
- Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
+ Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s);
+ if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type
- AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i,
- mom[0], mom[1], mom[2], tEtot,
- pos[0], pos[1], pos[2], 0.);
- //copy probabilitis of other species (if not zero)
- for (Int_t k = 0; k<kNSpecies; k++)
- {
- if (k == s) continue;
- if (w[k] == 0.0) continue;
- track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]);
- }
- if(Pass(track))//check if meets all criteria of any of our cuts
- //if it does not delete it and take next good track
- {
- delete track;
- continue;
- }
- tracks->AddParticle(totalNevents,track);
- if (AliHBTParticle::fgDebug > 4 )
- {
- Info("Read","\n\nAdding Particle with incarnation %d",pdgcode);
- track->Print();
- }
- }//for (Int_t s = 0; s<kNSpecies; s++)
+ Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
+ Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
- }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
- totalNevents++;
- }//for(Int_t currentEvent = 0; ;currentEvent++) -- for over ESD in single file
+ AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i,
+ mom[0], mom[1], mom[2], tEtot,
+ pos[0], pos[1], pos[2], 0.);
+ //copy probabilitis of other species (if not zero)
+ for (Int_t k = 0; k<kNSpecies; k++)
+ {
+ if (k == s) continue;
+ if (w[k] == 0.0) continue;
+ track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]);
+ }
+ if(Pass(track))//check if meets all criteria of any of our cuts
+ //if it does not delete it and take next good track
+ {
+ delete track;
+ continue;
+ }
+
+ fTracksEvent->AddParticle(track);
+ if (particle) fParticlesEvent->AddParticle(particle);
+ keeppart = kTRUE;
+
+ if (AliHBTParticle::fgDebug > 4 )
+ {
+ Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
+ track->Print();
+ }
+ }//for (Int_t s = 0; s<kNSpecies; s++)
+
+ if (keeppart == kFALSE) delete particle;//particle was not stored in event
+
+ }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
+
+ Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
+ fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
+ fNEventsRead,fCurrentEvent,fCurrentDir);
- }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
+ fCurrentEvent++;
+ fNEventsRead++;
+ delete esd;
+ return 0;//success -> read one event
+ }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
+
+ return 1; //no more directories to read
+}
+/**********************************************************/
+void AliHBTReaderESD::Rewind()
+{
+ delete fFile;
+ fFile = 0;
+ delete fRunLoader;
+ fCurrentDir = 0;
+ fNEventsRead = 0;
+ fCurrentEvent++;
- fIsRead = kTRUE;
- return 0;
-}
+}
/**********************************************************/
-
+
TFile* AliHBTReaderESD::OpenFile(Int_t n)
{
-//opens file with kine tree
+//opens fFile with kine tree
const TString& dirname = GetDirName(n);
if (dirname == "")
if ( ret == 0x0)
{
- Error("OpenFiles","Can't open file %s",filename.Data());
+ Error("OpenFiles","Can't open fFile %s",filename.Data());
return 0x0;
}
if (!ret->IsOpen())
{
- Error("OpenFiles","Can't open file %s",filename.Data());
+ Error("OpenFiles","Can't open fFile %s",filename.Data());
return 0x0;
}
+ if (fReadParticles)
+ {
+ fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
+ if (fRunLoader == 0x0)
+ {
+ Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
+ delete ret;
+ return 0x0;
+ }
+ fRunLoader->LoadHeader();
+ if (fRunLoader->LoadKinematics())
+ {
+ Error("Next","Error occured while loading kinematics.");
+ delete fRunLoader;
+ delete ret;
+ return 0x0;
+ }
+ }
return ret;
}
/**********************************************************/
#include <TString.h>
class TFile;
+class AliRunLoader;
class AliHBTReaderESD: public AliHBTReader
{
public:
- AliHBTReaderESD(const Char_t* esdfilename = "AliESDs.root");
+ AliHBTReaderESD(const Char_t* esdfilename = "AliESDs.root", const Char_t* galfilename = "galice.root");
- AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename = "AliESDs.root");
+ AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename = "AliESDs.root", const Char_t* galfilename = "galice.root");
virtual ~AliHBTReaderESD();
- Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
+ void Rewind();
- AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
- AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles
- Int_t GetNumberOfPartEvents();//returns number of particle events
- Int_t GetNumberOfTrackEvents();//returns number of track events
+ void ReadParticles(Bool_t flag){fReadParticles = flag;}
+ Bool_t ReadsTracks() const {return kTRUE;}
+ Bool_t ReadsParticles() const {return fReadParticles;}
enum ESpecies {kESDElectron = 0, kESDMuon, kESDPion, kESDKaon, kESDProton, kNSpecies};
- static Int_t GetSpeciesPdgCode(ESpecies spec);//skowron
- protected:
-
- TFile* OpenFile(Int_t evno);//opens files to be read for given event
- void CloseFiles(TFile*);//close files
+ static Int_t GetSpeciesPdgCode(ESpecies spec);//skowron
- AliHBTRun* fParticles; //!simulated particles
- AliHBTRun* fTracks; //!reconstructed tracks (particles)
-
- TString fESDFileName;//name of the file with tracks
- Bool_t fIsRead;//!flag indicating if the data are already read
+ protected:
+ Int_t ReadNext();
+ TFile* OpenFile(Int_t evno);//opens files to be read for given event
+ void CloseFiles(TFile*);//close files
+
+ TString fESDFileName;//name of the file with tracks
+ TString fGAlFileName;//name of the file with tracks
+ TFile* fFile;//! pointer to current ESD file
+ AliRunLoader* fRunLoader;//!Run Loader
+ Bool_t fReadParticles;//flag indicating wether to read particles from kinematics
private:
- ClassDef(AliHBTReaderESD,1)
+ ClassDef(AliHBTReaderESD,2)
};
-
#include "AliHBTReaderITSv2.h"
-#include <Riostream.h>
-#include <Riostream.h>
+#include <Varargs.h>
+
#include <TString.h>
#include <TObjString.h>
#include <TTree.h>
-#include <TFile.h>
#include <TParticle.h>
#include <AliRun.h>
#include <AliRunLoader.h>
#include <AliLoader.h>
+#include <AliStack.h>
#include <AliMagF.h>
-#include <AliITS.h>
#include <AliITStrackV2.h>
-#include <AliITStrackerV2.h>
-#include <AliITSgeom.h>
#include "AliHBTRun.h"
#include "AliHBTEvent.h"
ClassImp(AliHBTReaderITSv2)
-AliHBTReaderITSv2::AliHBTReaderITSv2():fFileName("galice.root")
+AliHBTReaderITSv2::AliHBTReaderITSv2():
+ fFileName("galice.root"),
+ fRunLoader(0x0),
+ fITSLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
{
//constructor,
//Defaults:
// galicefilename = "galice.root"
- fParticles = 0x0;
- fTracks = 0x0;
- fIsRead = kFALSE;
}
/********************************************************************/
-AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* galicefilename):fFileName(galicefilename)
+AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* galicefilename):
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fITSLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
{
//constructor,
//Defaults:
// galicefilename = "galice.root"
- fParticles = new AliHBTRun();
- fTracks = new AliHBTRun();
- fIsRead = kFALSE;
}
/********************************************************************/
AliHBTReaderITSv2::AliHBTReaderITSv2(TObjArray* dirs, const Char_t* galicefilename):
- AliHBTReader(dirs), fFileName(galicefilename)
+ AliHBTReader(dirs),
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fITSLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
{
//constructor,
//Defaults:
// galicefilename = "galice.root"
- fParticles = new AliHBTRun();
- fTracks = new AliHBTRun();
- fIsRead = kFALSE;
}
/********************************************************************/
-
-AliHBTReaderITSv2::~AliHBTReaderITSv2()
- {
- if (fParticles) delete fParticles;
- if (fTracks) delete fTracks;
- }
-/********************************************************************/
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderITSv2::GetParticleEvent(Int_t n)
- {
- //returns Nth event with simulated particles
- if (!fIsRead)
- {
- if (fParticles == 0x0) fParticles = new AliHBTRun();
- if (fTracks == 0x0) fTracks = new AliHBTRun();
- if(Read(fParticles,fTracks))
- {
- Error("GetParticleEvent","Error in reading");
- return 0x0;
- }
- }
- return (fParticles)?fParticles->GetEvent(n):0x0;
+
+void AliHBTReaderITSv2::Rewind()
+{
+ //rewinds reading
+ delete fRunLoader;
+ fRunLoader = 0x0;
+ fCurrentDir = 0;
+ fNEventsRead= 0;
}
/********************************************************************/
-AliHBTEvent* AliHBTReaderITSv2::GetTrackEvent(Int_t n)
- {
- //returns Nth event with reconstructed tracks
- if (!fIsRead)
- {
- if (fParticles == 0x0) fParticles = new AliHBTRun();
- if (fTracks == 0x0) fTracks = new AliHBTRun();
- if(Read(fParticles,fTracks))
- {
- Error("GetTrackEvent","Error in reading");
- return 0x0;
- }
- }
- return (fTracks)?fTracks->GetEvent(n):0x0;
- }
-/********************************************************************/
-
-Int_t AliHBTReaderITSv2::GetNumberOfPartEvents()
- {
- //returns number of events of particles
- if (!fIsRead)
- {
- if (fParticles == 0x0) fParticles = new AliHBTRun();
- if (fTracks == 0x0) fTracks = new AliHBTRun();
- if(Read(fParticles,fTracks))
- {
- Error("GetNumberOfPartEvents","Error in reading");
- return 0;
- }
- }
- return (fParticles)?fParticles->GetNumberOfEvents():0;
- }
-
-/********************************************************************/
-Int_t AliHBTReaderITSv2::GetNumberOfTrackEvents()
- {
- //returns number of events of tracks
- if (!fIsRead)
- {
- if (fParticles == 0x0) fParticles = new AliHBTRun();
- if (fTracks == 0x0) fTracks = new AliHBTRun();
- if(Read(fParticles,fTracks))
- {
- Error("GetNumberOfTrackEvents","Error in reading");
- return 0;
- }
- }
- return (fTracks)?fTracks->GetNumberOfEvents():0;
- }
-/********************************************************************/
+AliHBTReaderITSv2::~AliHBTReaderITSv2()
+{
+ //dtor
+ delete fRunLoader;
+}
/********************************************************************/
-
-Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
+Int_t AliHBTReaderITSv2::ReadNext()
{
-//reads data
- Int_t Nevents = 0; //number of events found in given directory
- Int_t Ndirs; //number of the directories to be read
- Int_t Ntracks; //number of tracks in current event
- Int_t currentdir = 0; //number of events in the current directory
- Int_t totalNevents = 0; //total number of events read from all directories up to now
+//reads data from next event
+
register Int_t i = 0; //iterator
// AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
- TTree *tracktree; // tree for tracks
+ TTree *tracktree = 0x0; // tree for tracks
+ AliITStrackV2 *iotrack = 0x0;
Double_t xk;
Double_t par[5]; //Kalman track parameters
Float_t phi, lam, pt;//angles and transverse momentum
Int_t label; //label of the current track
-
- AliITStrackV2 *iotrack = 0x0; //buffer track for reading data from tree
-
- if (!particles) //check if an object is instatiated
- {
- Error("Read"," particles object must instatiated before passing it to the reader");
- }
- if (!tracks) //check if an object is instatiated
- {
- Error("Read"," tracks object must instatiated before passing it to the reader");
- }
- particles->Reset();//clear runs == delete all old events
- tracks->Reset();
-
- if (fDirs) //if array with directories is supplied by user
- {
- Ndirs = fDirs->GetEntries(); //get the number if directories
- }
- else
- {
- Ndirs = 0; //if the array is not supplied read only from current directory
- }
-// cout<<"Found "<<Ndirs<<" directory entries"<<endl;
+ if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
+ if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
+
+ fParticlesEvent->Reset();
+ fTracksEvent->Reset();
do //do while is good even if Ndirs==0 (than read from current directory)
{
- TString filename = GetDirName(currentdir);
- if (filename.IsNull())
+ if (fRunLoader == 0x0)
+ if (OpenNextFile()) continue;//directory counter is increased inside in case of error
+
+ if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
{
- Error("Read","Can not get directory name");
- currentdir++;
- continue;
+ //read next directory
+ delete fRunLoader;//close current session
+ fRunLoader = 0x0;//assure pointer is null
+ fCurrentDir++;//go to next dir
+ continue;//directory counter is increased inside in case of error
}
- filename = filename +"/"+ fFileName;
- AliRunLoader* rl = AliRunLoader::Open(filename);
- if( rl == 0x0)
+
+ Info("ReadNext","Reading Event %d",fCurrentEvent);
+
+ fRunLoader->GetEvent(fCurrentEvent);
+
+ tracktree=fITSLoader->TreeT();
+ if (!tracktree)
{
- Error("Read","Exiting due to problems with opening files.");
- currentdir++;
+ Error("ReadNext","Can't get a tree with ITS tracks");
+ fCurrentEvent++;
continue;
}
-
- rl->LoadHeader();
- rl->LoadKinematics();
- rl->LoadgAlice();
- gAlice = rl->GetAliRun();
- AliITS* its = (AliITS*)gAlice->GetModule("ITS");
-
- AliLoader* itsl = rl->GetLoader("ITSLoader");
-
- if ((its == 0x0) || ( itsl== 0x0))
+
+ TBranch *tbranch=tracktree->GetBranch("tracks");
+ if (!tbranch)
{
- Error("Read","Can not found ITS in this run");
- delete rl;
- rl = 0x0;
- currentdir++;
+ Error("ReadNext","Can't get a branch with ITS tracks");
+ fCurrentEvent++;
continue;
}
- Nevents = rl->GetNumberOfEvents();
-
- if (Nevents > 0)//check if tree E exists
+
+ AliStack* stack = fRunLoader->Stack();
+ if (stack == 0x0)
{
- Info("Read","________________________________________________________");
- Info("Read","Found %d event(s) in directory %s",Nevents,GetDirName(currentdir).Data());
- Float_t mf;
- if (fUseMagFFromRun)
- {
- mf = gAlice->Field()->SolenoidField();
- Info("Read","Setting Magnetic Field from run: B=%fT",mf/10.);
- }
- else
- {
- Info("Read","Setting Own Magnetic Field: B=%fT",fMagneticField);
- mf = fMagneticField*10.;
- }
- AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
- if (iotrack == 0x0) iotrack = new AliITStrackV2();
- }
- else
- {//if not return an error
- Error("Read","No events in this run");
- delete rl;
- rl = 0x0;
- currentdir++;
+ Error("ReadNext","Can not get stack for current event",fCurrentEvent);
+ fCurrentEvent++;
continue;
}
+
+ //must be here because on the beginning conv. const. is not set yet
+ if (iotrack == 0x0) iotrack = new AliITStrackV2(); //buffer track for reading data from tree
- AliITSgeom *geom= its->GetITSgeom();
- if (!geom)
+ Int_t ntr = (Int_t)tracktree->GetEntries();
+
+ for (i=0; i < ntr; i++) //loop over all tpc tracks
{
- Error("Read","Can't get the ITS geometry!");
- delete rl;
- rl = 0x0;
- currentdir++;
- continue;
- }
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
- itsl->LoadTracks();
+ label=iotrack->GetLabel();
+ if (label < 0)
+ {
+ continue;
+ }
- for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
- {
- cout<<"Reading Event "<<currentEvent<<endl;
- rl->GetEvent(currentEvent);
- tracktree=itsl->TreeT();
-
- if (!tracktree)
- {
- Error("Read","Can't get a tree with ITS tracks");
- continue;
- }
- TBranch *tbranch=tracktree->GetBranch("tracks");
- Ntracks=(Int_t)tracktree->GetEntries();
+ TParticle *p = stack->Particle(label);
+ if(p == 0x0) continue; //if returned pointer is NULL
+ if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
- Int_t accepted = 0;
- Int_t tpcfault = 0;
- Int_t itsfault = 0;
- for (i=0; i<Ntracks; i++) //loop over all tpc tracks
- {
- if(i%100 == 0)cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<"\r";
-
- tbranch->SetAddress(&iotrack);
- tracktree->GetEvent(i);
+ if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
+ //if not take next partilce
- label=iotrack->GetLabel();
- if (label < 0)
- {
- tpcfault++;
- continue;
- }
+ AliHBTParticle* part = new AliHBTParticle(*p,i);
+ if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
+ //if it does not delete it and take next good track
- TParticle *p = (TParticle*)gAlice->Particle(label);
- if(p == 0x0) continue; //if returned pointer is NULL
- if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
+ iotrack->PropagateTo(3.,0.0028,65.19);
+ iotrack->PropagateToVertex();
- if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
- //if not take next partilce
-
- AliHBTParticle* part = new AliHBTParticle(*p,i);
- if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
- //if it does not delete it and take next good track
+ iotrack->GetExternalParameters(xk,par); //get properties of the track
+ phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
+ if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+ if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+ lam=par[3];
+ pt=1.0/TMath::Abs(par[4]);
- iotrack->PropagateTo(3.,0.0028,65.19);
- iotrack->PropagateToVertex();
-
- iotrack->GetExternalParameters(xk,par); //get properties of the track
- phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
- if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
- if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
- lam=par[3];
- pt=1.0/TMath::Abs(par[4]);
-
- Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
- Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
- Double_t tpz = pt * lam; //track z coordinate of momentum
-
- Double_t mass = p->GetMass();
- Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
-
- AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
- if(Pass(track))//check if meets all criteria of any of our cuts
- //if it does not delete it and take next good track
- {
- delete track;
- delete part;
- continue;
- }
- particles->AddParticle(totalNevents,part);//put track and particle on the run
- tracks->AddParticle(totalNevents,track);
- accepted++;
- }//end of loop over tracks in the event
+ Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
+ Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
+ Double_t tpz = pt * lam; //track z coordinate of momentum
+
+ Double_t mass = p->GetMass();
+ Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
+
+ AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+ if(Pass(track))//check if meets all criteria of any of our cuts
+ //if it does not delete it and take next good track
+ {
+ delete track;
+ delete part;
+ continue;
+ }
+
+ fParticlesEvent->AddParticle(part);
+ fTracksEvent->AddParticle(track);
+ }//end of loop over tracks in the event
- totalNevents++;
- cout<<"all: "<<i<<" accepted: "<<accepted<<" tpc faults: "<<tpcfault<<" its faults: "<<itsfault<<endl;
+ Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
+ fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
+ fNEventsRead,fCurrentEvent,fCurrentDir);
- }//end of loop over events in current directory
- delete rl;
- currentdir++;
- }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
+ fCurrentEvent++;
+ fNEventsRead++;
+ delete iotrack;
+ return 0;
+ }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
delete iotrack;
- fIsRead = kTRUE;
- return 0;
+ return 1;
}
/********************************************************************/
+Int_t AliHBTReaderITSv2::OpenNextFile()
+{
+ //opens next file
+ TString filename = GetDirName(fCurrentDir);
+ if (filename.IsNull())
+ {
+ DoOpenError("Can not get directory name");
+ return 1;
+ }
+ filename = filename +"/"+ fFileName;
+ fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
+ if( fRunLoader == 0x0)
+ {
+ DoOpenError("Can not open session.");
+ return 1;
+ }
+
+ if (fRunLoader->GetNumberOfEvents() <= 0)
+ {
+ DoOpenError("There is no events in this directory.");
+ return 1;
+ }
+
+ if (fRunLoader->LoadKinematics())
+ {
+ DoOpenError("Error occured while loading kinematics.");
+ return 1;
+ }
+ fITSLoader = fRunLoader->GetLoader("ITSLoader");
+ if ( fITSLoader == 0x0)
+ {
+ DoOpenError("Exiting due to problems with opening files.");
+ return 1;
+ }
+
+ Info("OpenNextSession","________________________________________________________");
+ Info("OpenNextSession","Found %d event(s) in directory %s",
+ fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
+ Float_t mf;
+ if (fUseMagFFromRun)
+ {
+ if (fRunLoader->LoadgAlice())
+ {
+ DoOpenError("Error occured while loading AliRun.");
+ return 1;
+ }
+ mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
+ Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
+ fRunLoader->UnloadgAlice();
+ }
+ else
+ {
+ Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
+ if (fMagneticField == 0x0)
+ {
+ Fatal("OpenNextSession","Magnetic field can not be 0.");
+ return 1;//pro forma
+ }
+ mf = fMagneticField*10.;
+ }
+ AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
+
+ if (fITSLoader->LoadTracks())
+ {
+ DoOpenError("Error occured while loading TPC tracks.");
+ return 1;
+ }
+
+ fCurrentEvent = 0;
+ return 0;
+}
/********************************************************************/
+void AliHBTReaderITSv2::DoOpenError( const char *va_(fmt), ...)
+{
+ // Does error display and clean-up in case error caught on Open Next Session
+ va_list ap;
+ va_start(ap,va_(fmt));
+ Error("OpenNextFile", va_(fmt), ap);
+ va_end(ap);
+
+ delete fRunLoader;
+ fRunLoader = 0x0;
+ fITSLoader = 0x0;
+ fCurrentDir++;
+}
+
+/********************************************************************/
#include "AliHBTReader.h"
-#include <TString.h>
-class TFile;
+class AliLoader;
+class AliRunLoader;
+class TString;
class AliHBTReaderITSv2: public AliHBTReader
{
virtual ~AliHBTReaderITSv2();
- Int_t Read(AliHBTRun*, AliHBTRun*);//reads tracks and particles and puts them in runs
+ void Rewind();
- AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
- AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles
- Int_t GetNumberOfPartEvents();//returns number of particle events
- Int_t GetNumberOfTrackEvents();//returns number of track events
-
- void SetMagneticField(Float_t mf){fMagneticField=mf;}
- void UseMagneticFieldFromRun(Bool_t flag = kTRUE){fUseMagFFromRun=flag;}
+ Bool_t ReadsTracks() const {return kTRUE;}
+ Bool_t ReadsParticles() const {return kTRUE;}
- protected:
+ void SetMagneticField(Float_t mf){fMagneticField=mf;}
+ void UseMagneticFieldFromRun(Bool_t flag = kTRUE){fUseMagFFromRun=flag;}
- AliHBTRun* fParticles; //!simulated particles
- AliHBTRun* fTracks; //!reconstructed tracks (particles)
+ protected:
- TString fFileName;//name of the file with galice.root
-
- Bool_t fIsRead;//!flag indicating if the data are already read
-
- Float_t fMagneticField;//magnetic field value that was enforced while reading
- Bool_t fUseMagFFromRun;//flag indicating if using field specified in gAlice (kTRUE)
+ Int_t ReadNext();//reads tracks and particles and puts them in runs
+ Int_t OpenNextFile();
+ void DoOpenError( const char *va_(fmt), ...);
+
+ TString fFileName;//name of the file with galice.root
+ AliRunLoader* fRunLoader;//!Run Loader
+ AliLoader* fITSLoader;//! ITS Loader
+
+ Float_t fMagneticField;//magnetic field value that was enforced while reading
+ Bool_t fUseMagFFromRun;//flag indicating if using field specified in gAlice (kTRUE)
// or enforece other defined by fMagneticField
- ClassDef(AliHBTReaderITSv2,1)
+
+ ClassDef(AliHBTReaderITSv2,2)
};
#endif
#include "AliHBTParticle.h"
#include "AliHBTParticleCut.h"
+AliHBTReaderInternal test;
ClassImp(AliHBTReaderInternal)
/********************************************************************/
-AliHBTReaderInternal::AliHBTReaderInternal()
+AliHBTReaderInternal::AliHBTReaderInternal():
+ fFileName(),
+ fPartBranch(0x0),
+ fTrackBranch(0x0),
+ fTree(0x0),
+ fFile(0x0)
{
//Defalut constructor
- fParticles = 0x0;
- fTracks = 0x0;
- fIsRead = kFALSE;
}
/********************************************************************/
-AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):fFileName(filename)
+AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):
+ fFileName(filename),
+ fPartBranch(0x0),
+ fTrackBranch(0x0),
+ fTree(0x0),
+ fFile(0x0)
{
//constructor
//filename - name of file to open
- fParticles = new AliHBTRun();
- fTracks = new AliHBTRun();
- fIsRead = kFALSE;
}
/********************************************************************/
AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
- AliHBTReader(dirs),fFileName(filename)
+ AliHBTReader(dirs),
+ fFileName(filename),
+ fPartBranch(0x0),
+ fTrackBranch(0x0),
+ fTree(0x0),
+ fFile(0x0)
{
//ctor
//dirs contains strings with directories to look data in
//filename - name of file to open
- fParticles = new AliHBTRun();
- fTracks = new AliHBTRun();
- fIsRead = kFALSE;
}
/********************************************************************/
AliHBTReaderInternal::~AliHBTReaderInternal()
{
//desctructor
- delete fParticles;
- delete fTracks;
+ delete fTree;
+ delete fFile;
}
/********************************************************************/
-
-AliHBTEvent* AliHBTReaderInternal::GetParticleEvent(Int_t n)
- {
- //returns nth event with simulated particles
- if (!fIsRead)
- if(Read(fParticles,fTracks))
- {
- Error("GetParticleEvent","Error in reading");
- return 0x0;
- }
- return fParticles->GetEvent(n);
- }
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderInternal::GetTrackEvent(Int_t n)
- {
- //returns nth event with reconstructed tracks
- if (!fIsRead)
- if(Read(fParticles,fTracks))
- {
- Error("GetTrackEvent","Error in reading");
- return 0x0;
- }
- return fTracks->GetEvent(n);
- }
-/********************************************************************/
-
-Int_t AliHBTReaderInternal::GetNumberOfPartEvents()
- {
- //returns number of events of particles
- if (!fIsRead)
- if ( Read(fParticles,fTracks))
- {
- Error("GetNumberOfPartEvents","Error in reading");
- return 0;
- }
- return fParticles->GetNumberOfEvents();
- }
-
-/********************************************************************/
-Int_t AliHBTReaderInternal::GetNumberOfTrackEvents()
- {
- //returns number of events of tracks
- if (!fIsRead)
- if(Read(fParticles,fTracks))
- {
- Error("GetNumberOfTrackEvents","Error in reading");
- return 0;
- }
- return fTracks->GetNumberOfEvents();
- }
+
+void AliHBTReaderInternal::Rewind()
+{
+ delete fTree;
+ fTree = 0x0;
+ delete fFile;
+ fFile = 0x0;
+ fCurrentDir = 0;
+ fNEventsRead= 0;
+}
/********************************************************************/
-Int_t AliHBTReaderInternal::Read(AliHBTRun* particles, AliHBTRun *tracks)
+Int_t AliHBTReaderInternal::ReadNext()
{
//reads data and puts put to the particles and tracks objects
//reurns 0 if everything is OK
//
- Info("Read","");
+ Info("ReadNext","");
Int_t i; //iterator and some temprary values
- Int_t totalnevents = 0; //total number of read events
- Int_t nevents = 0;
- Int_t currentdir = 0;
- Int_t Ndirs;
Int_t counter;
- TFile *aFile;//file with tracks
AliHBTParticle* tpart = 0x0, *ttrack = 0x0;
+
TDatabasePDG* pdgdb = TDatabasePDG::Instance();
if (pdgdb == 0x0)
{
- Error("Read","Can not get PDG Particles Data Base");
- return 1;
- }
- if (!particles) //check if an object is instatiated
- {
- Error("Read"," particles object must instatiated before passing it to the reader");
- return 1;
- }
- if (!tracks) //check if an object is instatiated
- {
- Error("Read"," tracks object must instatiated before passing it to the reader");
+ Error("ReadNext","Can not get PDG Particles Data Base");
return 1;
}
- particles->Reset();//clear runs == delete all old events
- tracks->Reset();
-
- if (fDirs) //if array with directories is supplied by user
- {
- Ndirs = fDirs->GetEntries(); //get the number if directories
- }
- else
- {
- Ndirs = 0; //if the array is not supplied read only from current directory
- }
-
- TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
- TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
-// pbuffer->BypassStreamer(kFALSE);
-// tbuffer->BypassStreamer(kFALSE);
+
+ if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
+ if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
+
+ fParticlesEvent->Reset();
+ fTracksEvent->Reset();
do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
{
-
- if( (i=OpenFile(aFile, currentdir)) )
+ if (fTree == 0x0)
+ if( (i=OpenNextFile()) )
+ {
+ Error("ReadNext","Skippimg directory due to problems with opening files. Errorcode %d",i);
+ fCurrentDir++;
+ continue;
+ }
+ if (fCurrentEvent == (Int_t)fTree->GetEntries())
{
- Error("Read","Skippimg directory due to problems with opening files. Errorcode %d",i);
- currentdir++;
+ delete fTree;
+ fTree = 0x0;
+ delete fFile;
+ fFile = 0x0;
+ fPartBranch = 0x0;
+ fTrackBranch= 0x0;
+ fCurrentDir++;
continue;
}
/***************************/
/***************************/
/***************************/
- TTree* tree = (TTree*)aFile->Get("data");
- if (tree == 0x0)
- {
- Error("Read","Can not get the tree");
- currentdir++;
- continue;
- }
-
- TBranch *trackbranch=tree->GetBranch("tracks");//get the branch with tracks
- if (trackbranch == 0x0) ////check if we got the branch
- {//if not return with error
- Warning("Read","Can't find a branch with tracks !\n");
- }
- else
- {
- trackbranch->SetAddress(&tbuffer);
- }
+
+ Info("ReadNext","Event %d",fCurrentEvent);
+ fTree->GetEvent(fCurrentEvent);
- TBranch *partbranch=tree->GetBranch("particles");//get the branch with particles
- if (partbranch == 0x0) ////check if we got the branch
- {//if not return with error
- Warning("Read","Can't find a branch with particles !\n");
- }
- else
- {
- partbranch->SetAddress(&pbuffer);
- }
-
- nevents = (Int_t)tree->GetEntries();
- Info("Read","________________________________________________________");
- Info("Read","Found %d event(s) in directory %s",nevents,GetDirName(currentdir).Data());
-
- for (Int_t currentEvent =0; currentEvent<nevents;currentEvent++)
- {
- Info("Read","Event %d",currentEvent);
- tree->GetEvent(currentEvent);
-
- counter = 0;
- if (partbranch && trackbranch)
- {
- for(i = 0; i < pbuffer->GetEntries(); i++)
+ counter = 0;
+ if (fPartBranch && fTrackBranch)
+ {
+ for(i = 0; i < fPartBuffer->GetEntries(); i++)
+ {
+ tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
+ ttrack = dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
+
+ if( tpart == 0x0 ) continue; //if returned pointer is NULL
+
+ if (ttrack->GetUID() != tpart->GetUID())
{
- tpart = dynamic_cast<AliHBTParticle*>(pbuffer->At(i));
- ttrack = dynamic_cast<AliHBTParticle*>(tbuffer->At(i));
-
- if( tpart == 0x0 ) continue; //if returned pointer is NULL
-
- if (ttrack->GetUID() != tpart->GetUID())
- {
- Error("Read","Sth. is wrong: Track and Particle has different UID.");
- Error("Read","They probobly do not correspond to each other.");
- }
-
- for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+ Error("ReadNext","Sth. is wrong: Track and Particle has different UID.");
+ Error("ReadNext","They probobly do not correspond to each other.");
+ }
+
+ for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+ {
+ if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
+ if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
+ //if not take next partilce
+ AliHBTParticle* part = new AliHBTParticle(*tpart);
+ part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
+ if( Pass(part) )
{
- if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
- if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
- //if not take next partilce
- AliHBTParticle* part = new AliHBTParticle(*tpart);
- part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
- if( Pass(part) )
- {
- delete part;
- continue;
- }
- AliHBTParticle* track = new AliHBTParticle(*ttrack);
-
- particles->AddParticle(totalnevents,part);//put track and particle on the run
- tracks->AddParticle(totalnevents,track);
- counter++;
+ delete part;
+ continue;
}
- }
- Info("Read"," Read: %d particles and tracks.",counter);
- }
- else
- {
- if (partbranch)
- {
- Info("Read","Found %d particles in total.",pbuffer->GetEntries());
- for(i = 0; i < pbuffer->GetEntries(); i++)
- {
- tpart = dynamic_cast<AliHBTParticle*>(pbuffer->At(i));
- if(tpart == 0x0) continue; //if returned pointer is NULL
-
- for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+ AliHBTParticle* track = new AliHBTParticle(*ttrack);
+
+ fParticlesEvent->AddParticle(part);
+ fTracksEvent->AddParticle(track);
+
+ counter++;
+ }
+ }
+ Info("ReadNext"," Read: %d particles and tracks.",counter);
+ }
+ else
+ {
+ if (fPartBranch)
+ {
+ Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
+ for(i = 0; i < fPartBuffer->GetEntries(); i++)
+ {
+ tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
+ if(tpart == 0x0) continue; //if returned pointer is NULL
+
+ for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+ {
+ if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
+ if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
+
+ AliHBTParticle* part = new AliHBTParticle(*tpart);
+ part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
+ if( Pass(part) )
{
- if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
- if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
-
- AliHBTParticle* part = new AliHBTParticle(*tpart);
- part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
- if( Pass(part) )
- {
- delete part;
- continue;
- }
- particles->AddParticle(totalnevents,part);//put track and particle on the run
- counter++;
+ delete part;
+ continue;
}
- }
- Info("Read"," Read: %d particles.",counter);
- }
- else Info("Read"," Read: 0 particles.");
-
- if (trackbranch)
- {
- for(i = 0; i < tbuffer->GetEntries(); i++)
- {
- tpart = dynamic_cast<AliHBTParticle*>(tbuffer->At(i));
- if(tpart == 0x0) continue; //if returned pointer is NULL
- for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+ fParticlesEvent->AddParticle(part);
+ counter++;
+ }
+ }
+ Info("ReadNext"," Read: %d particles.",counter);
+ }
+ else Info("ReadNext"," Read: 0 particles.");
+
+ if (fTrackBranch)
+ {
+ for(i = 0; i < fTrackBuffer->GetEntries(); i++)
+ {
+ tpart = dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
+ if(tpart == 0x0) continue; //if returned pointer is NULL
+ for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
+ {
+ if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
+ if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
+
+ AliHBTParticle* part = new AliHBTParticle(*tpart);
+ part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
+ if( Pass(part) )
{
- if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
- if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
-
- AliHBTParticle* part = new AliHBTParticle(*tpart);
- part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
- if( Pass(part) )
- {
- delete part;
- continue;
- }
- tracks->AddParticle(totalnevents,part);//put track and particle on the run
- counter++;
+ delete part;
+ continue;
}
- }
- Info("Read"," Read: %d tracks",counter);
- }
- else Info("Read"," Read: 0 tracks.");
- }
- totalnevents++;
- }
-
- /***************************/
- /***************************/
- /***************************/
- currentdir++;
- delete tree;
- aFile->Close();
- delete aFile;
- aFile = 0x0;
-
- }while(currentdir < Ndirs);
+ fTracksEvent->AddParticle(part);
+ counter++;
+ }
+ }
+ Info("ReadNext"," Read: %d tracks",counter);
+ }
+ else Info("ReadNext"," Read: 0 tracks.");
+ }
+ fPartBuffer->Delete();
+ fTrackBuffer->Delete();
+
+ fCurrentEvent++;
+ fNEventsRead++;
+
+ return 0;
- delete pbuffer;
- delete tbuffer;
- fIsRead = kTRUE;
- return 0;
+ }while(fCurrentEvent < GetNumberOfDirs());
+
+ return 1;//no more directories to read
}
/********************************************************************/
-Int_t AliHBTReaderInternal::OpenFile(TFile*& aFile,Int_t event)
+Int_t AliHBTReaderInternal::OpenNextFile()
{
-
- const TString& dirname = GetDirName(event);
+ //open file in current directory
+ const TString& dirname = GetDirName(fCurrentEvent);
if (dirname == "")
{
- Error("OpenFile","Can not get directory name");
+ Error("OpenNextFile","Can not get directory name");
return 4;
}
TString filename = dirname +"/"+ fFileName;
- aFile = TFile::Open(filename.Data());
- if ( aFile == 0x0 )
+ fFile = TFile::Open(filename.Data());
+ if ( fFile == 0x0 )
{
- Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
+ Error("OpenNextFile","Can't open file named %s",filename.Data());
return 1;
}
- if (!aFile->IsOpen())
+ if (fFile->IsOpen() == kFALSE)
{
- Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
+ Error("OpenNextFile","Can't open filenamed %s",filename.Data());
return 1;
}
+
+ fTree = (TTree*)fFile->Get("data");
+ if (fTree == 0x0)
+ {
+ Error("OpenNextFile","Can not get the tree.");
+ return 1;
+ }
+
+
+ fTrackBranch = fTree->GetBranch("tracks");//get the branch with tracks
+ if (fTrackBranch == 0x0) ////check if we got the branch
+ {//if not return with error
+ Info("OpenNextFile","Can't find a branch with tracks !\n");
+ }
+ else
+ {
+ if (fTrackBuffer == 0x0)
+ fTrackBuffer = new TClonesArray("AliHBTParticle",15000);
+ fTrackBranch->SetAddress(&fTrackBuffer);
+ }
+
+ fPartBranch = fTree->GetBranch("particles");//get the branch with particles
+ if (fPartBranch == 0x0) ////check if we got the branch
+ {//if not return with error
+ Info("OpenNextFile","Can't find a branch with particles !\n");
+ }
+ else
+ {
+ if (fPartBuffer == 0x0)
+ fPartBuffer = new TClonesArray("AliHBTParticle",15000);
+ fPartBranch->SetAddress(&fPartBuffer);
+ }
+
+ Info("OpenNextFile","________________________________________________________");
+ Info("OpenNextFile","Found %d event(s) in directory %s",
+ (Int_t)fTree->GetEntries(),GetDirName(fCurrentEvent).Data());
+
+ fCurrentEvent = 0;
+
return 0;
}
/********************************************************************/
TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
-// pbuffer->BypassStreamer(kFALSE);
-// tbuffer->BypassStreamer(kFALSE);
TClonesArray &particles = *pbuffer;
TClonesArray &tracks = *tbuffer;
#include <TString.h>
class TFile;
+class TTree;
+class TBranch;
class TClonesArray;
class AliHBTReaderInternal: public AliHBTReader
AliHBTReaderInternal(TObjArray* dirs, const char *filename);
virtual ~AliHBTReaderInternal();
- Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
- static Int_t Write(AliHBTReader* reader,const char* outfile = "data.root", Bool_t multcheck = kFALSE);//reads tracks from runs and writes them to file
+ void Rewind();
+
+ Bool_t ReadsTracks() const {return kTRUE;}
+ Bool_t ReadsParticles() const {return kTRUE;}
- AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
- AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles
- Int_t GetNumberOfPartEvents();//returns number of particle events
- Int_t GetNumberOfTrackEvents();//returns number of track events
+ static Int_t Write(AliHBTReader* reader,const char* outfile = "data.root", Bool_t multcheck = kFALSE);//reads tracks from runs and writes them to file
protected:
- AliHBTRun* fParticles; //!simulated particles
- AliHBTRun* fTracks; //!reconstructed tracks (particles)
- Bool_t fIsRead;//!flag indicating if the data are already read
- TString fFileName;//name of the file with tracks
-
- Int_t OpenFile(TFile*& aFile,Int_t event);//opens file to be read for given event
+
+ TString fFileName;//name of the file with tracks
+ TBranch* fPartBranch;//!branch with particles
+ TBranch* fTrackBranch;//!branch with tracks
+ TTree* fTree;//!tree
+ TFile* fFile;//!file
+ TClonesArray* fPartBuffer;//!buffer array that tree is read to
+ TClonesArray* fTrackBuffer;//!
+
+ Int_t ReadNext();//reads next event
+ Int_t OpenNextFile();//opens file to be read for given event
static Bool_t FindIndex(TClonesArray* arr,Int_t idx);
- ClassDef(AliHBTReaderInternal,1)
+ ClassDef(AliHBTReaderInternal,2)
};
#endif
#include "AliHBTReaderKineTree.h"
-#include <Riostream.h>
#include <TString.h>
#include <TObjString.h>
#include <TTree.h>
/**********************************************************/
const TString AliHBTReaderKineTree::fgkEventFolderName("HBTReaderKineTree");
-AliHBTReaderKineTree::AliHBTReaderKineTree():fFileName("galice.root")
+AliHBTReaderKineTree::AliHBTReaderKineTree():
+ fFileName("galice.root"),
+ fRunLoader(0x0)
{
- fParticles = new AliHBTRun();
- fIsRead = kFALSE;
+ //ctor
}
+/**********************************************************/
-AliHBTReaderKineTree::
-AliHBTReaderKineTree(TString& fname):fFileName(fname)
+AliHBTReaderKineTree::AliHBTReaderKineTree(TString& fname):
+ fFileName(fname),
+ fRunLoader(0x0)
{
- fParticles = new AliHBTRun();
- fIsRead = kFALSE;
+ //ctor
}
-
-
/**********************************************************/
-AliHBTReaderKineTree::
-AliHBTReaderKineTree(TObjArray* dirs,const Char_t *filename):AliHBTReader(dirs),fFileName(filename)
+
+AliHBTReaderKineTree::AliHBTReaderKineTree(TObjArray* dirs,const Char_t *filename):
+ AliHBTReader(dirs),
+ fFileName(filename),
+ fRunLoader(0x0)
{
- fParticles = new AliHBTRun();
- fIsRead = kFALSE;
+ //ctor
}
/**********************************************************/
-AliHBTEvent* AliHBTReaderKineTree::GetParticleEvent(Int_t n)
- {
- //returns Nth event with simulated particles
- if (!fIsRead)
- if(Read(fParticles,0x0))
- {
- Error("GetParticleEvent","Error in reading");
- return 0x0;
- }
-
- return fParticles->GetEvent(n);
- }
-
-/********************************************************************/
-
-Int_t AliHBTReaderKineTree::GetNumberOfPartEvents()
- {
- //returns number of events of particles
- if (!fIsRead)
- if(Read(fParticles,0x0))
- {
- Error("GetNumberOfPartEvents","Error in reading");
- return 0;
- }
- return fParticles->GetNumberOfEvents();
- }
-
+AliHBTReaderKineTree::~AliHBTReaderKineTree()
+{
+ //dtor
+ delete fRunLoader;
+}
+/**********************************************************/
+void AliHBTReaderKineTree::Rewind()
+{
+ delete fRunLoader;
+ fRunLoader = 0x0;
+ fCurrentDir = 0;
+ fNEventsRead= 0;
+}
/**********************************************************/
-Int_t AliHBTReaderKineTree::Read(AliHBTRun* particles, AliHBTRun* /*tracks*/)
+
+Int_t AliHBTReaderKineTree::ReadNext()
{
//Reads Kinematics Tree
Info("Read","");
-
- if (!particles) //check if an object is instatiated
- {
- Error("Read"," particles object must instatiated before passing it to the reader");
- return 1;
- }
- particles->Reset();//clear runs == delete all old events
-
- Int_t Ndirs;//total number of directories specified in fDirs array
- Int_t Nevents; //number of events read in current directory
- Int_t totalNevents = 0; //total number of read events
- Int_t currentdir = 0; //number of current directory name is fDirs array
-
- if (fDirs) //if array with directories is supplied by user
- {
- Ndirs = fDirs->GetEntries(); //get the number if directories
- }
- else
- {
- Ndirs = 0; //if the array is not supplied read only from current directory
- }
+ if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
+ fParticlesEvent->Reset();
do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
{
- Info("Read","________________________________________________________");
- AliRunLoader * rl = OpenFile(currentdir);
-
- if (rl == 0x0)
+ if (fRunLoader == 0x0)
+ if (OpenNextFile())
+ {
+ fCurrentDir++;
+ continue;
+ }
+
+ if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
{
- Error("Read","Cannot open session");
- return 2;
+ //read next directory
+ delete fRunLoader;//close current session
+ fRunLoader = 0x0;//assure pointer is null
+ fCurrentDir++;//go to next dir
+ continue;//directory counter is increased inside in case of error
}
- rl->LoadHeader();
- rl->LoadKinematics("READ");
- Nevents = rl->GetNumberOfEvents();
+
+ Info("ReadNext","Reading Event %d",fCurrentEvent);
+
+ fRunLoader->GetEvent(fCurrentEvent);
- for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
+ AliStack* stack = fRunLoader->Stack();
+ if (!stack)
{
-// cout<<"processing event "<<currentEvent<<" in current dir."<<endl;
- rl->GetEvent(currentEvent);
- AliStack* stack = rl->Stack();
- if (!stack)
- {
- Error("Read","Can not get stack for event %d",currentEvent);
- continue;
- }
- Int_t npart = stack->GetNtrack();
- Int_t nnn = 0;
- for (Int_t i = 0;i<npart; i++)
- {
-
+ Error("ReadNext","Can not get stack for event %d",fCurrentEvent);
+ continue;
+ }
+ Int_t npart = stack->GetNtrack();
+ for (Int_t i = 0;i<npart; i++)
+ {
TParticle * p = stack->Particle(i);
// if (p->GetFirstMother() >= 0) continue; do not apply with pythia etc
AliHBTParticle* part = new AliHBTParticle(*p,i);
if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
//if it does not delete it and take next good track
- particles->AddParticle(totalNevents,part);//put track and particle on the run
-
- if ( (nnn%100) == 0)
- {
- cout<<nnn<<"\r";
- fflush(0);
- }
- nnn++;
- }
- Info("Read","Total read %d",nnn);
- totalNevents++;
- }
- delete rl;
- currentdir++;
- }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
- fIsRead = kTRUE;
- return 0;
+ fParticlesEvent->AddParticle(part);//put particle in event
+ }
+ Info("ReadNext","Read %d particles from event %d (event %d in dir %d).",
+ fParticlesEvent->GetNumberOfParticles(),
+ fNEventsRead,fCurrentEvent,fCurrentDir);
+
+ fCurrentEvent++;
+ fNEventsRead++;
+ return 0;
+ }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
+
+ return 1;
}
+/**********************************************************/
-AliRunLoader* AliHBTReaderKineTree::OpenFile(Int_t n)
+Int_t AliHBTReaderKineTree::OpenNextFile()
{
//opens file with kine tree
-
- const TString& dirname = GetDirName(n);
+ Info("OpenNextFile","________________________________________________________");
+
+ const TString& dirname = GetDirName(fCurrentEvent);
if (dirname == "")
{
- Error("OpenFiles","Can not get directory name");
- return 0x0;
+ Error("OpenNextFile","Can not get directory name");
+ return 1;
}
TString filename = dirname +"/"+ fFileName;
- AliRunLoader *ret = AliRunLoader::Open(filename.Data(),fgkEventFolderName,"READ");
+ fRunLoader = AliRunLoader::Open(filename.Data(),fgkEventFolderName,"READ");
- if ( ret == 0x0)
+ if ( fRunLoader == 0x0)
{
- Error("OpenFiles","Can't open session from file %s",filename.Data());
+ Error("OpenNextFile","Can't open session from file %s",filename.Data());
return 0x0;
}
-
- return ret;
+
+ if (fRunLoader->GetNumberOfEvents() <= 0)
+ {
+ Error("OpenNextFile","There is no events in this directory.");
+ delete fRunLoader;
+ fRunLoader = 0x0;
+ return 1;
+ }
+
+ if (fRunLoader->LoadKinematics())
+ {
+ Error("OpenNextFile","Error occured while loading kinematics.");
+ return 1;
+ }
+
+ fCurrentEvent = 0;
+ return 0;
}
AliHBTReaderKineTree(TString&);
AliHBTReaderKineTree(TObjArray*,const Char_t *filename="galice.root");
- virtual ~AliHBTReaderKineTree(){}
+ virtual ~AliHBTReaderKineTree();
- Int_t Read(AliHBTRun* particles, AliHBTRun* /*tracks*/);//reads tracks and particles and puts them in runs
-
- AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
- AliHBTEvent* GetTrackEvent(Int_t){return 0x0;}//returns pointer to event with particles
- Int_t GetNumberOfPartEvents();//returns number of particle events
- Int_t GetNumberOfTrackEvents(){return 0;}//returns number of track events
-
+ void Rewind();
+
+ Bool_t ReadsTracks() const {return kFALSE;}
+ Bool_t ReadsParticles() const {return kTRUE;}
protected:
- TString fFileName;
- AliHBTRun* fParticles; //!simulated particles
-
- AliRunLoader* OpenFile(Int_t);
-
- Bool_t fIsRead;//!flag indicating if the data are already read
+ Int_t ReadNext();//reads tracks and particles and puts them in runs
+ Int_t OpenNextFile();
+
+ TString fFileName;//file name
+ AliRunLoader* fRunLoader;
static const TString fgkEventFolderName;
private:
public:
- ClassDef(AliHBTReaderKineTree,1)
+ ClassDef(AliHBTReaderKineTree,2)
};
#endif
#include <TFile.h>
#include <TParticle.h>
+#include <Varargs.h>
+
#include <AliRun.h>
#include <AliLoader.h>
#include <AliStack.h>
#include <AliTPCtrack.h>
#include <AliTPCParam.h>
#include <AliTPCtracker.h>
+#include <AliTPCLoader.h>
#include "AliHBTRun.h"
#include "AliHBTEvent.h"
#include "AliHBTParticle.h"
#include "AliHBTParticleCut.h"
-
ClassImp(AliHBTReaderTPC)
//______________________________________________
//
//
//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
//Piotr.Skowronski@cern.ch
-AliHBTReaderTPC::AliHBTReaderTPC():fFileName("galice.root")
+AliHBTReaderTPC::AliHBTReaderTPC():
+ fFileName("galice.root"),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
{
//constructor,
//Defaults:
// galicefilename = "" - this means: Do not open gAlice file -
// just leave the global pointer untouched
- fParticles = 0x0;
- fTracks = 0x0;
- fIsRead = kFALSE;
}
-AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):fFileName(galicefilename)
+AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
{
//constructor,
//Defaults:
// galicefilename = "" - this means: Do not open gAlice file -
// just leave the global pointer untouched
- fParticles = new AliHBTRun();
- fTracks = new AliHBTRun();
- fIsRead = kFALSE;
}
/********************************************************************/
AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
- AliHBTReader(dirs), fFileName(galicefilename)
-
+ AliHBTReader(dirs),
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE)
{
//constructor,
//Defaults:
// galicefilename = "" - this means: Do not open gAlice file -
// just leave the global pointer untached
- fParticles = new AliHBTRun();
- fTracks = new AliHBTRun();
- fIsRead = kFALSE;
}
/********************************************************************/
AliHBTReaderTPC::~AliHBTReaderTPC()
- {
+{
//desctructor
- delete fParticles;
- delete fTracks;
- }
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
- {
- //returns Nth event with simulated particles
- if (!fIsRead)
- {
- if (fParticles == 0x0) fParticles = new AliHBTRun();
- if (fTracks == 0x0) fTracks = new AliHBTRun();
-
- if(Read(fParticles,fTracks))
- {
- Error("GetParticleEvent","Error in reading");
- return 0x0;
- }
- }
- return (fParticles)?fParticles->GetEvent(n):0x0;
- }
-/********************************************************************/
-
-AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
- {
- //returns Nth event with reconstructed tracks
- if (!fIsRead)
- {
- if (fParticles == 0x0) fParticles = new AliHBTRun();
- if (fTracks == 0x0) fTracks = new AliHBTRun();
- if(Read(fParticles,fTracks))
- {
- Error("GetTrackEvent","Error in reading");
- return 0x0;
- }
- }
- return (fTracks)?fTracks->GetEvent(n):0x0;
- }
-/********************************************************************/
-
-Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
- {
- //returns number of events of particles
- if (!fIsRead)
- {
- if (fParticles == 0x0) fParticles = new AliHBTRun();
- if (fTracks == 0x0) fTracks = new AliHBTRun();
- if ( Read(fParticles,fTracks))
- {
- Error("GetNumberOfPartEvents","Error in reading");
- return 0;
- }
- }
- return (fParticles)?fParticles->GetNumberOfEvents():0;
- }
-
+ delete fRunLoader;
+}
/********************************************************************/
-Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
- {
- //returns number of events of tracks
- if (!fIsRead)
- {
- if (fParticles == 0x0) fParticles = new AliHBTRun();
- if (fTracks == 0x0) fTracks = new AliHBTRun();
- if(Read(fParticles,fTracks))
- {
- Error("GetNumberOfTrackEvents","Error in reading");
- return 0;
- }
- }
- return (fTracks)?fTracks->GetNumberOfEvents():0;
- }
+
+void AliHBTReaderTPC::Rewind()
+{
+ delete fRunLoader;
+ fRunLoader = 0x0;
+ fCurrentDir = 0;
+ fNEventsRead= 0;
+}
/********************************************************************/
-
-Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
+Int_t AliHBTReaderTPC::ReadNext()
{
//reads data and puts put to the particles and tracks objects
//reurns 0 if everything is OK
//
Info("Read","");
- Int_t Nevents = 0;
- Int_t totalNevents = 0;
- if (!particles) //check if an object is instatiated
- {
- Error("Read"," particles object must instatiated before passing it to the reader");
- }
- if (!tracks) //check if an object is instatiated
- {
- Error("Read"," tracks object must instatiated before passing it to the reader");
- }
- particles->Reset();//clear runs == delete all old events
- tracks->Reset();
TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
tarray->SetOwner(); //set the ownership of the objects it contains
//when array is is deleted or cleared all objects
//that it contains are deleted
- Int_t currentdir = 0;
- Int_t Ndirs;
- if (fDirs) //if array with directories is supplied by user
- {
- Ndirs = fDirs->GetEntries(); //get the number if directories
- }
- else
- {
- Ndirs = 0; //if the array is not supplied read only from current directory
- }
+ if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
+ if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
+
+ fParticlesEvent->Reset();
+ fTracksEvent->Reset();
do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
{
- TString filename = GetDirName(currentdir);
- if (filename.IsNull())
- {
- Error("Read","Can not get directory name");
- return 4;
- }
- filename = filename +"/"+ fFileName;
- AliRunLoader* rl = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
- if( rl == 0x0)
+
+ if (fRunLoader == 0x0)
+ if (OpenNextSession()) continue;//directory counter is increased inside in case of error
+
+ if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
{
- Error("Read","Can not open session.");
- currentdir++;
- continue;
+ //read next directory
+ delete fRunLoader;//close current session
+ fRunLoader = 0x0;//assure pointer is null
+ fCurrentDir++;//go to next dir
+ continue;//directory counter is increased inside in case of error
}
+
+ Info("ReadNext","Reading Event %d",fCurrentEvent);
- rl->LoadHeader();
- rl->LoadKinematics();
- AliLoader* tpcl = rl->GetLoader("TPCLoader");
-
- if ( tpcl== 0x0)
+ fRunLoader->GetEvent(fCurrentEvent);
+ TTree *tracktree = fTPCLoader->TreeT();//get the tree
+ if (!tracktree) //check if we got the tree
+ {//if not return with error
+ Error("ReadNext","Can't get a tree with TPC tracks !\n");
+ continue;
+ }
+
+ TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
+ if (!trackbranch) ////check if we got the branch
+ {//if not return with error
+ Error("ReadNext","Can't get a branch with TPC tracks !\n");
+ continue;
+ }
+ Int_t ntpctracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
+ Info("ReadNext","Found %d TPC tracks.",ntpctracks);
+ //Copy tracks to array
+
+ AliTPCtrack *iotrack=0;
+
+ for (Int_t i=0; i<ntpctracks; i++) //loop over all tpc tracks
{
- Error("Read","Exiting due to problems with opening files.");
- currentdir++;
- continue;
+ iotrack=new AliTPCtrack; //create new tracks
+ trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
+ tracktree->GetEvent(i); //stream track i to the iotrack
+ tarray->AddLast(iotrack); //put the track in the array
}
- Nevents = rl->GetNumberOfEvents();
-
- if (Nevents > 0)//check if tree E exists
+
+ Double_t xk;
+ Double_t par[5];
+ Float_t phi, lam, pt;//angles and transverse momentum
+ Int_t label; //label of the current track
+
+ AliStack* stack = fRunLoader->Stack();
+ if (stack == 0x0)
{
- Info("Read","________________________________________________________");
- Info("Read","Found %d event(s) in directory %s",Nevents,GetDirName(currentdir).Data());
- rl->LoadgAlice();
- Info("Read","Setting Magnetic Field: B=%fT",rl->GetAliRun()->Field()->SolenoidField());
- AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
- rl->UnloadgAlice();
- }
- else
- {//if not return an error
- Error("Read","Can not find Header tree (TreeE) in gAlice");
- currentdir++;
+ Error("ReadNext","Can not get stack for current event",fCurrentEvent);
+ fCurrentEvent++;
continue;
}
+ stack->Particles();
+
+ for (Int_t i=0; i<ntpctracks; i++) //loop over all good tracks
+ {
+ iotrack = (AliTPCtrack*)tarray->At(i);
+ label = iotrack->GetLabel();
+
+ if (label < 0) continue;
+
+ TParticle *p = (TParticle*)stack->Particle(label);
+
+ if(p == 0x0) continue; //if returned pointer is NULL
+ if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
+
+ if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
+ //if not take next partilce
+
+ AliHBTParticle* part = new AliHBTParticle(*p,i);
+ if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
+ //if it does not delete it and take next good track
+
+// iotrack->PropagateTo(3.,0.0028,65.19);
+// iotrack->PropagateToVertex(36.66,1.2e-3);//orig
+ iotrack->PropagateToVertex(50.,0.0353);
+
+ iotrack->GetExternalParameters(xk,par); //get properties of the track
+ phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
+ if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+ if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+ lam=par[3];
+ pt=1.0/TMath::Abs(par[4]);
+
+ Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
+ Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
+ Double_t tpz = pt * lam; //track z coordinate of momentum
+
+ Double_t mass = p->GetMass();
+ Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
+
+ AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+ if(Pass(track))//check if meets all criteria of any of our cuts
+ //if it does not delete it and take next good track
+ {
+ delete track;
+ delete part;
+ continue;
+ }
+ fParticlesEvent->AddParticle(part);
+ fTracksEvent->AddParticle(track);
+ }
+
+ Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
+ fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
+ fNEventsRead,fCurrentEvent,fCurrentDir);
- rl->CdGAFile();
- AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
-
- if (!TPCParam)
- {
- TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
- if (!TPCParam)
- {
- Error("Read","TPC parameters have not been found !\n");
- delete rl;
- rl = 0x0;
- currentdir++;
- continue;
+ fCurrentEvent++;
+ fNEventsRead++;
+ delete tarray;
+ return 0;
+ }while(fCurrentDir < GetNumberOfDirs());
+
+ delete tarray;
+ return 1;
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderTPC::OpenNextSession()
+{
+ TString filename = GetDirName(fCurrentDir);
+ if (filename.IsNull())
+ {
+ DoOpenError("Can not get directory name");
+ return 1;
+ }
+ filename = filename +"/"+ fFileName;
+ fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
+ if( fRunLoader == 0x0)
+ {
+ DoOpenError("Can not open session.");
+ return 1;
+ }
+
+ fRunLoader->LoadHeader();
+ if (fRunLoader->GetNumberOfEvents() <= 0)
+ {
+ DoOpenError("There is no events in this directory.");
+ return 1;
+ }
+
+ if (fRunLoader->LoadKinematics())
+ {
+ DoOpenError("Error occured while loading kinematics.");
+ return 1;
+ }
+ fTPCLoader = (AliTPCLoader*)fRunLoader->GetLoader("TPCLoader");
+ if ( fTPCLoader == 0x0)
+ {
+ DoOpenError("Exiting due to problems with opening files.");
+ return 1;
+ }
+ Info("OpenNextSession","________________________________________________________");
+ Info("OpenNextSession","Found %d event(s) in directory %s",
+ fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
+ Float_t mf;
+ if (fUseMagFFromRun)
+ {
+ fRunLoader->LoadgAlice();
+ mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
+ Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
+ fRunLoader->UnloadgAlice();
+ }
+ else
+ {
+ Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
+ if (fMagneticField == 0x0)
+ {
+ Fatal("OpenNextSession","Magnetic field can not be 0.");
+ return 1;//pro forma
}
- }
+ mf = fMagneticField*10.;
+ }
+ AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
- tpcl->LoadTracks();
-
- for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
- {
- Info("Read","Reading Event %d",currentEvent);
- /**************************************/
- /**************************************/
- /**************************************/
- rl->GetEvent(currentEvent);
- TTree *tracktree = tpcl->TreeT();//get the tree
- if (!tracktree) //check if we got the tree
- {//if not return with error
- Error("Read","Can't get a tree with TPC tracks !\n");
- continue;
- }
-
- TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
- if (!trackbranch) ////check if we got the branch
- {//if not return with error
- Error("Read","Can't get a branch with TPC tracks !\n");
- continue;
- }
- Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
- Info("Read","Found %d TPC tracks.",NTPCtracks);
- //Copy tracks to array
-
- AliTPCtrack *iotrack=0;
-
-printf("This method is not converted to the NewIO !\n"); //I.B.
- //AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent,AliConfig::fgkDefaultEventFolderName);//create the tacker for this event
- AliTPCtracker *tracker = new AliTPCtracker(TPCParam); //I.B.
- if (!tracker) //check if it has created succeffuly
- {//if not return with error
- Error("Read","Can't get a tracker !\n");
- continue;
- }
- tracker->LoadClusters(0);//I.Belikov, "0" must be a pointer to a tree
-
- for (Int_t i=0; i<NTPCtracks; i++) //loop over all tpc tracks
- {
- iotrack=new AliTPCtrack; //create new tracks
- trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
- tracktree->GetEvent(i); //stream track i to the iotrack
- tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track
- //which is the label of corresponding simulated particle
- tarray->AddLast(iotrack); //put the track in the array
- }
-
- delete tracker; //delete tracker
-
- tracker = 0x0;
- trackbranch = 0x0;
- tracktree = 0x0;
-
- Double_t xk;
- Double_t par[5];
- Float_t phi, lam, pt;//angles and transverse momentum
- Int_t label; //label of the current track
-
- rl->Stack()->Particles();
-
- for (Int_t i=0; i<NTPCtracks; i++) //loop over all good tracks
- {
- iotrack = (AliTPCtrack*)tarray->At(i);
- label = iotrack->GetLabel();
-
- if (label < 0) continue;
-
- TParticle *p = (TParticle*)rl->Stack()->Particle(label);
-
- if(p == 0x0) continue; //if returned pointer is NULL
- if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
-
- if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
- //if not take next partilce
-
- AliHBTParticle* part = new AliHBTParticle(*p,i);
- if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
- //if it does not delete it and take next good track
-
- iotrack->PropagateToVertex();
-
- iotrack->GetExternalParameters(xk,par); //get properties of the track
- phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
- if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
- if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
- lam=par[3];
- pt=1.0/TMath::Abs(par[4]);
-
- Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
- Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
- Double_t tpz = pt * lam; //track z coordinate of momentum
-
- Double_t mass = p->GetMass();
- Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
-
- AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
- if(Pass(track))//check if meets all criteria of any of our cuts
- //if it does not delete it and take next good track
- {
- delete track;
- delete part;
- continue;
- }
- particles->AddParticle(totalNevents,part);//put track and particle on the run
- tracks->AddParticle(totalNevents,track);
-
- }
- tarray->Clear(); //clear the array
-
- /**************************************/
- /**************************************/
- /**************************************/
- totalNevents++;
- }
-
- //save environment (resouces) --
- //clean your place after the work
- delete rl;
- currentdir++;
- }while(currentdir < Ndirs);
- delete tarray;
- fIsRead = kTRUE;
+ fRunLoader->CdGAFile();
+ AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
+ if (!TPCParam)
+ {
+ TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
+ if (!TPCParam)
+ {
+ DoOpenError("TPC parameters have not been found !\n");
+ return 1;
+ }
+ }
+
+ if (fTPCLoader->LoadTracks())
+ {
+ DoOpenError("Error occured while loading TPC tracks.");
+ return 1;
+ }
+ fCurrentEvent = 0;
return 0;
- }
+}
+/********************************************************************/
+void AliHBTReaderTPC::DoOpenError( const char *va_(fmt), ...)
+{
+ // Does error display and clean-up in case error caught on Open Next Session
+
+ va_list ap;
+ va_start(ap,va_(fmt));
+ Error("OpenNextSession", va_(fmt), ap);
+ va_end(ap);
+
+ delete fRunLoader;
+ fRunLoader = 0x0;
+ fTPCLoader = 0x0;
+ fCurrentDir++;
+}
/********************************************************************/
/********************************************************************/
#include <TString.h>
class TFile;
+class AliRunLoader;
+class AliTPCLoader;
class AliHBTReaderTPC: public AliHBTReader
{
virtual ~AliHBTReaderTPC();
- Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
+ void Rewind();
- AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
- AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles
- Int_t GetNumberOfPartEvents();//returns number of particle events
- Int_t GetNumberOfTrackEvents();//returns number of track events
+ Bool_t ReadsTracks() const {return kTRUE;}
+ Bool_t ReadsParticles() const {return kTRUE;}
+
+ void SetMagneticField(Float_t mf){fMagneticField=mf;}
+ void UseMagneticFieldFromRun(Bool_t flag = kTRUE){fUseMagFFromRun=flag;}
protected:
//in the future this class is will read global tracking
+ Int_t ReadNext();
+ Int_t OpenNextSession();
+ void DoOpenError(const char* msgfmt, ...);
- AliHBTRun* fParticles; //!simulated particles
- AliHBTRun* fTracks; //!reconstructed tracks (particles)
-
- TString fFileName;//name of the file with galice.root
+ TString fFileName;//name of the file with galice.root
+ AliRunLoader* fRunLoader;//!RL
+ AliTPCLoader* fTPCLoader;//!TPCLoader
+ Float_t fMagneticField;//magnetic field value that was enforced while reading
+ Bool_t fUseMagFFromRun;//flag indicating if using field specified in gAlice (kTRUE)
+ // or enforece other defined by fMagneticField
- Bool_t fIsRead;//!flag indicating if the data are already read
private:
public:
- ClassDef(AliHBTReaderTPC,2)
+ ClassDef(AliHBTReaderTPC,3)
};
if(!GetEvent(event)) fEvents->AddAtAndExpand(new AliHBTEvent, event);
GetEvent(event)->AddParticle(pdg,idx,px,py,pz,etot,vx,vy,vz,time);
}
+/**************************************************************************/
+void AliHBTRun::SetEvent(Int_t number, AliHBTEvent* event)
+{
+ //adds an event to the run
+ // makes an own copy of the event!
+ if (event == 0x0)
+ {
+ delete fEvents->RemoveAt(number);
+ return;
+ }
+ AliHBTEvent* ev = GetEvent(number);
+ if (ev) *ev = *event;
+ else fEvents->AddAtAndExpand(new AliHBTEvent(*event), number);
+}
/**************************************************************************/
+
Double_t px, Double_t py, Double_t pz, Double_t etot,
Double_t vx, Double_t vy, Double_t vz, Double_t time);
+ void SetEvent(Int_t number, AliHBTEvent* event);
AliHBTParticle* GetParticle(Int_t event, Int_t n); //returns nth particle from event
AliHBTEvent* GetEvent(Int_t event) const; //returns AliHBTEvent number "event"
AliHBTTwoTrackEffFctn3D();
virtual ~AliHBTTwoTrackEffFctn3D(){}
- void ProcessSameEventParticles(AliHBTPair* pair){}
- void ProcessDiffEventParticles(AliHBTPair* pair){}
+ void ProcessSameEventParticles(AliHBTPair* /*pair*/){}
+ void ProcessDiffEventParticles(AliHBTPair* /*pair*/){}
protected:
void GetValues(AliHBTPair*,Double_t&, Double_t&,Double_t&);
AliHBTReader* reader;
Int_t kine = strcmp(datatype,"Kine");
+ Int_t ESD = strcmp(datatype,"ESD");
Int_t TPC = strcmp(datatype,"TPC");
Int_t ITSv1 = strcmp(datatype,"ITSv1");
Int_t ITSv2 = strcmp(datatype,"ITSv2");
reader = new AliHBTReaderKineTree();
processopt="Particles"; //this reader by definition reads only simulated particles
}
+ else if(!ESD)
+ {
+ AliHBTReaderESD* esdreader = new AliHBTReaderESD();
+ esdreader->ReadParticles(kTRUE);
+ reader = esdreader;
+ }
else if(!TPC)
{
cout<<"AliHBTWriteInternFormat.C: Creating Reader TPC .....\n";
- reader = new AliHBTReaderTPC("tpc.tracks.root","tpc.clusters.root");
+ reader = new AliHBTReaderTPC();
cout<<"AliHBTWriteInternFormat.C: ..... Created\n";
}
else if(!ITSv1)
else if(!ITSv2)
{
cout<<"AliHBTWriteInternFormat.C: Creating Reader ITSv2 .....\n";
- reader = new AliHBTReaderITSv2("its.tracksv2.root","its.clustersv2.root");
+ reader = new AliHBTReaderITSv2();
cout<<"AliHBTWriteInternFormat.C: ..... Created\n";
}
else if(!intern)
{
- reader = new AliHBTReaderInternal("TPC.root");
+ reader = new AliHBTReaderInternal("data.root");
}
else
{
cout<<"AliHBTWriteInternFormat.C: P R O C S E S S I N G .....\n\n";
AliHBTReaderInternal::Write(reader,outfile);
- cout<<"\n\nAliHBTWriteInternFormat.C: F I N I S H E D";
+ cout<<"\n\nAliHBTWriteInternFormat.C: F I N I S H E D\n";
if (dirs) delete dirs;
delete reader;
#pragma link C++ class AliHBTParticle+;
#pragma link C++ class AliHBTEvent+;
#pragma link C++ class AliHBTRun+;
+#pragma link C++ class AliHBTEventBuffer+;
#pragma link C++ class AliHBTFunction+;
#pragma link C++ class AliHBTMonitorFunction+;
AliHBTMonDistributionFctns.cxx AliHBTMonResolutionFctns.cxx \
AliHBTLLWeights.cxx AliHBTLLWeightFctn.cxx \
AliHBTLLWeightsPID.cxx AliHBTLLWeightTheorFctn.cxx \
-AliHBTPositionRandomizer.cxx
+AliHBTPositionRandomizer.cxx AliHBTEventBuffer.cxx
FSRCS = fsiini.F fsiw.F led_bldata.F ltran12.F
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff