]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HBTAN/AliHBTReaderInternal.cxx
Bug corrections
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderInternal.cxx
index 333828278e60e1f599d4d802632ab16fd3fc506d..8480c4201e16f801af4c369f9d41ab97528f9705 100644 (file)
@@ -1,4 +1,19 @@
 #include "AliHBTReaderInternal.h"
+//_________________________________________________________________________
+///////////////////////////////////////////////////////////////////////////
+//                                                                       //
+// class AliHBTReaderInternal                                            //
+//                                                                       //
+// Multi file reader for Internal Data Format                            //
+//                                                                       //
+// This reader reads data created by itself                              //
+//   (method AliHBTReaderInternal::Write)                                //
+// Data are stored in form of tree of TClonesArray of AliHBTParticle's   //
+//                                                                       //
+// Piotr.Skowronski@cern.ch                                              //
+// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html  //
+//                                                                       //
+///////////////////////////////////////////////////////////////////////////
 
 #include <TTree.h>
 #include <TFile.h>
@@ -12,8 +27,6 @@
 #include "AliHBTParticle.h"
 #include "AliHBTParticleCut.h"
 
-AliHBTReaderInternal test;
-
 ClassImp(AliHBTReaderInternal)
 /********************************************************************/
 
@@ -22,7 +35,9 @@ AliHBTReaderInternal::AliHBTReaderInternal():
  fPartBranch(0x0),
  fTrackBranch(0x0),
  fTree(0x0),
- fFile(0x0)
+ fFile(0x0),
+ fPartBuffer(0x0),
+ fTrackBuffer(0x0)
 {
 //Defalut constructor
 }
@@ -33,7 +48,9 @@ AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):
  fPartBranch(0x0),
  fTrackBranch(0x0),
  fTree(0x0),
- fFile(0x0)
+ fFile(0x0),
+ fPartBuffer(0x0),
+ fTrackBuffer(0x0)
 { 
 //constructor 
 //filename - name of file to open
@@ -46,7 +63,9 @@ AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename
  fPartBranch(0x0),
  fTrackBranch(0x0),
  fTree(0x0),
- fFile(0x0)
+ fFile(0x0),
+ fPartBuffer(0x0),
+ fTrackBuffer(0x0)
 { 
 //ctor
 //dirs contains strings with directories to look data in
@@ -54,6 +73,30 @@ AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename
 }
 /********************************************************************/
 
+AliHBTReaderInternal::AliHBTReaderInternal(const AliHBTReaderInternal& in):
+ AliHBTReader(in),
+ fFileName(in.fFileName),
+ fPartBranch(0x0),
+ fTrackBranch(0x0),
+ fTree(0x0),
+ fFile(0x0),
+ fPartBuffer(0x0),
+ fTrackBuffer(0x0)
+{
+  //cpy constructor
+}
+/********************************************************************/
+
+AliHBTReaderInternal& AliHBTReaderInternal::operator=(const AliHBTReaderInternal& in)
+{
+  //Assigment operator
+  if (this == &in) return *this;
+  Rewind();//close current session
+  AliHBTReader::operator=((const AliHBTReader&)in);
+  fFileName = in.fFileName;
+  return *this;
+}
+/********************************************************************/
 AliHBTReaderInternal::~AliHBTReaderInternal()
  {
  //desctructor
@@ -126,33 +169,70 @@ Int_t AliHBTReaderInternal::ReadNext()
 
     counter = 0;  
     if (fPartBranch && fTrackBranch)
-    {
+     {
+       Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());        
+       Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
        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 == 0x0 ) continue; //if returned pointer is NULL
+           
+           if (AliHBTParticle::GetDebug() > 9)
+            {
+              Info("ReadNext","Particle:");
+              tpart->Print();
+              Info("ReadNext","Track:");
+              ttrack->Print();
+            }
            if (ttrack->GetUID() != tpart->GetUID())
              {
                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++)
+           
+           for (Int_t s = 0; s < ttrack->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) )
+              //check if we are intersted with particles of this type
+              //if not take next partilce
+              if( Pass(ttrack->GetNthPid(s)) ) 
+               {
+                 if (AliHBTParticle::GetDebug() > 9)
+                  Info("ReadNext","Track Incarnation %d did not pass PID cut.",ttrack->GetNthPid(s));
+                 continue; 
+               }
+              TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
+              if (pdgp == 0x0)//PDG part corresponding to new incarnation
+               {
+                 Error("ReadNext","Particle code unknown to PDG DB.");
+                 continue;
+               }
+              
+              AliHBTParticle* track = new AliHBTParticle(*ttrack);
+              
+              //apart of setting PDG code of an incarnation
+              //it is necessary tu recalculate energy on the basis of
+              //new PDG code (mass) hypothesis
+              Double_t mass = pdgp->Mass();//mass of new incarnation
+              Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() + 
+                                            ttrack->Py()*ttrack->Py() + 
+                                            ttrack->Pz()*ttrack->Pz() + 
+                                            mass*mass);//total energy of the new incarnation
+              //update Energy and Calc Mass 
+              track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
+              track->SetCalcMass(mass);
+              track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
+              
+              if( Pass(track) )
                 {
-                  delete part;
+                  if (AliHBTParticle::GetDebug() > 9)
+                   Info("ReadNext","Track Incarnation %d did not pass cut.",ttrack->GetNthPid(s));
+                  delete track;
                   continue; 
                 }
-              AliHBTParticle* track = new AliHBTParticle(*ttrack);
+              AliHBTParticle* part = new AliHBTParticle(*tpart);//particle has only 1 incarnation (real)
 
               fParticlesEvent->AddParticle(part);
               fTracksEvent->AddParticle(track);
@@ -176,7 +256,6 @@ Int_t AliHBTReaderInternal::ReadNext()
             {
               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) )
@@ -194,23 +273,44 @@ Int_t AliHBTReaderInternal::ReadNext()
 
      if (fTrackBranch)
       {
+        Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());       
         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
+           ttrack =  dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
+           if( ttrack == 0x0 ) continue; //if returned pointer is NULL
 
-              AliHBTParticle* part = new AliHBTParticle(*tpart);
-              part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
-              if( Pass(part) )
+           for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
+            {
+              if( Pass(ttrack->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
+                                                         //if not take next partilce
+              TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
+              if (pdgp == 0x0)//PDG part corresponding to new incarnation
+               {
+                 Error("ReadNext","Particle code unknown to PDG DB.");
+                 continue;
+               }
+              AliHBTParticle* track = new AliHBTParticle(*ttrack);
+              
+              //apart of setting PDG code of an incarnation
+              //it is necessary tu recalculate energy on the basis of
+              //new PDG code (mass) hypothesis
+              Double_t mass = pdgp->Mass();//mass of new incarnation
+              Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() + 
+                                            ttrack->Py()*ttrack->Py() + 
+                                            ttrack->Pz()*ttrack->Pz() + 
+                                            mass*mass);//total energy of the new incarnation
+              //update Energy and Calc Mass 
+              track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
+              track->SetCalcMass(mass);
+              track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
+              
+              if( Pass(track) )
                 {
-                  delete part;
+                  delete track;
                   continue; 
                 }
-              fTracksEvent->AddParticle(part);
+              fTracksEvent->AddParticle(track);
+
               counter++;
             }
          }
@@ -290,7 +390,7 @@ Int_t AliHBTReaderInternal::OpenNextFile()
 
    Info("OpenNextFile","________________________________________________________");
    Info("OpenNextFile","Found %d event(s) in directory %s",
-         (Int_t)fTree->GetEntries(),GetDirName(fCurrentEvent).Data());
+         (Int_t)fTree->GetEntries(),GetDirName(fCurrentDir).Data());
    
    fCurrentEvent = 0;