Centrality Task (A. Toia)
[u/mrichter/AliRoot.git] / ANALYSIS / AliCentralitySelectionTask.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, 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
16 //*****************************************************
17 //   Class AliCentralitySelectionTask
18 //   Class to analyze determine centrality            
19 //   author: Alberica Toia
20 //*****************************************************
21
22 #include <TTree.h>
23 #include <TList.h>
24 #include <TH1F.h>
25 #include <TH2F.h>
26 #include <TProfile.h>
27 #include <TFile.h>
28 #include <TString.h>
29 #include <TCanvas.h>
30 #include <iostream>
31
32 #include "AliAnalysisManager.h"
33 #include "AliVEvent.h"
34 #include "AliESD.h"
35 #include "AliESDEvent.h"
36 #include "AliESDHeader.h"
37 #include "AliESDInputHandler.h"
38 #include "AliESDZDC.h"
39 #include "AliESDFMD.h"
40 #include "AliESDVZERO.h"
41 #include "AliESDCentrality.h"
42 #include "AliMultiplicity.h"
43 #include "AliAODHandler.h"
44 #include "AliAODEvent.h"
45 #include "AliAODVertex.h"
46 #include "AliAODMCHeader.h"
47 #include "AliMCEvent.h"
48 #include "AliMCEventHandler.h"
49 #include "AliMCParticle.h"
50 #include "AliStack.h"
51 #include "AliHeader.h"
52 #include "AliAODMCParticle.h"
53 #include "AliAnalysisTaskSE.h"
54 #include "AliGenEventHeader.h"
55 #include "AliGenHijingEventHeader.h"
56 #include "AliPhysicsSelectionTask.h"
57 #include "AliPhysicsSelection.h"
58 #include "AliBackgroundSelection.h"
59 #include "AliCentralitySelectionTask.h"
60
61 ClassImp(AliCentralitySelectionTask)
62
63
64 //________________________________________________________________________
65 AliCentralitySelectionTask::AliCentralitySelectionTask():
66 AliAnalysisTaskSE(),
67   fDebug(0),
68   fAnalysisInput("ESD"),
69   fIsMCInput(kFALSE),
70   fFile(0),
71   fCentfilename(0),
72   fMethod(0),
73   fCent(0),
74   fHtemp(0)
75 {   
76   // Default constructor
77 }   
78
79 //________________________________________________________________________
80 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
81   AliAnalysisTaskSE(name),
82   fDebug(0),
83   fAnalysisInput("ESD"),
84   fIsMCInput(kFALSE),
85   fFile(0),
86   fCentfilename(0),
87   fMethod(0),
88   fCent(0),
89   fHtemp(0)
90 {
91   // Default constructor
92 }
93
94 //________________________________________________________________________
95 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
96 {
97   // Assignment operator
98   if (this!=&c) {
99     AliAnalysisTaskSE::operator=(c);
100   }
101   return *this;
102 }
103
104 //________________________________________________________________________
105 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
106   AliAnalysisTaskSE(ana),
107   fDebug(ana.fDebug),     
108   fAnalysisInput(ana.fDebug),
109   fIsMCInput(ana.fIsMCInput),
110   fFile(ana.fFile),
111   fCentfilename(ana.fCentfilename),
112   fMethod(ana.fMethod),
113   fCent(ana.fCent),
114   fHtemp(ana.fHtemp)
115 {
116   // Copy Constructor   
117 }
118  
119 //________________________________________________________________________
120  AliCentralitySelectionTask::~AliCentralitySelectionTask()
121  {
122    // Destructor
123  }  
124
125 //________________________________________________________________________
126 void AliCentralitySelectionTask::UserCreateOutputObjects()
127 {  
128   // Create the output containers
129   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
130 }
131
132 //________________________________________________________________________
133 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
134
135   // Execute analysis for current event:
136   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
137   
138   Float_t  zncEnergy;                 //  ZNC Energy
139   Float_t  zpcEnergy;                 //  ZPC Energy
140   Float_t  znaEnergy;                 //  ZNA Energy
141   Float_t  zpaEnergy;                 //  ZPA Energy
142   Float_t  zem1Energy;                //  ZEM1 Energy
143   Float_t  zem2Energy;                //  ZEM2 Energy
144   
145   Int_t    nTracks    = 0;            //  no. tracks
146   Int_t    nTracklets = 0;            //  no. tracklets
147   Int_t    nClusters[6];              //  no. clusters on 6 ITS layers
148   Int_t    nChips[2];                 //  no. chips on 2 SPD layers
149
150   Float_t  multV0A    = 0;            //  multiplicity from V0 reco side A
151   Float_t  multV0C    = 0;            //  multiplicity from V0 reco side C
152   Float_t  multFMDA   = 0;            //  multiplicity from FMD on detector A
153   Float_t  multFMDC   = 0;            //  multiplicity from FMD on detector C
154
155   AliESDCentrality *esdCent = 0;
156
157   if(fAnalysisInput.CompareTo("ESD")==0){
158
159     AliVEvent* event = InputEvent();
160     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
161     
162     esdCent = esd->GetCentrality();
163
164     // ***** V0 info    
165     AliESDVZERO* esdV0 = esd->GetVZEROData();
166     multV0A=esdV0->GetMTotV0A();
167     multV0C=esdV0->GetMTotV0C();
168     
169     // ***** CB info (tracklets, clusters, chips)
170     nTracks    = event->GetNumberOfTracks();     
171
172     const AliMultiplicity *mult = esd->GetMultiplicity();
173
174     nTracklets = mult->GetNumberOfTracklets();
175
176     for(Int_t ilay=0; ilay<6; ilay++){
177       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
178     }
179
180     for(Int_t ilay=0; ilay<2; ilay++){
181       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
182     }
183     
184
185     // ***** FMD info
186     AliESDFMD *fmd = esd->GetFMDData();
187     Float_t totalMultA = 0;
188     Float_t totalMultC = 0;
189     const Float_t fFMDLowCut = 0.4;
190     
191     for(UShort_t det=1;det<=3;det++) {
192       Int_t nRings = (det==1 ? 1 : 2);
193       for (UShort_t ir = 0; ir < nRings; ir++) {          
194         Char_t   ring = (ir == 0 ? 'I' : 'O');
195         UShort_t nsec = (ir == 0 ? 20  : 40);
196         UShort_t nstr = (ir == 0 ? 512 : 256);
197         for(UShort_t sec =0; sec < nsec;  sec++)  {
198           for(UShort_t strip = 0; strip < nstr; strip++) {
199             
200             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
201             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
202             
203             Float_t nParticles=0;
204             
205             if(FMDmult > fFMDLowCut) {
206               nParticles = 1.;
207             }
208             
209             if (det<3) totalMultA = totalMultA + nParticles;
210             else totalMultC = totalMultC + nParticles;
211             
212           }
213         }
214       }
215     }
216     multFMDA = totalMultA;
217     multFMDC = totalMultC;
218     
219     // ***** ZDC info
220     AliESDZDC *esdZDC = esd->GetESDZDC();
221     zncEnergy  = (Float_t) (esdZDC->GetZDCN1Energy());
222     zpcEnergy  = (Float_t) (esdZDC->GetZDCP1Energy());
223     znaEnergy  = (Float_t) (esdZDC->GetZDCN2Energy());
224     zpaEnergy  = (Float_t) (esdZDC->GetZDCP2Energy());
225     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
226     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
227     
228   }   
229   else if(fAnalysisInput.CompareTo("AOD")==0){
230     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
231     // to be implemented
232     printf("  AOD analysis not yet implemented!!!\n\n");
233     return;
234   }
235
236    // ***** Centrality Selection
237   if(fMethod.CompareTo("V0")==0){
238     fCent = fHtemp->GetBinContent(fHtemp->FindBin((multV0A+multV0C)));
239   }
240   if(fMethod.CompareTo("FMD")==0){
241     fCent = fHtemp->GetBinContent(fHtemp->FindBin((multFMDA+multFMDC)));
242   }
243   if(fMethod.CompareTo("TRACKS")==0) {
244     fCent = fHtemp->GetBinContent(fHtemp->FindBin(nTracks));
245   }
246   if(fMethod.CompareTo("TRACKLETS")==0){
247     fCent = fHtemp->GetBinContent(fHtemp->FindBin(nTracklets));
248   }
249   if(fMethod.CompareTo("CLUSTERS")==0) {
250     fCent = fHtemp->GetBinContent(fHtemp->FindBin(nClusters[0]));
251   }
252   printf(" **** centrality is %3.2f \n", fCent);
253   
254   esdCent->SetCentrality(fCent);
255 }
256
257 //________________________________________________________________________
258 void AliCentralitySelectionTask::SetCentralityMethod(const char* x) 
259 {
260   fMethod = x;
261
262   fFile  = new TFile(fCentfilename);
263   
264   if(fMethod.CompareTo("V0")==0)
265     fHtemp  = (TH1D*) (fFile->Get("hmultV0_percentile")); 
266   
267   if(fMethod.CompareTo("FMD")==0)
268     fHtemp  = (TH1D*) (fFile->Get("hmultFMD_percentile")); 
269   
270   if(fMethod.CompareTo("TRACKS")==0) 
271     fHtemp  = (TH1D*) (fFile->Get("hNtracks_percentile")); 
272   
273   if(fMethod.CompareTo("TRACKLETS")==0)
274     fHtemp  = (TH1D*) (fFile->Get("hNtracklets_percentile")); 
275   
276   if(fMethod.CompareTo("CLUSTERS")==0) 
277     fHtemp  = (TH1D*) (fFile->Get("hNclusters0_percentile")); 
278   
279 }
280
281 //________________________________________________________________________
282 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
283 {
284   // Terminate analysis
285   fFile->Close();  
286 }
287 //________________________________________________________________________
288