JETAN/findJets.C

   1 // $Id$
   2
   3 #ifndef __CINT__
   4 #include <stdio.h>
   5 #include <iostream.h>
   6 #include <time.h>
   7 #include <TROOT.h>
   8 #include <TCanvas.h>
   9 #include <TRandom.h>
  10 #include <TSystem.h>
  11 #include <TFile.h>
  12 #include <TChain.h>
  13 #include <TStopwatch.h>
  14 #include <TObjArray.h>
  15 #include <TObjString.h>
  16 #include <TString.h>
  17 #include <AliStack.h>
  18 #include <AliRunLoader.h>
  19 #include <AliHeader.h>
  20 #include <AliGenPythiaEventHeader.h>
  21 #include <AliJetParticle.h>
  22 #include <AliJetParticlesReader.h>
  23 #include <AliJetParticlesReaderKine.h>
  24 #include <AliJetParticlesReaderESD.h>
  25 #include <AliJetParticlesReaderHLT.h>
  26 #include <AliJetEventParticles.h>
  27 #include <AliTkConeJetEvent.h>
  28 #include <AliTkConeJetFinderV2.h>
  29 #include <AliTkChargedJetFinder.h>
  30 #endif
  31
  32 void findJets(Char_t *filename,Float_t ptcut=1.0,Float_t radius=0.3,
  33               Int_t seedpt=-1,Char_t *savename=0,Int_t nMaxEvents=-1)
  34 {
  35   //connect to jets
  36   TChain *theTree = new TChain("AJEPtree");
  37   theTree->Add(filename);
  38   AliJetEventParticles *ev=new AliJetEventParticles();
  39   theTree->SetBranchAddress("particles",&ev);
  40
  41   Int_t treeentries=(Int_t)theTree->GetEntries();
  42   if((nMaxEvents<0) || (nMaxEvents>treeentries))
  43     nMaxEvents=treeentries;
  44   //cout << "Found " << nMaxEvents << " in " << filename << endl;
  45
  46   // create the jet finder
  47   Char_t buffer[8096];
  48   AliTkConeJetFinderV2 *ConeFinder = new AliTkConeJetFinderV2();
  49   ConeFinder->defaultSettings();
  50   //ConeFinder->setSettings(120,40);
  51   if(gErrorIgnoreLevel<=kWarning)
  52     ConeFinder->setOutput(kTRUE);
  53   if(seedpt<5)
  54     ConeFinder->setEtMinJet(5.);
  55   else
  56     ConeFinder->setEtMinJet(seedpt);
  57   ConeFinder->setPtCut(ptcut);
  58   ConeFinder->setEtCut(seedpt);
  59   ConeFinder->setRadius(radius);
  60   if(savename) sprintf(buffer,"%s",savename);
  61   else {
  62     Char_t buffer2[8096];
  63     strncpy(buffer2,filename,strlen(filename)-5);
  64     buffer2[strlen(filename)-5]='\0';
  65     sprintf(buffer,"%s-thresh-%.1f-rad-%.1f.cone.evout.root",buffer2,ptcut,radius);
  66   }
  67   ConeFinder->setEvOutFilename(buffer);
  68   ConeFinder->init();
  69
  70   //=========================================================================
  71   // start the event loop
  72   //=========================================================================
  73   TStopwatch *stopwatch=new TStopwatch();
  74   Int_t nEvent = 0;
  75   while(nEvent<nMaxEvents){
  76
  77     //get the event
  78     theTree->GetEvent(nEvent);
  79
  80     //const AliJetEventParticles *ev=reader->GetEventParticles();
  81     if(gErrorIgnoreLevel<=kWarning){
  82       cout << "Read event: " << nEvent << endl;
  83       ev->Print();
  84     }
  85     TString dummy="Counter: ";
  86     dummy+=nEvent;
  87
  88     stopwatch->Reset();
  89     stopwatch->Start();
  90     ConeFinder->initEvent(ev,dummy);
  91     ConeFinder->run();
  92     ConeFinder->finishEvent();
  93     if(gErrorIgnoreLevel<=kWarning){
  94       cout << "JetFinding done:  CPU time " << stopwatch->CpuTime()
  95            << " Real Time " << stopwatch->RealTime() << endl;
  96     }
  97     nEvent++;
  98     ev->Reset();
  99   } //end of nev loop
 100
 101   ConeFinder->finish();
 102
 103   delete stopwatch;
 104   delete ev;
 105   delete theTree;
 106   delete ConeFinder;
 107 }
 108
 109 void findChargedJets(Char_t *filename,Float_t ptcut=1.0,Float_t radius=0.3,
 110                      Int_t seedpt=10,Char_t *savename=0,Int_t nMaxEvents=-1)
 111 {
 112   //connect to jets
 113   TChain *theTree = new TChain("AJEPtree");
 114   theTree->Add(filename);
 115   AliJetEventParticles *ev=new AliJetEventParticles();
 116   theTree->SetBranchAddress("particles",&ev);
 117
 118   Int_t treeentries=(Int_t)theTree->GetEntries();
 119   if((nMaxEvents<0) || (nMaxEvents>treeentries))
 120     nMaxEvents=treeentries;
 121   //cout << "Found " << nMaxEvents << " in " << filename << endl;
 122
 123   // create the jet finder
 124   Char_t buffer[8096];
 125   AliTkChargedJetFinder *ChFinder = new AliTkChargedJetFinder();
 126   ChFinder->defaultSettings();
 127   //  if(gErrorIgnoreLevel<=kWarning)
 128     ChFinder->setOutput(kTRUE);
 129   ChFinder->setEtMinJet(seedpt);
 130   ChFinder->setPtCut(ptcut);
 131   ChFinder->setPtSeed(seedpt);
 132   ChFinder->setRadius(radius);
 133   if(savename) sprintf(buffer,"%s",savename);
 134   else {
 135     Char_t buffer2[8096];
 136     strncpy(buffer2,filename,strlen(filename)-5);
 137     buffer2[strlen(filename)-5]='\0';
 138     sprintf(buffer,"%s-thresh-%.1f-rad-%.1f.charged.evout.root",buffer2,ptcut,radius);
 139   }
 140   ChFinder->setEvOutFilename(buffer);
 141   ChFinder->init();
 142
 143   //=========================================================================
 144   // start the event loop
 145   //=========================================================================
 146   TStopwatch *stopwatch=new TStopwatch();
 147   Int_t nEvent = 0;
 148   while(nEvent<nMaxEvents){
 149
 150     //get the event
 151     theTree->GetEvent(nEvent);
 152
 153     //  const AliJetEventParticles *ev=reader->GetEventParticles();
 154     if(gErrorIgnoreLevel<=kWarning){
 155       cout << "Read event: " << nEvent << endl;
 156       ev->Print();
 157     }
 158     TString dummy="Counter: ";
 159     dummy+=nEvent;
 160
 161     stopwatch->Reset();
 162     stopwatch->Start();
 163     ChFinder->initEvent(ev,dummy);
 164     ChFinder->run();
 165     ChFinder->finishEvent();
 166     stopwatch->Stop();
 167     if(gErrorIgnoreLevel<=kWarning){
 168       cout << "JetFinding done: CPU time:" << stopwatch->CpuTime()
 169            << " Real Time:" << stopwatch->RealTime() << endl;
 170     }
 171     nEvent++;
 172     ev->Reset();
 173   } //end of nev loop
 174
 175   ChFinder->finish();
 176
 177   delete stopwatch;
 178   delete ev;
 179   delete theTree;
 180   delete ChFinder;
 181 }