.so cleanup: more gSystem->Load()
[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 "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 }