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 "AliPHOSGetter.h"
  30 #include "AliPHOSGeometry.h"
  31 #include "Riostream.h"
  32 #include "AliESD.h"
  33 #include "AliESDtrack.h"
  34 #include "AliESDCaloTrack.h"
  35 #include "AliEMCALRecParticle.h"
  36 #include "AliPHOSRecParticle.h"
  37 #include "AliKalmanTrack.h"
  38 #include "AliPHOSGridFile.h"
  39 #endif
  40
  41 void Match(TParticle * pp, AliESDtrack * cp, Double_t * dist) ;
  42 TH1D * heta = new TH1D("heta", "Eta correlation", 100, 0., 360.) ;
  43 TH1D * hphi = new TH1D("hphi", "Phi correlation", 100, 0., 360.) ;
  44 // Data Challenge identification
  45 const TString kYear("2004") ;
  46 const TString kProd("02") ;
  47 const TString kVers("V4.01.Rev.00") ;
  48
  49 Bool_t Ana(const TString type = "per5", const Int_t run = 1, const Int_t nOfEvt = 1)
  50 {
  51   Double_t dist[3] ;
  52   // get the LFN file name in the Grid catalogue ;
  53   AliPHOSGridFile lfn ;
  54   if (!lfn.IsConnected())
  55     return kFALSE ;
  56   lfn.SetPath(kYear, kProd, kVers, type) ;
  57   lfn.SetRun(run) ;
  58
  59   //loop over the events
  60   Int_t nevt, evt = 0 ;
  61   for (nevt = 0 ; nevt < nOfEvt ; nevt++) {
  62     evt++ ;
  63     lfn.SetEvt(evt) ;
  64     TString fileName = lfn.GetLFN() ;
  65
  66     if (fileName.IsNull()) {
  67       nevt-- ;
  68       continue ;
  69     }
  70
  71     printf(">>>>>>>>>>>> Processing %s-%s/%s/%s : run # %d event # %d \n",
  72            kYear.Data(), kProd.Data(), kVers.Data(), type.Data(), run, evt) ;
  73     AliPHOSGetter * gime = AliPHOSGetter::Instance(fileName) ;
  74
  75     Int_t nEvent = gime->MaxEvent() ;
  76     Int_t event ;
  77     AliESD * esd = 0 ;
  78     for (event = 0 ; event < nEvent; event++) {
  79       esd = gime->ESD(event) ;
  80       if (!esd)
  81         return kFALSE ;
  82
  83       //esd->Print();
  84       Int_t caloindex ;
  85       // Calorimeter tracks
  86       AliESDCaloTrack * ct ;
  87       for (caloindex = 0 ; caloindex < esd->GetNumberOfCaloTracks() ; caloindex++) {
  88         // get the calorimeter type of particles (PHOS or EMCAL)
  89         ct = esd->GetCaloTrack(caloindex) ;
  90         TParticle * part = ct->GetRecParticle() ;
  91
  92         AliESDtrack * cp ;
  93         Int_t cpindex ;
  94         for (cpindex = 0 ; cpindex < esd->GetNumberOfTracks() ; cpindex++) {
  95           // get the charged tracks from central tracking
  96           cp = esd->GetTrack(cpindex) ;
  97           Match(part, cp, dist) ;
  98         }
  99         heta->Fill( dist[1] ) ;
 100         hphi->Fill( dist[2] ) ;
 101       }
 102     }
 103   }
 104   heta->Draw() ;
 105   //hphi->Draw() ;
 106   return kTRUE ;
 107 }
 108
 109 void Match(TParticle * part, AliESDtrack * cp, Double_t * dist)
 110 {
 111   // Calculates the distance (x,z) between  the particle detected by PHOS and
 112   // the charged particle reconstructed by the global tracking
 113
 114
 115   AliPHOSRecParticle * pp  = dynamic_cast<AliPHOSRecParticle*>(part) ;
 116   AliEMCALRecParticle * ep = dynamic_cast<AliEMCALRecParticle*>(part) ;
 117
 118   Int_t phN ;
 119   Double_t phZ, phX ;
 120
 121   if (pp) { // it is a PHOS particle
 122     Double_t cpTheta,  cpPhi ;
 123     Double_t phTheta,  phPhi ;
 124     cpTheta = cpPhi = phTheta = phPhi = 0. ;
 125    //    cout << "PHOS particle # " << " pos ("
 126     //   << pp->GetPos().X() << ", " << pp->GetPos().Y() << ", " << pp->GetPos().Z() << ")" << endl ;
 127
 128     AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 129     gime->PHOSGeometry()->ImpactOnEmc(*pp, phN, phZ, phX) ;
 130     Double_t xyzAtPHOS[3] ;
 131     cp->GetOuterXYZ(xyzAtPHOS) ;
 132     if ( (xyzAtPHOS[0] +  xyzAtPHOS[1] + xyzAtPHOS[2]) != 0.) { //it has reached PHOS
 133       //the next check are only if we want high quality tracks
 134       //       ULong_t status = cp->GetStatus() ;
 135       //       if ((status & AliESDtrack::kTRDput)==0)
 136       //        do not continue;
 137       //       if ((status & AliESDtrack::kTRDStop)!=0)
 138       //        do not continue;
 139       //       cout << "Charged particle # " << " pos ("
 140       //           << xyzAtPHOS[0] << ", " << xyzAtPHOS[1] << ", " << xyzAtPHOS[2] << ")" <<  endl ;
 141       TVector3 poscp(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2]) ;
 142       Int_t cpN ;
 143       Double_t cpZ,cpX ;
 144       gime->PHOSGeometry()->ImpactOnEmc(poscp, cpN, cpZ, cpX) ;
 145       if (cpN) {// we are inside the PHOS acceptance
 146         //      cout << "Charged Matching 1: " << cpN << " " << cpZ << " " << cpX << endl ;
 147         //      cout << "Charged Matching 2: " << phN << " " << phZ << " " << phX << endl ;
 148         dist[0] = TMath::Sqrt( (cpZ-phZ)*(cpZ-phZ) + (cpX-phX)*(cpX-phX)) ;
 149       }
 150       phTheta = pp->Theta() ;
 151       phPhi   = pp->Phi() ;
 152       TParticle tempo ;
 153       tempo.SetMomentum(xyzAtPHOS[0], xyzAtPHOS[1], xyzAtPHOS[2], 0.) ;
 154       cpTheta = tempo.Theta() ;
 155       cpPhi   = tempo.Phi() ;
 156       //cout << phTheta << " " << phPhi << " " << endl
 157       //cout <<  cpTheta << " " << cpPhi-phPhi << " " << endl ;
 158     }
 159     dist[1] = (phTheta - cpTheta)*TMath::RadToDeg() ;
 160     dist[2] = (phPhi - cpPhi)*TMath::RadToDeg() ;
 161   }
 162
 163   if (ep) {
 164     //cout << "EMCAL particle # " << endl ;
 165   }
 166 }