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