JETAN/jetan2004/findJetsMixed.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 <TStopwatch.h>
  18 #include <AliStack.h>
  19 #include <AliRunLoader.h>
  20 #include <AliHeader.h>
  21 #include <AliGenPythiaEventHeader.h>
  22 #include <AliJetParticle.h>
  23 #include <AliJetParticlesReader.h>
  24 #include <AliJetParticlesReaderKine.h>
  25 #include <AliJetParticlesReaderESD.h>
  26 #include <AliJetParticlesReaderHLT.h>
  27 #include <AliJetEventParticles.h>
  28 #include <AliTkConeJetEvent.h>
  29 #include <AliTkConeJetFinderV2.h>
  30 #endif
  31
  32 void findJetsMixed(Char_t *fname,Char_t *hname,Float_t ptcut=2.0,Float_t radius=0.3,Char_t *savename=0,Int_t nMaxEvents=-1)
  33 {
  34   //connect to jets
  35   TChain *theTree = new TChain("AJEPtree");
  36   theTree->Add(fname);
  37   AliJetEventParticles *ev=new AliJetEventParticles();
  38   theTree->SetBranchAddress("particles",&ev);
  39
  40   Int_t treeentries=(Int_t)theTree->GetEntries();
  41   if((nMaxEvents<0) || (nMaxEvents>treeentries))
  42     nMaxEvents=treeentries;
  43   //cout << "Found " << nMaxEvents << " in " << fname << endl;
  44
  45
  46   TChain *backtheTree = new TChain("AJEPtree");
  47   backtheTree->Add(hname);
  48   AliJetEventParticles *backev=new AliJetEventParticles();
  49   backtheTree->SetBranchAddress("particles",&backev);
  50
  51   Int_t backtreeentries=(Int_t)backtheTree->GetEntries();
  52   Int_t nPerBackground=nMaxEvents/backtreeentries;
  53   if(nPerBackground==0) nPerBackground=1;
  54
  55   // create the jet finder
  56   Char_t buffer[8096];
  57   AliTkConeJetFinderV2 *ConeFinder = new AliTkConeJetFinderV2();
  58   ConeFinder->defaultSettings();
  59   //ConeFinder->setSettings(120,40);
  60   if(gErrorIgnoreLevel<=kWarning)
  61     ConeFinder->setOutput(kTRUE);
  62   ConeFinder->setEtMinJet(10);
  63   ConeFinder->setPtCut(ptcut);
  64   ConeFinder->setEtCut(ptcut);
  65   ConeFinder->setRadius(radius);
  66   if(savename) sprintf(buffer,"%s",savename);
  67   else {
  68     Char_t buffer2[8096];
  69     strncpy(buffer2,fname,strlen(fname)-5);
  70     buffer2[strlen(fname)-5]='\0';
  71     sprintf(buffer,"%s-mixed-thresh-%.1f-rad-%.1f.cone.evout.root",buffer2,ptcut,radius);
  72   }
  73   ConeFinder->setEvOutFilename(buffer);
  74   ConeFinder->init();
  75
  76   //=========================================================================
  77   // start the event loop
  78   //=========================================================================
  79   TStopwatch *stopwatch=new TStopwatch();
  80   Int_t nEvent = 0;
  81   Int_t nEventHijing = -1;
  82   Int_t nEventHijingCounter = nPerBackground;
  83   while(nEvent<nMaxEvents){
  84
  85     //get the event
  86     theTree->GetEvent(nEvent);
  87
  88     //load background if needed
  89     if(nEventHijingCounter==nPerBackground){
  90       backev->Reset();
  91       nEventHijing++;
  92       backtheTree->GetEvent(nEventHijing);
  93       if(nEventHijing==backtreeentries) nEventHijing=0;
  94       nEventHijingCounter=0;
  95     }
  96
  97     backev->AddSignal(*ev);
  98     TString dummy="Counter: ";
  99     dummy+=nEvent;
 100     dummy+=" ";
 101     //dummy+=ev->GetHeader();
 102     dummy+="(Pythia ";dummy+=nEvent;
 103     dummy+=" ";
 104     //dummy+=backev->GetHeader();
 105     dummy+=", Hijing ";dummy+=nEventHijing;
 106     dummy+=")";
 107     backev->SetHeader(dummy);
 108     if(gErrorIgnoreLevel<=kWarning){
 109       cout << "Read event: " << nEvent << " " << nEventHijing << endl;
 110       backev->Print();
 111     }
 112
 113     stopwatch->Reset();
 114     stopwatch->Start();
 115     ConeFinder->initEvent(backev,dummy);
 116     ConeFinder->run();
 117     ConeFinder->finishEvent();
 118     if(gErrorIgnoreLevel<=kWarning){
 119       cout << "JetFinding done:  CPU time " << stopwatch->CpuTime()
 120            << " Real Time " << stopwatch->RealTime() << endl;
 121     }
 122
 123     nEvent++;
 124     nEventHijingCounter++;
 125     ev->Reset();
 126   } //end of nev loop
 127
 128   ConeFinder->finish();
 129
 130   delete ev;
 131   delete theTree;
 132   delete backev;
 133   delete backtheTree;
 134   delete ConeFinder;
 135 }
 136