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 <TLorentzVector.h>
  18 #include <AliRunLoader.h>
  19 #include <AliStack.h>
  20 #include <AliHeader.h>
  21 #include <AliGenEventHeader.h>
  22 #include <AliGenPythiaEventHeader.h>
  23 #include <AliGenHijingEventHeader.h>
  24 #include "AliJetParticle.h"
  25 #include "AliJetEventParticles.h"
  26 #include "AliJetParticlesReaderKine.h"
  27
  28 ClassImp(AliJetParticlesReaderKine)
  29
  30 AliJetParticlesReaderKine::AliJetParticlesReaderKine() :
  31   AliJetParticlesReader(),
  32   fFileName("galice.root"),
  33   fRunLoader(0),
  34   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
  35   fUseTracks(kFALSE)
  36 {
  37   //constructor
  38 }
  39
  40 AliJetParticlesReaderKine::AliJetParticlesReaderKine(TString& fname) :
  41   AliJetParticlesReader(),
  42   fFileName(fname),
  43   fRunLoader(0),
  44   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
  45   fUseTracks(kFALSE)
  46 {
  47   //constructor
  48 }
  49
  50 AliJetParticlesReaderKine::AliJetParticlesReaderKine(TObjArray* dirs, const Char_t *filename):
  51   AliJetParticlesReader(dirs),
  52   fFileName(filename),
  53   fRunLoader(0),
  54   fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
  55   fUseTracks(kFALSE)
  56 {
  57   //constructor
  58 }
  59
  60 AliJetParticlesReaderKine::~AliJetParticlesReaderKine()
  61 {
  62   //destructor
  63   if(fRunLoader) delete fRunLoader;
  64 }
  65
  66 void AliJetParticlesReaderKine::Rewind()
  67 {
  68   //Rewinds to the beginning
  69   if(fRunLoader) delete fRunLoader;
  70   fRunLoader   = 0;
  71   fCurrentDir  = 0;
  72   fNEventsRead = 0;
  73 }
  74
  75 Int_t AliJetParticlesReaderKine::ReadNext()
  76 {
  77   //Reads Kinematics Tree
  78
  79   if((!fOwner) || (fEventParticles == 0))
  80     fEventParticles = new AliJetEventParticles();
  81
  82   while(fCurrentDir < GetNumberOfDirs())
  83     {
  84
  85       if (!OpenFile(fCurrentDir))
  86       {
  87         delete fRunLoader; //close current session
  88         fRunLoader = 0;    //assure pointer is null
  89         fCurrentDir++;
  90         continue;
  91       }
  92
  93       if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
  94         {
  95           //read next directory
  96           delete fRunLoader; //close current session
  97           fRunLoader = 0;    //assure pointer is null
  98           fCurrentDir++;     //go to next dir
  99           continue;
 100         }
 101
 102       Info("ReadNext","Reading Event %d",fCurrentEvent);
 103       fRunLoader->GetEvent(fCurrentEvent);
 104       AliStack* stack = fRunLoader->Stack();
 105       if (!stack)
 106         {
 107           Error("ReadNext","Can't get stack for event %d",fCurrentEvent);
 108           continue;
 109         }
 110
 111       //clear old event
 112       Int_t nprim = stack->GetNprimary();
 113       Int_t npart = nprim;
 114       if(fUseTracks)
 115         npart = stack->GetNtrack();
 116       fEventParticles->Reset(npart);
 117
 118       TString headdesc="";
 119       Int_t gentype=-1;
 120       AliHeader *header=fRunLoader->GetHeader();
 121       if(!header) {
 122         Warning("ReadNext","Header not found in event %d",fCurrentEvent);
 123       } else {
 124         AliGenEventHeader *genheader=header->GenEventHeader();
 125         TString hname=genheader->GetName();
 126         if(hname.CompareTo("HIJINGparam")==0) {
 127           gentype=1;      ;
 128           headdesc+="HIJINGparam";
 129         } else if(hname.CompareTo("Hijing")==0) {
 130           gentype=2;
 131           AliGenHijingEventHeader *hheader=(AliGenHijingEventHeader*)header->GenEventHeader();
 132           headdesc+="Hijing";
 133           if(!hheader) {
 134             Warning("ReadNext","Hijing-Header not found in event %d",fCurrentEvent);
 135           } else {
 136             Int_t ntrials=hheader->Trials();
 137             fEventParticles->SetTrials(ntrials);
 138             fEventParticles->SetImpact(hheader->ImpactParameter());
 139             fEventParticles->SetNhard(hheader->HardScatters());
 140             fEventParticles->SetNpart(hheader->NwNw());
 141
 142             TLorentzVector  jet1;
 143             TLorentzVector  jet2;
 144             TLorentzVector  jetFsr1;
 145             TLorentzVector  jetFsr2;
 146             hheader->GetJets(jet1,jet2,jetFsr1,jetFsr2);
 147             fEventParticles->AddHard(0,jet1.X(),jet1.Y(),jet1.Z(),jet1.T(),0);
 148             fEventParticles->AddHard(1,jet2.X(),jet2.Y(),jet2.Z(),jet2.T(),0);
 149           }
 150         } else if(hname.CompareTo("Pythia")==0){
 151           gentype=3;
 152           headdesc+="Pythia";
 153           AliGenPythiaEventHeader *pheader=(AliGenPythiaEventHeader*)header->GenEventHeader();
 154           if(!pheader) {
 155             Warning("ReadNext","Pythia-Header not found in event %d",fCurrentEvent);
 156           } else {
 157             Int_t ntruq=0;
 158 #ifndef NOUQHEADERINFO
 159             ntruq=pheader->NUQTriggerJets();
 160             if(ntruq){
 161               Double_t x0=pheader->GetXJet();
 162               Double_t y0=pheader->GetYJet();
 163               if(x0==y0==-1){
 164                 x0=y0=0.;
 165               }
 166               Double_t zquench[4];
 167               pheader->GetZQuench(zquench);
 168               fEventParticles->SetXYJet(x0,y0);
 169               fEventParticles->SetZQuench(zquench);
 170               for(Int_t j=0;j<ntruq;j++){
 171                 Float_t pjet[4];
 172                 pheader->UQJet(j,pjet);
 173                 fEventParticles->AddUQJet(pjet);
 174               }
 175             }
 176 #endif
 177             //Int_t ptyp=pheader->ProcessType();
 178             Int_t ntrials=pheader->Trials();
 179             fEventParticles->SetTrials(ntrials);
 180
 181             Int_t ntr=pheader->NTriggerJets();
 182             if(ntr){
 183               for(Int_t j=0;j<ntr;j++){
 184                 Float_t pjet[4];
 185                 pheader->TriggerJet(j,pjet);
 186                 fEventParticles->AddJet(pjet);
 187                 if(!ntruq) fEventParticles->AddUQJet(pjet);
 188               }
 189             }
 190             for(Int_t i=6;i<=7;i++){
 191               TParticle *MP = stack->Particle(i);
 192               if(!MP) break;
 193               Int_t type=0;
 194               if(MP->GetPdgCode()==21) type=2;
 195               else type=1;
 196               fEventParticles->AddHard(i-6,MP->Px(),MP->Py(),MP->Pz(),MP->Energy(),type);
 197             }
 198           }
 199         }
 200         headdesc+="Run ";
 201         headdesc+=header->GetRun();
 202         headdesc+=": Ev ";
 203         headdesc+=header->GetEventNrInRun();
 204       }
 205       fEventParticles->SetHeader(headdesc);
 206
 207       //get vertex
 208       const TParticle *kv = stack->Particle(0);
 209       if(kv) {
 210         fEventParticles->SetVertex(kv->Vx(),kv->Vy(),kv->Vz());
 211       }
 212
 213       //loop over particles
 214       for (Int_t i = 0;i<npart; i++)
 215         {
 216           TParticle *p = stack->Particle(i);
 217           if(!p) continue;
 218           Int_t child1 = p->GetFirstDaughter();
 219           //Int_t child2 = p->GetLastDaughter();
 220           //Int_t mother = p->GetFirstMother();
 221           //cout << child1 << " " << child2 << " " << mother << endl;
 222           if((child1>=0) && (child1<nprim)) continue;
 223
 224           //check status code depending on gentype
 225           if(gentype==1){ // Hijing Param
 226             if(p->GetStatusCode()!=0) continue;
 227           } else if(gentype==2){ //Hijing
 228             //if( p->GetStatusCode()!=0){
 229             //  cout <<p->GetStatusCode() << " ";p->Print();
 230             //}
 231           } else if(gentype==3){ //Pythia
 232             if(p->GetStatusCode()!=1) continue;
 233           }
 234
 235           //kinematic cuts
 236           if(IsAcceptedParticle(p)){ //put particle in event
 237             //if(p->Pt()>20){cout <<p->GetStatusCode() << " ";p->Print(); }
 238             fEventParticles->AddParticle(p,i);
 239           }
 240         }
 241       fCurrentEvent++;
 242       fNEventsRead++;
 243       return kTRUE;
 244     }
 245
 246   //end of loop over directories specified in fDirs ObjArray
 247   return kFALSE;
 248 }
 249
 250 Int_t AliJetParticlesReaderKine::OpenFile(Int_t n)
 251 {
 252   //opens file with kine tree
 253
 254   if(fRunLoader){
 255     if(fCurrentEvent < fRunLoader->GetNumberOfEvents()) return kTRUE;
 256     else return kFALSE;
 257   }
 258
 259   const TString& dirname = GetDirName(n);
 260   if (dirname == "")
 261     {
 262       Error("OpenNextFile","Can't get directory name with index %d",n);
 263       return kFALSE;
 264     }
 265
 266   TString filename = dirname +"/"+ fFileName;
 267   fRunLoader = AliRunLoader::Open(filename.Data());
 268
 269   if (fRunLoader == 0)
 270     {
 271       Error("OpenNextFile","Can't open session from file %s",filename.Data());
 272       return kFALSE;
 273     }
 274
 275   if (fRunLoader->GetNumberOfEvents() <= 0)
 276     {
 277       Error("OpenNextFile","There is no events in this directory.");
 278       delete fRunLoader;
 279       fRunLoader = 0;
 280       return kFALSE;
 281     }
 282
 283   if (fRunLoader->LoadKinematics())
 284     {
 285       Error("OpenNextFile","Error occured while loading kinematics.");
 286       return kFALSE;
 287     }
 288
 289   fCurrentEvent = 0;
 290   return kTRUE;
 291 }