updating
[u/mrichter/AliRoot.git] / 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 }