JETAN/AliJetParticlesReaderKine.cxx

   1 // $Id$
   2
   3 //_______________________________________________________________________
   4 /////////////////////////////////////////////////////////////////////////
   5 //
   6 // class AliJetParticlesReaderKine
   7 //
   8 // Reader for Kinematics
   9 //
  10 // loizides@ikf.uni-frankfurt.de
  11 //
  12 /////////////////////////////////////////////////////////////////////////
  13
  14 #include <Riostream.h>
  15 #include <TString.h>
  16 #include <TParticle.h>
  17 #include <AliRunLoader.h>
  18 #include <AliStack.h>
  19
  20 #include "AliJetParticle.h"
  21 #include "AliJetEventParticles.h"
  22 #include "AliJetParticlesReaderKine.h"
  23 #include "AliHeader.h"
  24 #include "AliGenPythiaEventHeader.h"
  25
  26 ClassImp(AliJetParticlesReaderKine)
  27
  28 AliJetParticlesReaderKine::AliJetParticlesReaderKine() :
  29   AliJetParticlesReader(),
  30   fFileName("galice.root"),
  31   fRunLoader(0),
  32   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
  33   fUseTracks(kFALSE)
  34 {
  35   //constructor
  36 }
  37
  38 AliJetParticlesReaderKine::AliJetParticlesReaderKine(TString& fname) :
  39   AliJetParticlesReader(),
  40   fFileName(fname),
  41   fRunLoader(0),
  42   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
  43   fUseTracks(kFALSE)
  44 {
  45   //constructor
  46 }
  47
  48 AliJetParticlesReaderKine::AliJetParticlesReaderKine(TObjArray* dirs, const Char_t *filename):
  49   AliJetParticlesReader(dirs),
  50   fFileName(filename),
  51   fRunLoader(0),
  52   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
  53   fUseTracks(kFALSE)
  54 {
  55   //constructor
  56 }
  57
  58 AliJetParticlesReaderKine::~AliJetParticlesReaderKine()
  59 {
  60   //destructor
  61   if(fRunLoader) delete fRunLoader;
  62 }
  63
  64 void AliJetParticlesReaderKine::Rewind()
  65 {
  66   //Rewinds to the beginning
  67   if(fRunLoader) delete fRunLoader;
  68   fRunLoader   = 0;
  69   fCurrentDir  = 0;
  70   fNEventsRead = 0;
  71 }
  72
  73 Int_t AliJetParticlesReaderKine::ReadNext()
  74 {
  75   //Reads Kinematics Tree
  76
  77   if((!fOwner) || (fEventParticles == 0))
  78     fEventParticles = new AliJetEventParticles();
  79
  80   while(fCurrentDir < GetNumberOfDirs())
  81     {
  82
  83       if (!OpenFile(fCurrentDir))
  84       {
  85         delete fRunLoader; //close current session
  86         fRunLoader = 0;    //assure pointer is null
  87         fCurrentDir++;
  88         continue;
  89       }
  90
  91       if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
  92         {
  93           //read next directory
  94           delete fRunLoader; //close current session
  95           fRunLoader = 0;    //assure pointer is null
  96           fCurrentDir++;     //go to next dir
  97           continue;
  98         }
  99
 100       Info("ReadNext","Reading Event %d",fCurrentEvent);
 101       fRunLoader->GetEvent(fCurrentEvent);
 102       AliStack* stack = fRunLoader->Stack();
 103       if (!stack)
 104         {
 105           Error("ReadNext","Can't get stack for event %d",fCurrentEvent);
 106           continue;
 107         }
 108
 109       //clear old event
 110       Int_t nprim = stack->GetNprimary();
 111       Int_t npart = nprim;
 112       if(fUseTracks)
 113         npart = stack->GetNtrack();
 114       fEventParticles->Reset(npart);
 115
 116       TString headdesc="";
 117       AliHeader *header=fRunLoader->GetHeader();
 118       if(!header) {
 119         Warning("ReadNext","Header not found in event %d",fCurrentEvent);
 120       } else {
 121         headdesc+="Run ";
 122         headdesc+=header->GetRun();
 123         headdesc+=": Ev ";
 124         headdesc+=header->GetEventNrInRun();
 125         AliGenPythiaEventHeader *pheader=(AliGenPythiaEventHeader*)header->GenEventHeader();
 126         if(!pheader) {
 127           Warning("ReadNext","Pythia-Header not found in event %d",fCurrentEvent);
 128         } else {
 129           Int_t ntruq=0;
 130 #ifndef NOUQHEADERINFO
 131           ntruq=pheader->NUQTriggerJets();
 132           if(ntruq){
 133             Double_t x0=pheader->GetXJet();
 134             Double_t y0=pheader->GetYJet();
 135             Double_t zquench[4];
 136             pheader->GetZQuench(zquench);
 137             fEventParticles->SetXYJet(x0,y0);
 138             fEventParticles->SetZQuench(zquench);
 139             for(Int_t j=0;j<ntruq;j++){
 140               Float_t pjet[4];
 141               pheader->UQJet(j,pjet);
 142               fEventParticles->AddUQJet(pjet);
 143             }
 144           }
 145 #endif
 146           //Int_t ptyp=pheader->ProcessType();
 147           Int_t ntrials=pheader->Trials();
 148           headdesc+=": Tr ";
 149           headdesc+=ntrials;
 150           Int_t ntr=pheader->NTriggerJets();
 151           if(ntr){
 152             for(Int_t j=0;j<ntr;j++){
 153               Float_t pjet[4];
 154               pheader->TriggerJet(j,pjet);
 155               fEventParticles->AddJet(pjet);
 156               if(!ntruq) fEventParticles->AddUQJet(pjet);
 157             }
 158           }
 159         }
 160       }
 161       fEventParticles->SetHeader(headdesc);
 162
 163       //get vertex
 164       const TParticle *kv = stack->Particle(0);
 165       if(kv) {
 166         fEventParticles->SetVertex(kv->Vx(),kv->Vy(),kv->Vz());
 167       }
 168
 169       //loop over particles
 170       for (Int_t i = 0;i<npart; i++)
 171         {
 172           TParticle *p = stack->Particle(i);
 173           if(!p) continue;
 174           Int_t child1 = p->GetFirstDaughter();
 175           //Int_t child2 = p->GetLastDaughter();
 176           //Int_t mother = p->GetFirstMother();
 177           //cout << child1 << " " << child2 << " " << mother << endl;
 178           if((child1>=0) && (child1<nprim)) continue;
 179           if(p->GetStatusCode()!=1){
 180             //p->Print();
 181             continue;
 182           }
 183           if(IsAcceptedParticle(p)) //put particle in event
 184             fEventParticles->AddParticle(p,i);
 185         }
 186       fCurrentEvent++;
 187       fNEventsRead++;
 188       return kTRUE;
 189     }
 190
 191   //end of loop over directories specified in fDirs ObjArray
 192   return kFALSE;
 193 }
 194
 195 Int_t AliJetParticlesReaderKine::OpenFile(Int_t n)
 196 {
 197   //opens file with kine tree
 198
 199   if(fRunLoader){
 200     if(fCurrentEvent < fRunLoader->GetNumberOfEvents()) return kTRUE;
 201     else return kFALSE;
 202   }
 203
 204   const TString& dirname = GetDirName(n);
 205   if (dirname == "")
 206     {
 207       Error("OpenNextFile","Can't get directory name with index %d",n);
 208       return kFALSE;
 209     }
 210
 211   TString filename = dirname +"/"+ fFileName;
 212   fRunLoader = AliRunLoader::Open(filename.Data());
 213
 214   if (fRunLoader == 0)
 215     {
 216       Error("OpenNextFile","Can't open session from file %s",filename.Data());
 217       return kFALSE;
 218     }
 219
 220   if (fRunLoader->GetNumberOfEvents() <= 0)
 221     {
 222       Error("OpenNextFile","There is no events in this directory.");
 223       delete fRunLoader;
 224       fRunLoader = 0;
 225       return kFALSE;
 226     }
 227
 228   if (fRunLoader->LoadKinematics())
 229     {
 230       Error("OpenNextFile","Error occured while loading kinematics.");
 231       return kFALSE;
 232     }
 233
 234   fCurrentEvent = 0;
 235   return kTRUE;
 236 }