added low and high flux parameter options
[u/mrichter/AliRoot.git] / EMCAL / beamtest07 / TestESDCaloCluster.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
17 //AliRoot include files
18 #include "AliESD.h"
19 #include "AliESDEvent.h"
20 #include "AliEMCALLoader.h"
21 #include "AliESDCaloCluster.h"
22 #include "AliEMCALRecPoint.h"
23 #include "AliPID.h"
24 #include "AliLog.h" 
25
26 //macros
27
28 #endif
29
30 TChain * AliReadESDfromdisk(const UInt_t eventsToRead, 
31                                    const TString dirName, 
32                                    const TString esdTreeName, 
33                                    const char *  pattern) 
34 {
35   // Reads ESDs from Disk
36  TChain *  rv = 0; 
37   
38   // create a TChain of all the files 
39   TChain * cESDTree = new TChain(esdTreeName) ; 
40
41   // read from the directory file until the require number of events are collected
42   void * from = gSystem->OpenDirectory(dirName) ;
43    if (!from) 
44      rv = 0 ;   
45   else{ // reading file names from directory
46     const char * subdir ; 
47     // search all subdirectories witch matching pattern
48     while( (subdir = gSystem->GetDirEntry(from))  && 
49            (cESDTree->GetEntries() < eventsToRead)) {
50       if ( strstr(subdir, pattern) != 0 ) { 
51         char file[200] ; 
52         sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir);     
53         cESDTree->Add(file) ;
54       }
55     } // while file
56   
57     rv = cESDTree ; 
58     
59   } // reading file names from directory
60   return rv ; 
61 }
62
63 //======================================================================
64 TChain * AliReadESD(const UInt_t eventsToRead,
65                   const TString dirName, 
66                   const TString esdTreeName, 
67                   const char *  pattern)  
68 {
69   // Read AliESDs files and return a Chain of events
70  
71   if ( dirName == "" ) 
72     return 0 ; 
73   if ( esdTreeName == "" ) 
74     return AliReadESDfromdisk(eventsToRead, dirName,"","") ;//Last 2 arguments are not necessary but pdsf compiler complains "","'
75   else if ( strcmp(pattern, "") == 0 )
76     return AliReadESDfromdisk(eventsToRead, dirName, esdTreeName,"") ;//Last argument is not necessary but pdsf compiler complains "","'
77   else 
78     return AliReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ;        
79 }
80
81 //=====================================================================
82 //  Do:
83 //  .L TestESDCaloCluster.C++
84 //  TestESDCaloCluster(number of events to process)
85 //=====================================================================
86 void TestESDCaloCluster(const UInt_t eventsToProcess = 5, 
87                         TString dirName = "./", 
88                         const TString esdTreeName = "esdTree", 
89                         const char *  pattern = ".") 
90
91   
92   //Create chain of esd trees
93   //AliReadESD(eventsToProcess, directoryName,esdTreeName,patternOfDirectory) ;
94   //By default the root files are in the same directory 
95   TChain * t = AliReadESD(eventsToProcess, dirName,esdTreeName,pattern) ; 
96  
97   
98   // ESD
99   AliESDEvent * esd = new AliESDEvent();
100   esd->ReadFromTree(t);
101   
102   //Define few variables to be used in macro
103   TString alirunName = "" ; 
104
105   //Define example histograms
106   TH1F * hEnergy = new TH1F("hEnergy","Energy Distribution",100,0.,100.);   
107   TH1F * hEta    = new TH1F("hEta","Eta Distribution",100,-0.7,0.7);
108   TH1F * hPhi    = new TH1F("hPhi","Phi Distribution",100,0,2*TMath::Pi());
109
110   Int_t beg, end ;
111   UInt_t event ;
112   Float_t pos[3] ; 
113
114   for (event = 0; event < eventsToProcess; event++) {//event loop
115     //AliInfo( Form("Event %d \n",event) );  
116     Int_t nbytes = t->GetEntry(event); // store event in esd
117     //cout<<"nbytes "<<nbytes<<endl;
118     if ( nbytes == 0 ) //If nothing in ESD 
119       break ; 
120     
121     // Check that name of file is correct
122     if (alirunName != t->GetFile()->GetName()) {        
123       alirunName = t->GetFile()->GetName() ; 
124       alirunName.ReplaceAll("galice.root", "AliESDs.root") ;
125     }
126     
127     //get reconstructed vertex position 
128     
129     //Double_t vertex_position[3] ; 
130     //    esd->GetVertex()->GetXYZ(vertex_position) ; 
131     
132     cout<<"Event >>>>>>>>>>> "<<event<<endl;
133     
134
135     //select EMCAL tracks only 
136     end = esd->GetNumberOfCaloClusters() ;  
137     beg = 0;
138     Int_t nphP ;
139     //cout<<"begin "<<beg<<" end "<<end<<endl;
140     Int_t iclus = 0;
141     for (nphP =  beg; nphP <  end; nphP++) {//////////////EMCAL cluster loop
142       AliESDCaloCluster * clus = esd->GetCaloCluster(nphP) ; // retrieve cluster from esd
143       
144       //Get the cluster parameters
145       Float_t energy   = clus->E() ;
146       cout << "Cluster " << nphP << " energy = " << energy << endl;
147
148         //      Int_t mult       = clus->GetNumberOfDigits() ;
149         //      Float_t disp     = clus->GetClusterDisp() ;
150         //      UShort_t * amp   = clus->GetDigitAmplitude() ;
151         //      UShort_t * time  = clus->GetDigitTime() ;
152         //      UShort_t * index = clus->GetDigitIndex() ;
153         //      Int_t iprim = clus->GetPrimaryIndex();
154
155       /*
156       clus->GetPosition(pos) ;
157       TVector3 vpos(pos[0],pos[1],pos[2]) ;
158       //Print values on screen
159       cout<<"ESD cluster "<<iclus <<"; Energy "<<energy
160           <<"; Phi "<<vpos.Phi()*180/TMath::Pi()<<"; Eta "<<vpos.Eta()<<endl;
161       //cout<<"Dispersion "<<disp<<"; multiplicity "<<mult<<endl;
162       //Get the parent main values
163       */
164       //Fill histograms
165       hEnergy->Fill(energy) ;
166       //      hEta->Fill(vpos.Eta()) ;
167       //hPhi->Fill(vpos.Phi()) ;
168
169     }// track loop 
170   }//////////////////////////////////////////event loop
171   hEnergy->Draw();//Draw histogram on screen
172   //Write histograms in Root file 
173   TFile outf("histo.root","recreate") ;
174   hEnergy->Write() ;
175   hPhi->Write() ;
176   hEta->Write() ;
177   outf.Close() ;
178 }
179
180