JETAN/AliJetKineReader.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-------------------------------------------------------------------------
  17 // Jet kinematics reader
  18 // MC reader for jet analysis
  19 // Author: Andreas Morsch (andreas.morsch@cern.ch)
  20 //-------------------------------------------------------------------------
  21
  22 // From root ...
  23 #include <TClonesArray.h>
  24 #include <TPDGCode.h>
  25 #include <TParticle.h>
  26 #include <TParticlePDG.h>
  27 #include <TLorentzVector.h>
  28 #include <TSystem.h>
  29 #include <TDatabasePDG.h>
  30 #include <TRandom.h>
  31 // From AliRoot ...
  32 #include "AliJet.h"
  33 #include "AliJetKineReaderHeader.h"
  34 #include "AliJetKineReader.h"
  35 #include "AliRunLoader.h"
  36 #include "AliStack.h"
  37 #include "AliHeader.h"
  38 #include "AliGenPythiaEventHeader.h"
  39
  40 ClassImp(AliJetKineReader)
  41
  42 AliJetKineReader::AliJetKineReader():
  43     AliJetReader(),
  44     fRunLoader(0x0),
  45     fAliHeader(0x0)
  46 {
  47   // Default constructor
  48 }
  49
  50 //____________________________________________________________________________
  51
  52 AliJetKineReader::~AliJetKineReader()
  53 {
  54   // Destructor
  55   delete fAliHeader;
  56 }
  57
  58 //____________________________________________________________________________
  59
  60 void AliJetKineReader::OpenInputFiles()
  61 {
  62   // Opens the input file using the run loader
  63   const char* dirName = fReaderHeader->GetDirectory();
  64   char path[256];
  65   sprintf(path, "%s/galice.root",dirName);
  66   fRunLoader = AliRunLoader::Open(path);
  67   fRunLoader->LoadKinematics();
  68   fRunLoader->LoadHeader();
  69
  70   // get number of events
  71   Int_t nMax = fRunLoader->GetNumberOfEvents();
  72   printf("\nTotal number of events = %d", nMax);
  73
  74   // set number of events in header
  75   if (fReaderHeader->GetLastEvent() == -1)
  76     fReaderHeader->SetLastEvent(nMax);
  77   else {
  78     Int_t nUsr = fReaderHeader->GetLastEvent();
  79     fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
  80   }
  81 }
  82
  83 //____________________________________________________________________________
  84
  85 Bool_t AliJetKineReader::FillMomentumArray(Int_t event)
  86 {
  87   // Fill momentum array for event
  88   char path[256];
  89   const char* dirName   = fReaderHeader->GetDirectory();
  90   const char* bgdirName = fReaderHeader->GetBgDirectory();
  91   Int_t goodTrack = 0;
  92   // Clear array
  93   // temporary storage of signal and cut flags
  94   Int_t* sflag  = new Int_t[50000];
  95   Int_t* cflag  = new Int_t[50000];
  96   Int_t  spb = 0;
  97
  98   ClearArray();
  99   for (Int_t i = 0; i < 2; i++) {
 100       if (i == 1) {
 101 // Pythia Events
 102           if (fRunLoader) delete fRunLoader;
 103           sprintf(path, "%s/galice.root", dirName);
 104           fRunLoader = AliRunLoader::Open(path);
 105           fRunLoader->LoadKinematics();
 106           fRunLoader->LoadHeader();
 107           // Get event from runloader
 108           fRunLoader->GetEvent(event);
 109       } else {
 110 // Background events
 111           if ((spb = fReaderHeader->GetSignalPerBg()) > 0) {
 112               if (fRunLoader) delete fRunLoader;
 113               sprintf(path, "%s/galice.root", bgdirName);
 114               fRunLoader = AliRunLoader::Open(path);
 115               fRunLoader->LoadKinematics();
 116               fRunLoader->LoadHeader();
 117               // Get event from runloader
 118               fRunLoader->GetEvent(event / spb);
 119           } else {
 120               continue;
 121           }
 122       }
 123
 124       // Get the stack
 125       AliStack* stack = fRunLoader->Stack();
 126       // Number of primaries
 127       Int_t nt = stack->GetNprimary();
 128
 129       // Get cuts set by user and header
 130       Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
 131       Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 132       Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 133       fAliHeader = fRunLoader->GetHeader();
 134
 135
 136       TLorentzVector p4;
 137       // Loop over particles
 138       for (Int_t it = 0; it < nt; it++) {
 139           TParticle *part = stack->Particle(it);
 140           Int_t   status  = part->GetStatusCode();
 141           Int_t   pdg     = TMath::Abs(part->GetPdgCode());
 142           Float_t pt      = part->Pt();
 143
 144           // Skip non-final state particles, neutrinos
 145           // and particles with pt < pt_min
 146           if ((status != 1)
 147               || (pdg == 12 || pdg == 14 || pdg == 16)) continue;
 148
 149           Float_t p       = part->P();
 150           Float_t p0      = p;
 151           Float_t eta     = part->Eta();
 152           Float_t phi     = part->Phi();
 153
 154           // Fast simulation of TPC if requested
 155           if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
 156               // Charged particles only
 157               Float_t charge =
 158                   TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
 159               if (charge == 0)               continue;
 160               // Simulate efficiency
 161               if (!Efficiency(p0, eta, phi)) continue;
 162               // Simulate resolution
 163               p = SmearMomentum(4, p0);
 164           } // End fast TPC
 165
 166           // Fast simulation of EMCAL if requested
 167           if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
 168               Float_t charge =
 169                   TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
 170               // Charged particles only
 171               if (charge != 0){
 172                   // Simulate efficiency
 173                   if (!Efficiency(p0, eta, phi)) continue;
 174                   // Simulate resolution
 175                   p = SmearMomentum(4, p0);
 176               } // end "if" charged particles
 177               // Neutral particles (exclude K0L, n, nbar)
 178               if (pdg == kNeutron || pdg == kK0Long) continue;
 179           } // End fast EMCAL
 180
 181           // Fill momentum array
 182           Float_t r  = p/p0;
 183           Float_t px = r * part->Px();
 184           Float_t py = r * part->Py();
 185           Float_t pz = r * part->Pz();
 186           Float_t m  = part->GetMass();
 187           Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
 188           p4 = TLorentzVector(px, py, pz, e);
 189           if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
 190           new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
 191
 192           // all tracks are signal ... ?
 193           sflag[goodTrack]=1;
 194           cflag[goodTrack]=0;
 195           if (pt > ptMin) cflag[goodTrack]=1; // track surviving pt cut
 196           goodTrack++;
 197       }
 198   }
 199
 200   // set the signal flags
 201   fSignalFlag.Set(goodTrack, sflag);
 202   fCutFlag.Set(goodTrack, cflag);
 203   printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
 204
 205   delete[] sflag;
 206   delete[] cflag;
 207
 208   return kTRUE;
 209 }
 210
 211
 212 Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
 213 {
 214   //  The relative momentum error, i.e.
 215   //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
 216   //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
 217   //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
 218   //
 219   //  If we include information from TRD, b will be a factor 2/3 smaller.
 220   //
 221   //  ind = 1: high multiplicity
 222   //  ind = 2: low  multiplicity
 223   //  ind = 3: high multiplicity + TRD
 224   //  ind = 4: low  multiplicity + TRD
 225
 226   Float_t pSmeared;
 227   Float_t a = 0.75;
 228   Float_t b = 0.12;
 229
 230   if (ind == 1) b = 0.12;
 231   if (ind == 2) b = 0.08;
 232   if (ind == 3) b = 0.12;
 233   if (ind == 4) b = 0.08;
 234
 235   Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
 236   pSmeared = p + gRandom->Gaus(0., sigma);
 237   return pSmeared;
 238 }
 239
 240
 241 Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
 242 {
 243   // Fast simulation of geometrical acceptance and tracking efficiency
 244
 245   //  Tracking
 246   Float_t eff = 0.99;
 247   if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
 248   // Geometry
 249   if (p > 0.5) {
 250     phi *= 180. / TMath::Pi();
 251     // Sector number 0 - 17
 252     Int_t isec  = Int_t(phi / 20.);
 253     // Sector centre
 254     Float_t phi0 = isec * 20. + 10.;
 255     Float_t phir = TMath::Abs(phi-phi0);
 256     // 2 deg of dead space
 257     if (phir > 9.) eff = 0.;
 258   }
 259
 260   if (gRandom->Rndm() > eff) {
 261     return kFALSE;
 262   } else {
 263     return kTRUE;
 264   }
 265
 266 }
 267
 268 Bool_t AliJetKineReader::GetGenJets(AliJet* genJets)
 269 {
 270     // Get generated jets from mc header
 271     AliHeader* alih = GetAliHeader();
 272     if (alih == 0) return kFALSE;
 273     AliGenEventHeader * genh = alih->GenEventHeader();
 274     if (genh == 0) return kFALSE;
 275     Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets();
 276     Int_t* m = new Int_t[nj];
 277     Int_t* k = new Int_t[nj];
 278     for (Int_t i=0; i< nj; i++) {
 279         Float_t p[4];
 280         ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
 281         genJets->AddJet(p[0],p[1],p[2],p[3]);
 282         m[i]=1;
 283         k[i]=i;
 284     }
 285     genJets->SetNinput(nj);
 286     genJets->SetMultiplicities(m);
 287     genJets->SetInJet(k);
 288     delete[] k;
 289     delete[] m;
 290     return kTRUE;
 291 }