JETAN/jetan2004/testJets.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 <AliTkConeJetEvent.h>
  22 #include <AliTkConeJetFinderV2.h>
  23 #include <AliJetParticle.h>
  24 #include <AliJetParticlesReader.h>
  25 #include <AliJetParticlesReaderKine.h>
  26 #include <AliJetParticlesReaderESD.h>
  27 #include <AliJetParticlesReaderHLT.h>
  28 #include <AliJetEventParticles.h>
  29 #endif
  30
  31 #if 0
  32     kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
  33     kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCpid=0x0080,
  34     kTRDin=0x0100,kTRDout=0x0200,kTRDrefit=0x0400,kTRDpid=0x0800,
  35     kTOFin=0x1000,kTOFout=0x2000,kTOFrefit=0x4000,kTOFpid=0x8000,
  36     kPHOSpid=0x10000, kRICHpid=0x20000,
  37     kTRDStop=0x20000000,
  38     kESDpid=0x40000000,
  39     kTIME=0x80000000
  40 #endif
  41
  42 void testJets(Char_t *input,Int_t esd=0,Float_t ptcut=2.0,Float_t radius=0.3)
  43 {
  44
  45   TObjArray *dirs=new TObjArray(10);
  46   dirs->Add(new TObjString(input));
  47
  48   AliJetParticlesReader *reader=0;
  49   TString desc="";
  50   if (esd==0) {
  51     reader=new AliJetParticlesReaderKine(dirs);
  52     //((AliJetParticlesReaderKine*)reader)->SetUseTracks(kTRUE);
  53     ((AliJetParticlesReaderKine*)reader)->SetCharged(kTRUE);
  54     ((AliJetParticlesReaderKine*)reader)->SetEM(kFALSE);
  55     ((AliJetParticlesReaderKine*)reader)->SetNeutral(kFALSE);
  56     desc+="Kine";
  57   } else if(esd==1){
  58     reader=new AliJetParticlesReaderESD(0,dirs);
  59     ((AliJetParticlesReaderESD*)reader)->SetCompareFlag(0x0030);
  60     desc+="ESD";
  61   } else if (esd==2) {
  62     reader=new AliJetParticlesReaderHLT(kTRUE,dirs);
  63     desc+="HLTConf";
  64   } else if (esd==3) {
  65     reader=new AliJetParticlesReaderHLT(kFALSE,dirs);
  66     desc+="HLTHough";
  67    } else if (esd==10) {
  68     reader=new AliJetParticlesReaderKine(dirs);
  69     //((AliJetParticlesReaderKine*)reader)->SetUseTracks(kTRUE);
  70     ((AliJetParticlesReaderKine*)reader)->SetCharged(kTRUE);
  71     ((AliJetParticlesReaderKine*)reader)->SetEM(kTRUE);
  72     ((AliJetParticlesReaderKine*)reader)->SetNeutral(kTRUE);
  73     desc+="Kine-All";
  74    } else if (esd==11) {
  75     reader=new AliJetParticlesReaderKine(dirs);
  76     //((AliJetParticlesReaderKine*)reader)->SetUseTracks(kTRUE);
  77     ((AliJetParticlesReaderKine*)reader)->SetCharged(kTRUE);
  78     ((AliJetParticlesReaderKine*)reader)->SetEM(kTRUE);
  79     ((AliJetParticlesReaderKine*)reader)->SetNeutral(kFALSE);
  80     desc+="Kine-EM";
  81   } else {
  82     cout << "Parameter settings: " << endl;
  83     cout << "0  == Kine (ch)"  << endl;
  84     cout << "1  == TPC"  << endl;
  85     cout << "2  == Conformal"  << endl;
  86     cout << "3  == Hough"  << endl;
  87     cout << "10 == Kine (all)"  << endl;
  88     cout << "11 == Kine (ch+em)"  << endl;
  89     return;
  90   }
  91   //reader->ReadEventsFromTo(0,10);
  92   reader->SetPtCut(ptcut,100);
  93   reader->SetEtaCut(-0.9,0.9);
  94
  95   desc+=" Ptcut: ";
  96   desc+=ptcut;
  97
  98   // create the jet finder
  99   Char_t buffer[8096];
 100   AliTkConeJetFinderV2 *ConeFinder = new AliTkConeJetFinderV2();
 101   ConeFinder->defaultSettings();
 102   //ConeFinder->setSettings(120,40);
 103   if(gErrorIgnoreLevel<=kWarning)
 104     ConeFinder->setOutput(kTRUE);
 105   ConeFinder->setEtMinJet(1.);
 106   ConeFinder->setEtCut(-1);
 107   ConeFinder->setRadius(radius);
 108   sprintf(buffer,"./conefinder-type%d.evout.root",esd);
 109   //cout << "Cone Jet Finder init: " << endl;
 110   ConeFinder->setEvOutFilename(buffer);
 111   ConeFinder->init();
 112
 113   Int_t counter=0;
 114   while(reader->Next())
 115     {
 116       const AliJetEventParticles *ev=reader->GetEventParticles();
 117       if(gErrorIgnoreLevel<=kWarning){
 118         cout << "Read event: " << counter << endl;
 119         ev->Print();
 120       }
 121       TString dummy="Counter: ";
 122       dummy+=counter;
 123       dummy+=" ";
 124       dummy+=desc;
 125       ConeFinder->initEvent(ev,dummy);
 126       ConeFinder->run();
 127       ConeFinder->finishEvent();
 128     }
 129   ConeFinder->finish();
 130
 131   delete ConeFinder;
 132   delete reader;
 133 }
 134
 135 void readJets(Int_t nMaxEvents,Char_t *filename)
 136 {
 137   //connect to jets
 138   TChain *theTree = new TChain("jets");
 139   theTree->Add(filename);
 140   AliTkConeJetEvent *event = new AliTkConeJetEvent();
 141   theTree->SetBranchAddress("ConeFinder",&event);
 142
 143   Int_t treeentries=(Int_t)theTree->GetEntries();
 144   if((nMaxEvents<0) || (nMaxEvents>treeentries))
 145     nMaxEvents=treeentries;
 146
 147   cout << "Found " << nMaxEvents << " in " << filename << endl;
 148
 149   //=========================================================================
 150   // start the event loop
 151   //=========================================================================
 152   Int_t nEvent = 0;
 153   while(nEvent<nMaxEvents){
 154     if ((nEvent % 100) == 0) {
 155       cout << "Analysing event " << nEvent << endl;
 156     }
 157
 158     //connect the cone jets
 159     theTree->GetEvent(nEvent);
 160
 161     AliJetEventParticles *jetparts=event->getJetParticles();
 162     jetparts->Print();
 163 #if 0
 164     const TClonesArray *p=jetparts->GetParticles();
 165     for(Int_t i=0;i<p->GetEntriesFast();i++) {
 166       cout << i << endl;
 167       ((AliJetParticle*)p->At(i))->Print();
 168     }
 169 #endif
 170
 171     TClonesArray *tkjets=event->getJets();
 172     if(!tkjets){
 173       cerr << "No Cone jet found in event " << nEvent << ", exiting..." <<endl;
 174       continue;
 175     }
 176
 177     event->Print("");
 178     nEvent++;
 179     event->Clear();
 180   } //end of nev loop
 181
 182   delete event;
 183   delete theTree;
 184 }