EMCAL/macros/TestESDCaloClusterAndCell.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 #include <TRefArray.h>
  17 #include <TArrayS.h>
  18
  19 //AliRoot include files
  20 #include "AliRunLoader.h"
  21 #include "AliStack.h"
  22 #include "AliESDEvent.h"
  23 #include "AliESDVertex.h"
  24 #include "AliESDCaloCluster.h"
  25 #include "AliESDCaloCells.h"
  26 #include "AliPID.h"
  27 #include "AliLog.h"
  28
  29 #endif
  30
  31 //example macro to extract information from AliESDCaloClusters and AliESDCaloCells.
  32 // Author: Gustavo.Conesa.Balbastre@cern.ch (INFN-LNF)
  33
  34 TChain * AliReadESDfromdisk(const UInt_t eventsToRead,
  35                                    const TString dirName,
  36                                    const TString esdTreeName,
  37                                    const char *  pattern)
  38 {
  39   // Reads ESDs from Disk
  40  TChain *  rv = 0;
  41
  42   // create a TChain of all the files
  43   TChain * cESDTree = new TChain(esdTreeName) ;
  44
  45   // read from the directory file until the require number of events are collected
  46   void * from = gSystem->OpenDirectory(dirName) ;
  47    if (!from)
  48      rv = 0 ;
  49   else{ // reading file names from directory
  50     const char * subdir ;
  51     // search all subdirectories witch matching pattern
  52     while( (subdir = gSystem->GetDirEntry(from))  &&
  53            (cESDTree->GetEntries() < eventsToRead)) {
  54       if ( strstr(subdir, pattern) != 0 ) {
  55         char file[200] ;
  56         sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir);
  57         cESDTree->Add(file) ;
  58       }
  59     } // while file
  60
  61     rv = cESDTree ;
  62
  63   } // reading file names from directory
  64   return rv ;
  65 }
  66
  67 //======================================================================
  68 TChain * AliReadESD(const UInt_t eventsToRead,
  69                   const TString dirName,
  70                   const TString esdTreeName,
  71                   const char *  pattern)
  72 {
  73   // Read AliESDs files and return a Chain of events
  74
  75   if ( dirName == "" )
  76     return 0 ;
  77   if ( esdTreeName == "" )
  78     return AliReadESDfromdisk(eventsToRead, dirName,"","") ;//Last 2 arguments are not necessary but pdsf compiler complains "","'
  79   else if ( strcmp(pattern, "") == 0 )
  80     return AliReadESDfromdisk(eventsToRead, dirName, esdTreeName,"") ;//Last argument is not necessary but pdsf compiler complains "","'
  81   else
  82     return AliReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ;
  83 }
  84
  85 //=====================================================================
  86 //  Do:
  87 //  .L TestESDCaloClusterAndCell.C++
  88 //  TestESDCaloClusterAndCell(calo, number of events to process)
  89 //=====================================================================
  90 void TestESDCaloClusterAndCell(TString det = "EMCAL", const UInt_t eventsToProcess = 5,
  91                         TString dirName = "./",
  92                         const TString esdTreeName = "esdTree",
  93                         const char *  pattern = ".")
  94 {
  95
  96   if(det !="EMCAL" && det !="PHOS" )
  97  {
  98    cout << "Wrong detector name "<<det<<endl;
  99    return;
 100  }
 101
 102   //Create chain of esd trees
 103   //By default the root files are in the same directory
 104   TChain * t = AliReadESD(eventsToProcess, dirName,esdTreeName,pattern) ;
 105
 106
 107   // ESD
 108   AliESDEvent * esd = new AliESDEvent();
 109   esd->ReadFromTree(t); // This checks also if the branch address is already set
 110
 111
 112   //Define few variables to be used in macro
 113   TString alirunName = "" ;
 114
 115   UInt_t event ;
 116   Float_t pos[3] ;
 117
 118   for (event = 0; event < eventsToProcess; event++) {//event loop
 119     //AliInfo( Form("Event %d \n",event) );
 120     Int_t nbytes = t->GetEntry(event); // store event in esd
 121     if ( nbytes == 0 ) //If nothing in ESD
 122       break ;
 123
 124     // Check that name of file is correct
 125     if (alirunName != t->GetFile()->GetName()) {
 126       alirunName = t->GetFile()->GetName() ;
 127       alirunName.ReplaceAll("galice.root", "AliESDs.root") ;
 128     }
 129
 130     AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),  "read");
 131     rl->LoadKinematics();
 132     rl->GetEvent(event);
 133     AliStack *sta=rl->Stack();
 134
 135     //get reconstructed vertex position
 136     Double_t vertex_position[3] ;
 137     esd->GetVertex()->GetXYZ(vertex_position) ;
 138
 139     //cout<<"Event >>>>>>>>>>>  "<<event<<" vertex >>> "<<vertex_position[0]<<" "<<vertex_position[1]<<" "<<vertex_position[2]<<endl;
 140
 141     //Get CaloCells
 142
 143     AliESDCaloCells &cells= *(esd->GetEMCALCells());
 144     if(det == "PHOS")  AliESDCaloCells &cells = *(esd->GetPHOSCells());
 145     Int_t ncell = cells.GetNumberOfCells() ;
 146     Int_t type = cells.GetType();
 147     cout<<">>> Event "<<event<<" ncells "<<ncell<< " type "<<type<<endl;
 148     // Uncomment to see the full list of digit amplitudes and times.
 149 //     for (Int_t icell=  0; icell <  ncell; icell++) {
 150 //       cout<<"icell "<<icell<<
 151 //      " ID "<<cells.GetCellNumber(icell)<<
 152 //      " amp "<<cells.GetAmplitude(icell)<<
 153 //      " time "<<cells.GetTime(icell)<<endl;
 154
 155 //     }// cell loop
 156
 157     //Get the CaloClusters
 158
 159     //select EMCAL clusters only
 160     TRefArray * caloClusters  = new TRefArray();
 161
 162     if(det == "EMCAL") esd->GetEMCALClusters(caloClusters);
 163     else if(det == "PHOS") esd->GetPHOSClusters(caloClusters);
 164
 165     Int_t nclus = caloClusters->GetEntries() ;
 166
 167     cout<<"Event: "<<event<<"; Number of clusters "<<nclus<<endl;
 168     for (Int_t iclus =  0; iclus <  nclus; iclus++) {//////////////EMCAL cluster loop
 169       AliESDCaloCluster * clus = (AliESDCaloCluster *) caloClusters->At(iclus) ; // retrieve cluster from esd
 170       //cout<<">>>> Cluster >>>"<<iclus<<endl;
 171
 172       //Get the cluster info
 173
 174       Float_t energy   = clus->E() ;
 175       //Float_t disp     = clus->GetClusterDisp() ;
 176       Int_t iprim = clus->GetLabel();
 177
 178       clus->GetPosition(pos) ;
 179       TVector3 vpos(pos[0],pos[1],pos[2]) ;
 180       TLorentzVector p ;
 181       clus->GetMomentum(p,vertex_position);
 182
 183       Int_t mult       = clus->GetNCells() ;
 184       UShort_t * index = clus->GetCellsAbsId() ;
 185       Double_t * fraction = clus->GetCellsAmplitudeFraction() ;
 186       //Print cluster info
 187       cout<<"Cluster "<<iclus <<"; digits mult "<<mult<<"; type "<<(Int_t )clus->GetClusterType()
 188           <<"; Energy "<<energy
 189           <<"; Phi "<<vpos.Phi()*180/TMath::Pi()<<"; Eta "<<vpos.Eta()
 190           <<"; label "<<iprim<<endl;
 191       //Print primary info
 192       if(iprim>=0 && iprim < sta->GetNtrack()){
 193         TParticle * particle = sta->Particle(clus->GetLabel());
 194         //Print primary values
 195         cout<<" Primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl;
 196       }
 197       else if( iprim >= sta->GetNtrack()) cout <<"PROBLEM, label is too large : "<<iprim<<" >= particles in stack "<< sta->GetNtrack() <<endl;
 198       else cout<<"Negative label!!!  : "<<iprim<<endl;
 199
 200       //Get CaloCells of cluster and print
 201       for(Int_t i = 0; i < mult ; i++){
 202         Int_t absId =   index[i]; // or clus->GetCellNumber(i) ;
 203         Double_t ampFract =  fraction[i];
 204         Float_t amp = cells.GetCellAmplitude(absId) ;
 205         Float_t time = cells.GetCellTime(absId);
 206         cout<<"AbsID : "<< absId << "; Amplitude "<< amp << "; Fraction "<<ampFract<<"; Time " <<time<<endl;
 207       }
 208
 209     }// clusters
 210
 211   }// event loop
 212
 213 }
 214
 215