PHOS/AnaESD.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id$ */
  16 //_________________________________________________________________________
  17 // Macros analyzing the ESD file
  18 // Use Case :
  19 //          root> .L AnaESD.C++
  20 //          root> ana() --> prints the objects stored in ESD
  21 //
  22 // author  : Yves Schutz (CERN/SUBATECH)
  23 // February 2004
  24 //_________________________________________________________________________
  25 #if !defined(__CINT__) || defined(__MAKECINT__)
  26 #include "TFile.h"
  27 #include "TMath.h"
  28 #include "TH1D.h"
  29 #include "TROOT.h"
  30 #include "TBrowser.h"
  31 #include "AliPHOSGetter.h"
  32 #include "AliPHOSGeometry.h"
  33 #include "Riostream.h"
  34 #include "AliESD.h"
  35 #include "AliESDtrack.h"
  36 #include "AliESDCaloTrack.h"
  37 #include "AliEMCALRecParticle.h"
  38 #include "AliPHOSRecParticle.h"
  39 #include "AliKalmanTrack.h"
  40 #include "AliPHOSGridFile.h"
  41 #endif
  42
  43 Bool_t AnaESD(TString filename) ;
  44 void   Match(AliESDtrack * ct, AliESDtrack * cp, Double_t * dist) ;
  45
  46
  47 //__________________________________________________________________________
  48 Bool_t Ana(const TString type = "per5", const Int_t run = 1, const Int_t nOfEvt = 1)
  49 {
  50   // Analyzes data from the AliEn data catalog
  51
  52   // Data Challenge identification
  53   const TString kYear("2004") ;
  54   const TString kProd("02") ;
  55   const TString kVers("V4.01.Rev.00") ;
  56
  57   // get the LFN file name in the Grid catalogue ;
  58   AliPHOSGridFile lfn ;
  59   if (!lfn.IsConnected())
  60     return kFALSE ;
  61   lfn.SetPath(kYear, kProd, kVers, type) ;
  62   lfn.SetRun(run) ;
  63
  64   //loop over the events
  65   Int_t nevt, evt = 0 ;
  66   for (nevt = 0 ; nevt < nOfEvt ; nevt++) {
  67     evt++ ;
  68     lfn.SetEvt(evt) ;
  69     TString fileName = lfn.GetLFN() ;
  70
  71     if (fileName.IsNull()) {
  72       nevt-- ;
  73       continue ;
  74     }
  75
  76     printf(">>>>>>>>>>>> Processing %s-%s/%s/%s : run # %d event # %d \n",
  77            kYear.Data(), kProd.Data(), kVers.Data(), type.Data(), run, evt) ;
  78     AnaESD(fileName) ;
  79   }
  80   return kTRUE ;
  81 }
  82
  83 //__________________________________________________________________________
  84 Bool_t AnaESD(TString fileName)
  85 {
  86   // Analyze ESD from file fileName
  87   // calo track histograms ;
  88
  89   TFile * out = new TFile("AOD.root", "RECREATE") ;
  90   TH1D * hCaloEnergyA = 0 ;
  91   TH1D * hCaloEnergyE = 0 ;
  92   TH1D * hCaloEnergyG = 0 ;
  93   TH1D * hCaloEnergyP = 0 ;
  94   TH1D * heta = 0 ;
  95   TH1D * hphi = 0 ;
  96   Double_t dist[3] ;
  97
  98   AliPHOSGetter * gime = AliPHOSGetter::Instance(fileName) ;
  99
 100   Int_t nEvent = gime->MaxEvent() ;
 101   Int_t event ;
 102   AliESD * esd = 0 ;
 103   for (event = 0 ; event < nEvent; event++) {
 104     cout << "AnaESD : Processing event # " << event << endl ;
 105     esd = gime->ESD(event) ;
 106     if (!esd)
 107       return kFALSE ;
 108
 109     //esd->Print();
 110     // Calorimeter tracks
 111     AliESDtrack * ct ;
 112     Int_t caloindex ;
 113     for (caloindex = 0 ; caloindex < esd->GetNumberOfTracks() ; caloindex++) {
 114       // get the tracks and check if it is from PHOS
 115       ct = esd->GetTrack(caloindex) ;
 116       if ( !ct->IsPHOS() )
 117         continue ;
 118       Double_t energy = ct->GetPHOSsignal() ;
 119       // check if CPV bit is set (see AliPHOSFastRecParticle)
 120       Double_t type[AliESDtrack::kSPECIES+4] ;
 121       ct->GetPHOSpid(type) ;
 122
 123       if ( (type[AliESDtrack::kElectron] == 1.0) || (type[AliESDtrack::kPhoton] == 1.0) ) {
 124         if(!hCaloEnergyA) {
 125           out->cd() ;
 126           hCaloEnergyA = new TH1D("CaloEnergyA", "Energy in calorimeter Electron/Photon like", 500, 0., 50.) ;
 127         }
 128         hCaloEnergyA->Fill(energy) ;
 129       }
 130
 131       if ( type[AliESDtrack::kElectron] == 1.0 ) {
 132         if(!hCaloEnergyE) {
 133           out->cd() ;
 134           hCaloEnergyE = new TH1D("CaloEnergyE", "Energy in calorimeter Electron like", 500, 0., 50.) ;
 135         }
 136         hCaloEnergyE->Fill(energy) ;
 137       }
 138
 139       if ( type[AliESDtrack::kPhoton] == 1.0 ) {
 140         if(!hCaloEnergyG) {
 141           out->cd() ;
 142           hCaloEnergyG = new TH1D("CaloEnergyG", "Energy in calorimeter Gamma like", 500, 0., 50.) ;
 143         }
 144         hCaloEnergyG->Fill(energy) ;
 145       }
 146
 147       if ( type[AliESDtrack::kPi0] == 1.0 ) {
 148         if(!hCaloEnergyP) {
 149           out->cd() ;
 150           hCaloEnergyP = new TH1D("CaloEnergyP", "Energy in calorimeter Pi0 like", 500, 0., 100.) ;
 151         }
 152         hCaloEnergyP->Fill(energy) ;
 153       }
 154
 155       AliESDtrack * cp ;
 156       Int_t cpindex ;
 157       for (cpindex = 0 ; cpindex < esd->GetNumberOfTracks() ; cpindex++) {
 158         // get the charged tracks from central tracking
 159         cp = esd->GetTrack(cpindex) ;
 160         if ( cp->IsPHOS() )
 161           continue ;
 162         Match(ct, cp, dist) ;
 163         if (!heta && !hphi) {
 164           heta = new TH1D("Correta", "neutral-charged correlation in eta" , 100, 0., 360.) ;
 165           hphi = new TH1D("Corrphi", "neutral-charged correlation in phi" , 100, 0., 360.) ;
 166         }
 167         heta->Fill(dist[1]) ;
 168         hphi->Fill(dist[2]) ;
 169       }
 170     }
 171   }
 172   TBrowser * bs = new TBrowser("Root Memory Bowser", gROOT->FindObjectAny("AOD.root") ) ;
 173   bs->Show() ;
 174   out->Write() ;
 175   return kTRUE ;
 176 }
 177
 178 //__________________________________________________________________________
 179 void Match(AliESDtrack * ct, AliESDtrack * cp, Double_t * dist)
 180 {
 181   // Calculates the distance (x,z) between  the particle detected by PHOS and
 182   // the charged particle reconstructed by the global tracking
 183
 184
 185  //  Int_t phN ;
 186 //   Double_t phZ, phX ;
 187
 188 //   if (ct->IsPHOS()) { // it is a PHOS particle
 189 //     Double_t cpTheta,  cpPhi ;
 190 //     Double_t phTheta,  phPhi ;
 191 //     cpTheta = cpPhi = phTheta = phPhi = 0. ;
 192 //    //    cout << "PHOS particle # " << " pos ("
 193 //     //        << pp->GetPos().X() << ", " << pp->GetPos().Y() << ", " << pp->GetPos().Z() << ")" << endl ;
 194
 195 //     AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 196 //     gime->PHOSGeometry()->ImpactOnEmc(*pp, phN, phZ, phX) ;
 197 //     Double_t xyzAtPHOS[3] ;
 198 //     cp->GetOuterXYZ(xyzAtPHOS) ;
 199 //     if ( (xyzAtPHOS[0] +  xyzAtPHOS[1] + xyzAtPHOS[2]) != 0.) { //it has reached PHOS
 200 //       //the next check are only if we want high quality tracks
 201 //       //       ULong_t status = cp->GetStatus() ;
 202 //       //       if ((status & AliESDtrack::kTRDput)==0)
 203 //       //     do not continue;
 204 //       //       if ((status & AliESDtrack::kTRDStop)!=0)
 205 //       //     do not continue;
 206 //       //       cout << "Charged particle # " << " pos ("
 207 //       //        << xyzAtPHOS[0] << ", " << xyzAtPHOS[1] << ", " << xyzAtPHOS[2] << ")" <<  endl ;
 208 //       TVector3 poscp(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2]) ;
 209 //       Int_t cpN ;
 210 //       Double_t cpZ,cpX ;
 211 //       gime->PHOSGeometry()->ImpactOnEmc(poscp, cpN, cpZ, cpX) ;
 212 //       if (cpN) {// we are inside the PHOS acceptance
 213 //      //      cout << "Charged Matching 1: " << cpN << " " << cpZ << " " << cpX << endl ;
 214 //      //      cout << "Charged Matching 2: " << phN << " " << phZ << " " << phX << endl ;
 215 //      dist[0] = TMath::Sqrt( (cpZ-phZ)*(cpZ-phZ) + (cpX-phX)*(cpX-phX)) ;
 216 //       }
 217 //       phTheta = pp->Theta() ;
 218 //       phPhi   = pp->Phi() ;
 219 //       TParticle tempo ;
 220 //       tempo.SetMomentum(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2], 0.) ;
 221 //       cpTheta = tempo.Theta() ;
 222 //       cpPhi   = tempo.Phi() ;
 223 //       //cout << phTheta << " " << phPhi << " " << endl
 224 //       //cout <<       cpTheta << " " << cpPhi-phPhi << " " << endl ;
 225 //     }
 226 //     dist[1] = (phTheta - cpTheta)*TMath::RadToDeg() ;
 227 //     dist[2] = (phPhi - cpPhi)*TMath::RadToDeg() ;
 228 //   }
 229
 230 //   if (ep) {
 231 //     //cout << "EMCAL particle # " << endl ;
 232 //   }
 233 }