Obsolete code removed.
[u/mrichter/AliRoot.git] / JETAN / AliJetKineReader.cxx
index 363650c..a957076 100644 (file)
 #include <TPDGCode.h>
 #include <TParticle.h>
 #include <TParticlePDG.h>
-//#include <TVector3.h>
 #include <TLorentzVector.h>
 #include <TSystem.h>
 #include <TDatabasePDG.h>
 #include <TRandom.h>
 // From AliRoot ...
+#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()
+AliJetKineReader::AliJetKineReader():
+    AliJetReader(),  
+    fAliHeader(0),
+    fMCEvent(0),
+    fGenJets(0)
 {
-  // Constructor
-  fReaderHeader = 0x0;
-  fRunLoader    = 0x0;
-  fMass         = 0;
-  fPdgC         = 0;
+  // Default constructor
 }
 
 //____________________________________________________________________________
@@ -51,7 +56,10 @@ AliJetKineReader::AliJetKineReader()
 AliJetKineReader::~AliJetKineReader()
 {
   // Destructor
-  delete fReaderHeader;
+    if (fAliHeader) {
+       delete fAliHeader;
+       fAliHeader = 0;
+    }
 }
 
 //____________________________________________________________________________
@@ -59,123 +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   
 }
 
 //____________________________________________________________________________
 
-void AliJetKineReader::FillMomentumArray(Int_t event)
+Bool_t AliJetKineReader::FillMomentumArray()
 {
   // Fill momentum array for event
-  Int_t goodTrack = 0;
-  // Clear array
-  ClearArray();
-  // Get event from runloader
-  fRunLoader->GetEvent(event);
-  // 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();
-  
-  // temporary storage of signal and cut flags
-  Int_t* sflag  = new Int_t[nt];
-  Int_t* cflag  = new Int_t[nt];
+    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];
   
-  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; 
+    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();
+      
     
-    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
-      // Neutral particles (exclude K0L, n, nbar)
-      if (pdg == kNeutron || pdg == kK0Long) continue;
-    } // End fast EMCAL
-    
-    fMass = part->GetCalcMass();
-    fPdgC = part->GetPdgCode();
-
-    // 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++;
-  }
+    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
+         
+         // 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
+             // Neutral particles (exclude K0L, n, nbar)
+           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();
+       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);
+       }
+       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;
 }
 
 
@@ -235,3 +241,35 @@ Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
 
 }
 
+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 0;
+    AliGenEventHeader * genh = alih->GenEventHeader();
+    if (genh == 0) return 0;
+
+    Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets(); 
+    for (Int_t i = 0; i < nj; i++) {
+       Float_t p[4];
+       ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
+       new ((*fGenJets)[i]) AliAODJet(p[0], p[1], p[2], p[3]);
+    }
+
+    return fGenJets;
+}
+
+void AliJetKineReader::SetInputEvent(TObject* /*esd*/, TObject* /*aod*/, TObject* mc)
+{
+    // Connect the input event
+    fMCEvent = (AliMCEvent*) mc;
+}