Obsolete code removed.
[u/mrichter/AliRoot.git] / JETAN / AliJetKineReader.cxx
index 13280ea..a957076 100644 (file)
 #include <TDatabasePDG.h>
 #include <TRandom.h>
 // From AliRoot ...
-#include "AliJet.h"
+#include "AliAODJet.h"
+#include "AliPDG.h"
 #include "AliJetKineReaderHeader.h"
 #include "AliJetKineReader.h"
-#include "AliRunLoader.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
 #include "AliStack.h"
 #include "AliHeader.h"
 #include "AliGenPythiaEventHeader.h"
+#include "AliLog.h"
 
 ClassImp(AliJetKineReader)
 
 AliJetKineReader::AliJetKineReader():
     AliJetReader(),  
-    fRunLoader(0x0),
-    fAliHeader(0x0)
+    fAliHeader(0),
+    fMCEvent(0),
+    fGenJets(0)
 {
   // Default constructor
 }
@@ -52,7 +56,10 @@ AliJetKineReader::AliJetKineReader():
 AliJetKineReader::~AliJetKineReader()
 {
   // Destructor
-  delete fAliHeader;
+    if (fAliHeader) {
+       delete fAliHeader;
+       fAliHeader = 0;
+    }
 }
 
 //____________________________________________________________________________
@@ -60,152 +67,121 @@ AliJetKineReader::~AliJetKineReader()
 void AliJetKineReader::OpenInputFiles()
 {
   // Opens the input file using the run loader
-  const char* dirName = fReaderHeader->GetDirectory();
-  char path[256];
-  sprintf(path, "%s/galice.root",dirName);
-  fRunLoader = AliRunLoader::Open(path);
-  fRunLoader->LoadKinematics();
-  fRunLoader->LoadHeader(); 
-  
-  // get number of events
-  Int_t nMax = fRunLoader->GetNumberOfEvents();
-  printf("\nTotal number of events = %d", nMax);
-  
-  // set number of events in header
-  if (fReaderHeader->GetLastEvent() == -1)
-    fReaderHeader->SetLastEvent(nMax);
-  else {
-    Int_t nUsr = fReaderHeader->GetLastEvent();
-    fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
-  }
+  // Not used anymore   
 }
 
 //____________________________________________________________________________
 
-Bool_t AliJetKineReader::FillMomentumArray(Int_t event)
+Bool_t AliJetKineReader::FillMomentumArray()
 {
   // Fill momentum array for event
-  char path[256];
-  const char* dirName   = fReaderHeader->GetDirectory();
-  const char* bgdirName = fReaderHeader->GetBgDirectory();
-  Int_t goodTrack = 0;
-  // Clear array
-  // temporary storage of signal and cut flags
-  Int_t* sflag  = new Int_t[50000];
-  Int_t* cflag  = new Int_t[50000];
-  Int_t  spb = 0;
+    Int_t goodTrack = 0;
+    // Clear array
+    // temporary storage of signal and cut flags
+    Int_t* sflag  = new Int_t[50000];
+    Int_t* cflag  = new Int_t[50000];
   
-  ClearArray();
-  for (Int_t i = 0; i < 2; i++) {
-      if (i == 1) {
-// Pythia Events
-         if (fRunLoader) delete fRunLoader;
-         sprintf(path, "%s/galice.root", dirName);
-         fRunLoader = AliRunLoader::Open(path);
-         fRunLoader->LoadKinematics();
-         fRunLoader->LoadHeader(); 
-         // Get event from runloader
-         fRunLoader->GetEvent(event);
-      } else {
-// Background events
-         if ((spb = fReaderHeader->GetSignalPerBg()) > 0) {
-             if (fRunLoader) delete fRunLoader;
-             sprintf(path, "%s/galice.root", bgdirName);
-             fRunLoader = AliRunLoader::Open(path);
-             fRunLoader->LoadKinematics();
-             fRunLoader->LoadHeader(); 
-             // Get event from runloader
-             fRunLoader->GetEvent(event / spb);
-         } else {
-             continue;
-         }
-      }
-
-      // Get the stack
-      AliStack* stack = fRunLoader->Stack();
-      // Number of primaries
-      Int_t nt = stack->GetNprimary();
-      
-      // Get cuts set by user and header
-      Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
-      Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
-      Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
-      fAliHeader = fRunLoader->GetHeader();
+    ClearArray();
+    // Get the stack
+    AliStack* stack = fMCEvent->Stack();
+    // Number of primaries
+    Int_t nt = stack->GetNprimary();
       
+    // Get cuts set by user and header
+    Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
+    Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+    Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
+    fAliHeader = fMCEvent->Header();
       
-      TLorentzVector p4;
-      // Loop over particles
-      for (Int_t it = 0; it < nt; it++) {
-         TParticle *part = stack->Particle(it);
-         Int_t   status  = part->GetStatusCode();
-         Int_t   pdg     = TMath::Abs(part->GetPdgCode());
-         Float_t pt      = part->Pt(); 
+    
+    TLorentzVector p4;
+    // Loop over particles
+    for (Int_t it = 0; it < nt; it++) {
+       TParticle *part = stack->Particle(it);
+       Int_t   status  = part->GetStatusCode();
+       Int_t   pdg     = TMath::Abs(part->GetPdgCode());
+       Float_t pt      = part->Pt(); 
          
-         // Skip non-final state particles, neutrinos 
-         // and particles with pt < pt_min 
-         if ((status != 1)            
-             || (pdg == 12 || pdg == 14 || pdg == 16)) continue; 
-         
-         Float_t p       = part->P();
-         Float_t p0      = p;
-         Float_t eta     = part->Eta();
-         Float_t phi     = part->Phi();
-         
-         // Fast simulation of TPC if requested
-         if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
-             // Charged particles only
-             Float_t charge = 
-                 TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
-             if (charge == 0)               continue;
-             // Simulate efficiency
-             if (!Efficiency(p0, eta, phi)) continue;
-             // Simulate resolution
-             p = SmearMomentum(4, p0);
-         } // End fast TPC
+       // Skip non-final state particles, neutrinos 
+       // and particles with pt < pt_min 
+       if ((status != 1)            
+           || (pdg == 12 || pdg == 14 || pdg == 16)) continue; 
+       
+       Float_t p       = part->P();
+       Float_t p0      = p;
+       Float_t eta     = part->Eta();
+       Float_t phi     = part->Phi();
+       
+       // Fast simulation of TPC if requested
+       if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
+           // Charged particles only
+           Float_t charge = 
+               TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+           if (charge == 0)               continue;
+           // Simulate efficiency
+           if (!Efficiency(p0, eta, phi)) continue;
+           // Simulate resolution
+           p = SmearMomentum(4, p0);
+       } // End fast TPC
          
          // Fast simulation of EMCAL if requested
-         if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
-             Float_t charge = 
-                 TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
-             // Charged particles only
-             if (charge != 0){
-                 // Simulate efficiency
-                 if (!Efficiency(p0, eta, phi)) continue;
-                 // Simulate resolution
-                 p = SmearMomentum(4, p0);
-             } // end "if" charged particles
+       if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
+           Float_t charge = 
+               TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+           // Charged particles only
+           if (charge != 0){
+               // Simulate efficiency
+               if (!Efficiency(p0, eta, phi)) continue;
+               // Simulate resolution
+               p = SmearMomentum(4, p0);
+           } // end "if" charged particles
              // Neutral particles (exclude K0L, n, nbar)
-             if (pdg == kNeutron || pdg == kK0Long) continue;
-         } // End fast EMCAL
-         
+           if (pdg == kNeutron || pdg == kK0Long) continue;
+       } // End fast EMCAL
+       
          // Fill momentum array
-         Float_t r  = p/p0;
-         Float_t px = r * part->Px(); 
-         Float_t py = r * part->Py(); 
-         Float_t pz = r * part->Pz();
-         Float_t m  = part->GetMass();
-         Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
-         p4 = TLorentzVector(px, py, pz, e);
-         if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
-         new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
-         
-         // all tracks are signal ... ?
-         sflag[goodTrack]=1;
-         cflag[goodTrack]=0;
-         if (pt > ptMin) cflag[goodTrack]=1; // track surviving pt cut
-         goodTrack++;
-      }
-  }
+       Float_t r  = p/p0;
+       Float_t px = r * part->Px(); 
+       Float_t py = r * part->Py(); 
+       Float_t pz = r * part->Pz();
+       TParticlePDG *pPDG = part->GetPDG();
+       Float_t m = 0;
+       if(!pPDG){
+         // this is very rare...
+         // Is avoided by AliPDG::AddParticlesToPdgDataBase();
+         // but this should be called only once... (how in proof?)
+         // Calucluate mass with unsmeared momentum values
+         m  = part->Energy()*part->Energy() - 
+           (px * px + py * py + pz * pz)/r/r; 
+         if(m>0)m = TMath::Sqrt(m);
+         else m = 0;
+         AliInfo(Form("Unknown Particle using %d calculated mass m = %3.3f",part->GetPdgCode(),m));
 
-  // set the signal flags
-  fSignalFlag.Set(goodTrack, sflag);
-  fCutFlag.Set(goodTrack, cflag);
-  printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
-  
-  delete[] sflag;
-  delete[] cflag;
-  
-  return kTRUE;
+       }
+       else{
+         m  = pPDG->Mass();
+       }
+       Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
+       p4 = TLorentzVector(px, py, pz, e);
+       if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
+       new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
+       
+       // all tracks are signal ... ?
+       sflag[goodTrack]=1;
+       cflag[goodTrack]=0;
+       if (pt > ptMin) cflag[goodTrack] = 1; // track surviving pt cut
+       goodTrack++;
+    } // track loop
+
+// set the signal flags
+    fSignalFlag.Set(goodTrack, sflag);
+    fCutFlag.Set(goodTrack, cflag);
+    AliInfo(Form(" Number of good tracks %d \n", goodTrack));
+    
+    delete[] sflag;
+    delete[] cflag;
+    
+    return kTRUE;
 }
 
 
@@ -265,27 +241,35 @@ Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
 
 }
 
-Bool_t AliJetKineReader::GetGenJets(AliJet* genJets)
+TClonesArray*  AliJetKineReader::GetGeneratedJets()
 {
+    //
     // Get generated jets from mc header
+    //
+    if (!fGenJets) {
+       fGenJets = new TClonesArray("AliAODJets", 100);
+    } else {
+       fGenJets->Clear();
+    }
+    
+    
     AliHeader* alih = GetAliHeader(); 
-    if (alih == 0) return kFALSE;
+    if (alih == 0) return 0;
     AliGenEventHeader * genh = alih->GenEventHeader();
-    if (genh == 0) return kFALSE;
+    if (genh == 0) return 0;
+
     Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets(); 
-    Int_t* m = new Int_t[nj];
-    Int_t* k = new Int_t[nj];
-    for (Int_t i=0; i< nj; i++) {
+    for (Int_t i = 0; i < nj; i++) {
        Float_t p[4];
        ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
-       genJets->AddJet(p[0],p[1],p[2],p[3]);
-       m[i]=1;
-       k[i]=i;
+       new ((*fGenJets)[i]) AliAODJet(p[0], p[1], p[2], p[3]);
     }
-    genJets->SetNinput(nj);
-    genJets->SetMultiplicities(m);
-    genJets->SetInJet(k);
-    delete[] k;
-    delete[] m;
-    return kTRUE;
+
+    return fGenJets;
+}
+
+void AliJetKineReader::SetInputEvent(TObject* /*esd*/, TObject* /*aod*/, TObject* mc)
+{
+    // Connect the input event
+    fMCEvent = (AliMCEvent*) mc;
 }