Corrections Ch. Kleinn-Boesing
[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 <TObjString.h>
29 #include <TString.h>
30 #include <TCanvas.h>
31 #include <TROOT.h>
32 #include <TDirectory.h>
33 #include <iostream>
34
35 #include "AliAnalysisManager.h"
36 #include "AliVEvent.h"
37 #include "AliESD.h"
38 #include "AliESDEvent.h"
39 #include "AliESDHeader.h"
40 #include "AliESDInputHandler.h"
41 #include "AliESDZDC.h"
42 #include "AliESDFMD.h"
43 #include "AliESDVZERO.h"
44 #include "AliESDCentrality.h"
45 #include "AliMultiplicity.h"
46 #include "AliAODHandler.h"
47 #include "AliAODEvent.h"
48 #include "AliAODVertex.h"
49 #include "AliAODMCHeader.h"
50 #include "AliMCEvent.h"
51 #include "AliMCEventHandler.h"
52 #include "AliMCParticle.h"
53 #include "AliStack.h"
54 #include "AliHeader.h"
55 #include "AliAODMCParticle.h"
56 #include "AliAnalysisTaskSE.h"
57 #include "AliGenEventHeader.h"
58 #include "AliGenHijingEventHeader.h"
59 #include "AliPhysicsSelectionTask.h"
60 #include "AliPhysicsSelection.h"
61 #include "AliBackgroundSelection.h"
62 #include "AliCentralitySelectionTask.h"
63
64 ClassImp(AliCentralitySelectionTask)
65
66
67 //________________________________________________________________________
68 AliCentralitySelectionTask::AliCentralitySelectionTask():
69 AliAnalysisTaskSE(),
70   fDebug(0),
71   fAnalysisInput("ESD"),
72   fIsMCInput(kFALSE),
73   fFile(0),
74   fFile2(0),
75   fCentfilename(""),
76   fCentfilename2(""),
77   fFileList(new TList),
78   fFileList2(new TList),
79   fCurrentRun(-1),
80   fCentV0M(0),
81   fCentFMD(0),
82   fCentTRK(0),
83   fCentTKL(0),
84   fCentCL0(0),
85   fCentCL1(0),
86   fCentV0MvsFMD(0),
87   fCentTKLvsV0M(0),
88   fCentZEMvsZDC(0),
89   fHtempV0M(0),
90   fHtempFMD(0),
91   fHtempTRK(0),
92   fHtempTKL(0),
93   fHtempCL0(0),
94   fHtempCL1(0),
95   fHtempV0MvsFMD(0),
96   fHtempTKLvsV0M(0),
97   fHtempZEMvsZDC(0)
98 {   
99   // Default constructor
100   fFileList->SetOwner();
101   fFileList2->SetOwner();
102
103   AliInfo("Centrality Selection enabled.");
104 }   
105
106 //________________________________________________________________________
107 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
108   AliAnalysisTaskSE(name),
109   fDebug(0),
110   fAnalysisInput("ESD"),
111   fIsMCInput(kFALSE),
112   fFile(0),
113   fFile2(0),
114   fCentfilename(""),
115   fCentfilename2(""),
116   fFileList(new TList),
117   fFileList2(new TList),
118   fCurrentRun(-1),
119   fCentV0M(0),
120   fCentFMD(0),
121   fCentTRK(0),
122   fCentTKL(0),
123   fCentCL0(0),
124   fCentCL1(0),
125   fCentV0MvsFMD(0),
126   fCentTKLvsV0M(0),
127   fCentZEMvsZDC(0),
128   fHtempV0M(0),
129   fHtempFMD(0),
130   fHtempTRK(0),
131   fHtempTKL(0),
132   fHtempCL0(0),
133   fHtempCL1(0),
134   fHtempV0MvsFMD(0),
135   fHtempTKLvsV0M(0),
136   fHtempZEMvsZDC(0)
137 {
138   // Default constructor
139   fFileList->SetOwner();
140   fFileList2->SetOwner();
141
142   AliInfo("Centrality Selection enabled.");
143 }
144
145 //________________________________________________________________________
146 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
147 {
148   // Assignment operator
149   if (this!=&c) {
150     AliAnalysisTaskSE::operator=(c);
151   }
152   return *this;
153 }
154
155 //________________________________________________________________________
156 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
157   AliAnalysisTaskSE(ana),
158   fDebug(ana.fDebug),     
159   fAnalysisInput(ana.fDebug),
160   fIsMCInput(ana.fIsMCInput),
161   fFile(ana.fFile),
162   fFile2(ana.fFile2),
163   fCentfilename(ana.fCentfilename),
164   fCentfilename2(ana.fCentfilename2),
165   fFileList(ana.fFileList),
166   fFileList2(ana.fFileList2),
167   fCurrentRun(ana.fCurrentRun),
168   fCentV0M(ana.fCentV0M),
169   fCentFMD(ana.fCentFMD),
170   fCentTRK(ana.fCentTRK),
171   fCentTKL(ana.fCentTKL),
172   fCentCL0(ana.fCentCL0),
173   fCentCL1(ana.fCentCL1),
174   fCentV0MvsFMD(ana.fCentV0MvsFMD),
175   fCentTKLvsV0M(ana.fCentTKLvsV0M),
176   fCentZEMvsZDC(ana.fCentZEMvsZDC),
177   fHtempV0M(ana.fHtempV0M),
178   fHtempFMD(ana.fHtempFMD),
179   fHtempTRK(ana.fHtempTRK),
180   fHtempTKL(ana.fHtempTKL),
181   fHtempCL0(ana.fHtempCL0),
182   fHtempCL1(ana.fHtempCL1),
183   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
184   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
185   fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
186 {
187   // Copy Constructor   
188 }
189  
190 //________________________________________________________________________
191 AliCentralitySelectionTask::~AliCentralitySelectionTask()
192 {
193   // Destructor
194   
195   if (fFileList) {
196     fFileList->Clear();
197     delete fFileList;
198   }
199   fFileList = NULL;
200
201   if (fFileList2) {
202     fFileList2->Clear();
203     delete fFileList2;
204   }
205   fFileList2 = NULL;
206 }  
207
208 //________________________________________________________________________
209 void AliCentralitySelectionTask::UserCreateOutputObjects()
210 {  
211   // Create the output containers
212   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
213
214   if (fFileList->GetEntries() < 1) {
215     AliError("No Inputfiles Added");
216   }
217
218   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
219 }
220
221 //________________________________________________________________________
222 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
223
224   // Execute analysis for current event:
225   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
226   
227   Float_t  zncEnergy;          //  ZNC Energy
228   Float_t  zpcEnergy;          //  ZPC Energy
229   Float_t  znaEnergy;          //  ZNA Energy
230   Float_t  zpaEnergy;          //  ZPA Energy
231   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
232   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
233   
234   Int_t    nTracks = 0;             //  no. tracks
235   Int_t    nTracklets = 0;          //  no. tracklets
236   Int_t    nClusters[6];            //  no. clusters on 6 ITS layers
237   Int_t    nChips[2];               //  no. chips on 2 SPD layers
238
239   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
240   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
241   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
242   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
243
244   AliESDCentrality *esdCent = 0;
245
246   if(fAnalysisInput.CompareTo("ESD")==0){
247
248     AliVEvent* event = InputEvent();
249     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
250
251     if (SetupRun(esd))   
252       return;
253
254     esdCent = esd->GetCentrality();
255
256     // ***** V0 info    
257     AliESDVZERO* esdV0 = esd->GetVZEROData();
258     multV0A=esdV0->GetMTotV0A();
259     multV0C=esdV0->GetMTotV0C();
260     
261     // ***** CB info (tracklets, clusters, chips)
262     nTracks    = event->GetNumberOfTracks();     
263
264     const AliMultiplicity *mult = esd->GetMultiplicity();
265
266     nTracklets = mult->GetNumberOfTracklets();
267
268     for(Int_t ilay=0; ilay<6; ilay++){
269       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
270     }
271
272     for(Int_t ilay=0; ilay<2; ilay++){
273       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
274     }
275     
276
277     // ***** FMD info
278     AliESDFMD *fmd = esd->GetFMDData();
279     Float_t totalMultA = 0;
280     Float_t totalMultC = 0;
281     const Float_t fFMDLowCut = 0.4;
282     
283     for(UShort_t det=1;det<=3;det++) {
284       Int_t nRings = (det==1 ? 1 : 2);
285       for (UShort_t ir = 0; ir < nRings; ir++) {          
286         Char_t   ring = (ir == 0 ? 'I' : 'O');
287         UShort_t nsec = (ir == 0 ? 20  : 40);
288         UShort_t nstr = (ir == 0 ? 512 : 256);
289         for(UShort_t sec =0; sec < nsec;  sec++)  {
290           for(UShort_t strip = 0; strip < nstr; strip++) {
291             
292             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
293             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
294             
295             Float_t nParticles=0;
296             
297             if(FMDmult > fFMDLowCut) {
298               nParticles = 1.;
299             }
300             
301             if (det<3) totalMultA = totalMultA + nParticles;
302             else totalMultC = totalMultC + nParticles;
303             
304           }
305         }
306       }
307     }
308     multFMDA = totalMultA;
309     multFMDC = totalMultC;
310     
311     // ***** ZDC info
312     AliESDZDC *esdZDC = esd->GetESDZDC();
313     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
314     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
315     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
316     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
317     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
318     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
319     
320   }   
321   else if(fAnalysisInput.CompareTo("AOD")==0){
322     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
323     // to be implemented
324     printf("  AOD analysis not yet implemented!!!\n\n");
325     return;
326   }
327
328   // ***** Centrality Selection
329   if(fHtempV0M)  fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0A+multV0C)));
330   ///  else     printf("  Centrality by V0 not available!!!\n\n");
331   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
332   //  else     printf("  Centrality by FMD not available!!!\n\n");
333   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
334   //  else     printf("  Centrality by TRK not available!!!\n\n");
335   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
336   //  else     printf("  Centrality by TKL not available!!!\n\n");
337   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
338   //  else     printf("  Centrality by CL0 not available!!!\n\n");
339   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(nClusters[1]));
340   ///  else     printf("  Centrality by CL1 not available!!!\n\n");
341   
342   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
343   //  else     printf("  Centrality by V0 vs FMD not available!!!\n\n");
344   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
345   //  else     printf("  Centrality by V0 vs TKL not available!!!\n\n");
346   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
347   //  else     printf("  Centrality by ZEM vs ZDC not available!!!\n\n");
348
349   esdCent->SetCentralityV0M(fCentV0M);
350   esdCent->SetCentralityFMD(fCentFMD);
351   esdCent->SetCentralityTRK(fCentTRK);
352   esdCent->SetCentralityTKL(fCentTKL);
353   esdCent->SetCentralityCL0(fCentCL0);
354   esdCent->SetCentralityCL1(fCentCL1);
355   esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
356   esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
357   esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
358 }
359
360 //________________________________________________________________________
361 void AliCentralitySelectionTask::ReadCentralityHistos() 
362 {
363 //  Read centrality histograms
364     TDirectory *owd = gDirectory;
365     fFile  = TFile::Open(fCentfilename);
366     owd->cd();
367     fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
368     fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
369     fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
370     fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
371     fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
372     fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
373 }  
374
375 //________________________________________________________________________
376 void AliCentralitySelectionTask::ReadCentralityHistos2() 
377 {
378 //  Read centrality histograms
379     TDirectory *owd = gDirectory;
380     fFile2  = TFile::Open(fCentfilename2);
381     owd->cd();
382     fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
383     fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
384     fHtempZEMvsZDC  = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
385 }
386
387 //________________________________________________________________________
388 void AliCentralitySelectionTask::SetPercentileFile(TString filename) 
389 {
390 // Set the percentile file name
391   fCentfilename = filename;
392   ReadCentralityHistos();
393 }
394
395 //________________________________________________________________________
396 void AliCentralitySelectionTask::SetPercentileFile2(TString filename) 
397 {
398 // Set the percentile file name
399   fCentfilename2 = filename;
400   ReadCentralityHistos2();
401 }
402
403 //________________________________________________________________________
404 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
405 {
406   // Terminate analysis
407   if (fFile && fFile->IsOpen())
408     fFile->Close();  
409   if (fFile2 && fFile2->IsOpen())
410     fFile2->Close();  
411 }
412 //________________________________________________________________________
413 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)
414 {
415   // Setup files for run
416
417   if (!esd)
418     return -1;
419
420   // check if something to be done
421   if (fCurrentRun == esd->GetRunNumber())
422     return 0;
423   else
424     fCurrentRun = esd->GetRunNumber();
425   
426   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
427
428   Int_t runNo = fCurrentRun;
429
430   // CHANGE HERE FOR RUN RANGES
431   if ( runNo == 137162 ) runNo = 137161;
432   else if ( runNo == 137365 ) runNo = 137366;
433   // CHANGE HERE FOR RUN RANGES
434
435   TString runName(Form("%d", runNo));
436   TString fileName("");
437   Bool_t isRunKnown = kFALSE;
438
439   // Check if run is in fileList
440   // if not, take the last name in the list
441   for ( Int_t idx=0 ; idx < fFileList->GetEntries(); ++idx ) {
442
443     TString str((dynamic_cast<TObjString*>(fFileList->At(idx)))->GetString());
444     if (str.Contains(runName)) {
445       fileName += str;
446       isRunKnown = kTRUE;
447       break;
448     }
449   }
450
451   if (!isRunKnown) {
452     if (fFileList->Last()) {
453       fileName += (dynamic_cast<TObjString*>(fFileList->Last()))->GetString();
454       AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
455     }
456   }
457
458   if (fileName.Contains(".root")) {
459     AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
460     SetPercentileFile(fileName.Data());
461     return 0;
462   }
463   
464   return -1;
465 }