Possibility to read separate input streams for signal and background.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Jan 2007 12:39:58 +0000 (12:39 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Jan 2007 12:39:58 +0000 (12:39 +0000)
JETAN/AliJetKineReader.cxx
JETAN/AliJetReaderHeader.cxx
JETAN/AliJetReaderHeader.h

index 55bf580..fd04b6b 100644 (file)
@@ -24,7 +24,6 @@
 #include <TPDGCode.h>
 #include <TParticle.h>
 #include <TParticlePDG.h>
-//#include <TVector3.h>
 #include <TLorentzVector.h>
 #include <TSystem.h>
 #include <TDatabasePDG.h>
@@ -86,88 +85,116 @@ void AliJetKineReader::OpenInputFiles()
 Bool_t AliJetKineReader::FillMomentumArray(Int_t event)
 {
   // Fill momentum array for event
+  char path[256];
+  const char* dirName   = fReaderHeader->GetDirectory();
+  const char* bgdirName = fReaderHeader->GetBgDirectory();
   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* sflag  = new Int_t[50000];
+  Int_t* cflag  = new Int_t[50000];
+  Int_t  spb = 0;
   
-  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();
-    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++;
+  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();
+      
+      
+      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();
+         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++;
+      }
   }
 
   // set the signal flags
index 796437d..f172d84 100755 (executable)
@@ -29,11 +29,14 @@ AliJetReaderHeader::AliJetReaderHeader():
  TNamed("AliJetReaderHeader", "Jet Reader Header"),
  fFirst(0),
  fLast(-1),
+ fOption(0),
+ fSignalPerBg(0),
  fFiducialEtaMin(-0.9),
  fFiducialEtaMax(0.9),
  fPtCut(2.0),
  fComment("No comment"),
  fDir(""),
+ fBgDir(""),
  fPattern("")
 {
   // Default constructor
@@ -45,11 +48,14 @@ AliJetReaderHeader::AliJetReaderHeader(const char * name):
  TNamed(name, "Jet Reader Header"),
  fFirst(0),
  fLast(-1),
+ fOption(0),
+ fSignalPerBg(0),
  fFiducialEtaMin(-0.9),
  fFiducialEtaMax(0.9),
  fPtCut(2.0),
  fComment("No comment"),
  fDir(""),
+ fBgDir(""),
  fPattern("")
 {
   // Constructor
index 4d7c6f0..7c87dac 100755 (executable)
@@ -22,21 +22,26 @@ class AliJetReaderHeader : public TNamed
   virtual ~AliJetReaderHeader();
   
   // Getters
-  virtual const TString GetComment() {return fComment;}
-  virtual const char* GetDirectory() {return fDir.Data();}
+  virtual const TString GetComment()  {return fComment;}
+  virtual const char* GetDirectory()  {return fDir.Data();}
+  virtual const char* GetBgDirectory(){return fBgDir.Data();}
   virtual const char* GetPattern() {return fPattern.Data();}
   virtual Float_t     GetFiducialEtaMin() const {return fFiducialEtaMin;}
   virtual Float_t     GetFiducialEtaMax() const {return fFiducialEtaMax;}  
   virtual Float_t     GetPtCut() const {return fPtCut;}
-  Int_t   GetNEvents()    const {return fLast-fFirst;}
-  Int_t   GetFirstEvent() const {return fFirst;}
-  Int_t   GetLastEvent()  const {return fLast;}
-  Int_t   GetDetector()   const {return fOption;}
-  Int_t   GetDebug()      const {return fDebug;}
+  Int_t   GetNEvents()     const {return fLast-fFirst;}
+  Int_t   GetFirstEvent()  const {return fFirst;}
+  Int_t   GetLastEvent()   const {return fLast;}
+  Int_t   GetDetector()    const {return fOption;}
+  Int_t   GetDebug()       const {return fDebug;}
+  Int_t   GetSignalPerBg() const {return fSignalPerBg;}
+         
   // Setters
-  virtual void SetComment(const char* s) {fComment=TString(s);}
-  virtual void SetPattern(const char* s) {fPattern=TString(s);}
-  virtual void SetDirectory(const char* s) {fDir=TString(s);}
+  virtual void SetComment(const char* s)     {fComment=TString(s);}
+  virtual void SetPattern(const char* s)     {fPattern=TString(s);}
+  virtual void SetDirectory(const char* s)   {fDir=TString(s);}
+  virtual void SetBgDirectory(const char* s, Int_t n = 1)
+      {fDir=TString(s); fSignalPerBg = n;}
   virtual void SetFirstEvent(Int_t i=0) {fFirst=i;}
   virtual void SetLastEvent(Int_t i=-1) {fLast=i;}
   virtual void SetFiducialEta(Float_t etamin, Float_t etamax) 
@@ -50,11 +55,13 @@ class AliJetReaderHeader : public TNamed
   Int_t   fLast;           // in current set of files
   Int_t   fOption;         // detector used for jet reconstruction   
   Int_t   fDebug;          // debug option
+  Int_t   fSignalPerBg;
   Float_t fFiducialEtaMin; // Fiducial minimum eta
   Float_t fFiducialEtaMax; // Fiducial maximum eta
   Float_t fPtCut;          // pt cut
   TString fComment;        // a comment
-  TString fDir;            // directory with input files
+  TString fDir;            // directory with input files for signal
+  TString fBgDir;          // directory with input files for background
   TString fPattern;        // pattern to look for input files
   
   ClassDef(AliJetReaderHeader,2);