EMCAL/beamtest07/TestESDCaloCluster.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2
   3 //Root include files
   4 #include <Riostream.h>
   5 #include <TFile.h>
   6 #include <TChain.h>
   7 #include <TParticle.h>
   8 #include <TNtuple.h>
   9 #include <TCanvas.h>
  10 #include <TObjArray.h>
  11 #include <TSystem.h>
  12 #include <TString.h>
  13 #include <TH1F.h>
  14 #include <TVector.h>
  15 #include <TParticle.h>
  16
  17 //AliRoot include files
  18 #include "AliESD.h"
  19 #include "AliESDEvent.h"
  20 #include "AliEMCALLoader.h"
  21 #include "AliESDCaloCluster.h"
  22 #include "AliEMCALRecPoint.h"
  23 #include "AliPID.h"
  24 #include "AliLog.h"
  25
  26 //macros
  27
  28 #endif
  29
  30 TChain * AliReadESDfromdisk(const UInt_t eventsToRead,
  31                                    const TString dirName,
  32                                    const TString esdTreeName,
  33                                    const char *  pattern)
  34 {
  35   // Reads ESDs from Disk
  36  TChain *  rv = 0;
  37
  38   // create a TChain of all the files
  39   TChain * cESDTree = new TChain(esdTreeName) ;
  40
  41   // read from the directory file until the require number of events are collected
  42   void * from = gSystem->OpenDirectory(dirName) ;
  43    if (!from)
  44      rv = 0 ;
  45   else{ // reading file names from directory
  46     const char * subdir ;
  47     // search all subdirectories witch matching pattern
  48     while( (subdir = gSystem->GetDirEntry(from))  &&
  49            (cESDTree->GetEntries() < eventsToRead)) {
  50       if ( strstr(subdir, pattern) != 0 ) {
  51         char file[200] ;
  52         sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir);
  53         cESDTree->Add(file) ;
  54       }
  55     } // while file
  56
  57     rv = cESDTree ;
  58
  59   } // reading file names from directory
  60   return rv ;
  61 }
  62
  63 //======================================================================
  64 TChain * AliReadESD(const UInt_t eventsToRead,
  65                   const TString dirName,
  66                   const TString esdTreeName,
  67                   const char *  pattern)
  68 {
  69   // Read AliESDs files and return a Chain of events
  70
  71   if ( dirName == "" )
  72     return 0 ;
  73   if ( esdTreeName == "" )
  74     return AliReadESDfromdisk(eventsToRead, dirName,"","") ;//Last 2 arguments are not necessary but pdsf compiler complains "","'
  75   else if ( strcmp(pattern, "") == 0 )
  76     return AliReadESDfromdisk(eventsToRead, dirName, esdTreeName,"") ;//Last argument is not necessary but pdsf compiler complains "","'
  77   else
  78     return AliReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ;
  79 }
  80
  81 //=====================================================================
  82 //  Do:
  83 //  .L TestESDCaloCluster.C++
  84 //  TestESDCaloCluster(number of events to process)
  85 //=====================================================================
  86 void TestESDCaloCluster(const UInt_t eventsToProcess = 5,
  87                         TString dirName = "./",
  88                         const TString esdTreeName = "esdTree",
  89                         const char *  pattern = ".")
  90 {
  91
  92   //Create chain of esd trees
  93   //AliReadESD(eventsToProcess, directoryName,esdTreeName,patternOfDirectory) ;
  94   //By default the root files are in the same directory
  95   TChain * t = AliReadESD(eventsToProcess, dirName,esdTreeName,pattern) ;
  96
  97
  98   // ESD
  99   AliESDEvent * esd = new AliESDEvent();
 100   esd->ReadFromTree(t);
 101
 102   //Define few variables to be used in macro
 103   TString alirunName = "" ;
 104
 105   //Define example histograms
 106   TH1F * hEnergy = new TH1F("hEnergy","Energy Distribution",100,0.,100.);
 107   TH1F * hEta    = new TH1F("hEta","Eta Distribution",100,-0.7,0.7);
 108   TH1F * hPhi    = new TH1F("hPhi","Phi Distribution",100,0,2*TMath::Pi());
 109
 110   Int_t beg, end ;
 111   UInt_t event ;
 112   Float_t pos[3] ;
 113
 114   for (event = 0; event < eventsToProcess; event++) {//event loop
 115     //AliInfo( Form("Event %d \n",event) );
 116     Int_t nbytes = t->GetEntry(event); // store event in esd
 117     //cout<<"nbytes "<<nbytes<<endl;
 118     if ( nbytes == 0 ) //If nothing in ESD
 119       break ;
 120
 121     // Check that name of file is correct
 122     if (alirunName != t->GetFile()->GetName()) {
 123       alirunName = t->GetFile()->GetName() ;
 124       alirunName.ReplaceAll("galice.root", "AliESDs.root") ;
 125     }
 126
 127     //get reconstructed vertex position
 128
 129     //Double_t vertex_position[3] ;
 130     //    esd->GetVertex()->GetXYZ(vertex_position) ;
 131
 132     cout<<"Event >>>>>>>>>>> "<<event<<endl;
 133
 134
 135     //select EMCAL tracks only
 136     end = esd->GetNumberOfCaloClusters() ;
 137     beg = 0;
 138     Int_t nphP ;
 139     //cout<<"begin "<<beg<<" end "<<end<<endl;
 140     Int_t iclus = 0;
 141     for (nphP =  beg; nphP <  end; nphP++) {//////////////EMCAL cluster loop
 142       AliESDCaloCluster * clus = esd->GetCaloCluster(nphP) ; // retrieve cluster from esd
 143
 144       //Get the cluster parameters
 145       Float_t energy   = clus->E() ;
 146       cout << "Cluster " << nphP << " energy = " << energy << endl;
 147
 148         //      Int_t mult       = clus->GetNumberOfDigits() ;
 149         //      Float_t disp     = clus->GetClusterDisp() ;
 150         //      UShort_t * amp   = clus->GetDigitAmplitude() ;
 151         //      UShort_t * time  = clus->GetDigitTime() ;
 152         //      UShort_t * index = clus->GetDigitIndex() ;
 153         //      Int_t iprim = clus->GetPrimaryIndex();
 154
 155       /*
 156       clus->GetPosition(pos) ;
 157       TVector3 vpos(pos[0],pos[1],pos[2]) ;
 158       //Print values on screen
 159       cout<<"ESD cluster "<<iclus <<"; Energy "<<energy
 160           <<"; Phi "<<vpos.Phi()*180/TMath::Pi()<<"; Eta "<<vpos.Eta()<<endl;
 161       //cout<<"Dispersion "<<disp<<"; multiplicity "<<mult<<endl;
 162       //Get the parent main values
 163       */
 164       //Fill histograms
 165       hEnergy->Fill(energy) ;
 166       //      hEta->Fill(vpos.Eta()) ;
 167       //hPhi->Fill(vpos.Phi()) ;
 168
 169     }// track loop
 170   }//////////////////////////////////////////event loop
 171   hEnergy->Draw();//Draw histogram on screen
 172   //Write histograms in Root file
 173   TFile outf("histo.root","recreate") ;
 174   hEnergy->Write() ;
 175   hPhi->Write() ;
 176   hEta->Write() ;
 177   outf.Close() ;
 178 }
 179
 180