1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //*****************************************************
17 // Class AliCentralitySelectionTask
18 // Class to analyze determine centrality
19 // author: Alberica Toia
20 //*****************************************************
28 #include <TObjString.h>
33 #include "AliAnalysisManager.h"
34 #include "AliVEvent.h"
36 #include "AliESDEvent.h"
37 #include "AliESDHeader.h"
38 #include "AliESDInputHandler.h"
39 #include "AliESDZDC.h"
40 #include "AliESDFMD.h"
41 #include "AliESDVZERO.h"
42 #include "AliESDCentrality.h"
43 #include "AliMultiplicity.h"
44 #include "AliAODHandler.h"
45 #include "AliAODEvent.h"
46 #include "AliAODVertex.h"
47 #include "AliAODMCHeader.h"
48 #include "AliMCEvent.h"
49 #include "AliMCEventHandler.h"
50 #include "AliMCParticle.h"
52 #include "AliHeader.h"
53 #include "AliAODMCParticle.h"
54 #include "AliAnalysisTaskSE.h"
55 #include "AliGenEventHeader.h"
56 #include "AliGenHijingEventHeader.h"
57 #include "AliPhysicsSelectionTask.h"
58 #include "AliPhysicsSelection.h"
59 #include "AliBackgroundSelection.h"
60 #include "AliCentralitySelectionTask.h"
62 ClassImp(AliCentralitySelectionTask)
65 //________________________________________________________________________
66 AliCentralitySelectionTask::AliCentralitySelectionTask():
69 fAnalysisInput("ESD"),
76 fFileList2(new TList),
97 // Default constructor
98 fFileList->SetOwner();
99 fFileList2->SetOwner();
101 AliInfo("Centrality Selection enabled.");
104 //________________________________________________________________________
105 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
106 AliAnalysisTaskSE(name),
108 fAnalysisInput("ESD"),
114 fFileList(new TList),
115 fFileList2(new TList),
136 // Default constructor
137 fFileList->SetOwner();
138 fFileList2->SetOwner();
140 AliInfo("Centrality Selection enabled.");
143 //________________________________________________________________________
144 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
146 // Assignment operator
148 AliAnalysisTaskSE::operator=(c);
153 //________________________________________________________________________
154 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
155 AliAnalysisTaskSE(ana),
157 fAnalysisInput(ana.fDebug),
158 fIsMCInput(ana.fIsMCInput),
161 fCentfilename(ana.fCentfilename),
162 fCentfilename2(ana.fCentfilename2),
163 fFileList(ana.fFileList),
164 fFileList2(ana.fFileList2),
165 fCurrentRun(ana.fCurrentRun),
166 fCentV0M(ana.fCentV0M),
167 fCentFMD(ana.fCentFMD),
168 fCentTRK(ana.fCentTRK),
169 fCentTKL(ana.fCentTKL),
170 fCentCL0(ana.fCentCL0),
171 fCentCL1(ana.fCentCL1),
172 fCentV0MvsFMD(ana.fCentV0MvsFMD),
173 fCentTKLvsV0M(ana.fCentTKLvsV0M),
174 fCentZEMvsZDC(ana.fCentZEMvsZDC),
175 fHtempV0M(ana.fHtempV0M),
176 fHtempFMD(ana.fHtempFMD),
177 fHtempTRK(ana.fHtempTRK),
178 fHtempTKL(ana.fHtempTKL),
179 fHtempCL0(ana.fHtempCL0),
180 fHtempCL1(ana.fHtempCL1),
181 fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
182 fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
183 fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
188 //________________________________________________________________________
189 AliCentralitySelectionTask::~AliCentralitySelectionTask()
206 //________________________________________________________________________
207 void AliCentralitySelectionTask::UserCreateOutputObjects()
209 // Create the output containers
210 if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
212 if (fFileList->GetEntries() < 1) {
213 AliError("No Inputfiles Added");
216 AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
219 //________________________________________________________________________
220 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
222 // Execute analysis for current event:
223 if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
225 Float_t zncEnergy; // ZNC Energy
226 Float_t zpcEnergy; // ZPC Energy
227 Float_t znaEnergy; // ZNA Energy
228 Float_t zpaEnergy; // ZPA Energy
229 Float_t zem1Energy = 0.; // ZEM1 Energy
230 Float_t zem2Energy = 0.; // ZEM2 Energy
232 Int_t nTracks = 0; // no. tracks
233 Int_t nTracklets = 0; // no. tracklets
234 Int_t nClusters[6]; // no. clusters on 6 ITS layers
235 Int_t nChips[2]; // no. chips on 2 SPD layers
237 Float_t multV0A = 0; // multiplicity from V0 reco side A
238 Float_t multV0C = 0; // multiplicity from V0 reco side C
239 Float_t multFMDA = 0; // multiplicity from FMD on detector A
240 Float_t multFMDC = 0; // multiplicity from FMD on detector C
242 AliESDCentrality *esdCent = 0;
244 if(fAnalysisInput.CompareTo("ESD")==0){
246 AliVEvent* event = InputEvent();
247 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
252 esdCent = esd->GetCentrality();
255 AliESDVZERO* esdV0 = esd->GetVZEROData();
256 multV0A=esdV0->GetMTotV0A();
257 multV0C=esdV0->GetMTotV0C();
259 // ***** CB info (tracklets, clusters, chips)
260 nTracks = event->GetNumberOfTracks();
262 const AliMultiplicity *mult = esd->GetMultiplicity();
264 nTracklets = mult->GetNumberOfTracklets();
266 for(Int_t ilay=0; ilay<6; ilay++){
267 nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
270 for(Int_t ilay=0; ilay<2; ilay++){
271 nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
276 AliESDFMD *fmd = esd->GetFMDData();
277 Float_t totalMultA = 0;
278 Float_t totalMultC = 0;
279 const Float_t fFMDLowCut = 0.4;
281 for(UShort_t det=1;det<=3;det++) {
282 Int_t nRings = (det==1 ? 1 : 2);
283 for (UShort_t ir = 0; ir < nRings; ir++) {
284 Char_t ring = (ir == 0 ? 'I' : 'O');
285 UShort_t nsec = (ir == 0 ? 20 : 40);
286 UShort_t nstr = (ir == 0 ? 512 : 256);
287 for(UShort_t sec =0; sec < nsec; sec++) {
288 for(UShort_t strip = 0; strip < nstr; strip++) {
290 Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
291 if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
293 Float_t nParticles=0;
295 if(FMDmult > fFMDLowCut) {
299 if (det<3) totalMultA = totalMultA + nParticles;
300 else totalMultC = totalMultC + nParticles;
306 multFMDA = totalMultA;
307 multFMDC = totalMultC;
310 AliESDZDC *esdZDC = esd->GetESDZDC();
311 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
312 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
313 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
314 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
315 zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
316 zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
319 else if(fAnalysisInput.CompareTo("AOD")==0){
320 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
322 printf(" AOD analysis not yet implemented!!!\n\n");
326 // ***** Centrality Selection
327 if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0A+multV0C)));
328 /// else printf(" Centrality by V0 not available!!!\n\n");
329 if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
330 // else printf(" Centrality by FMD not available!!!\n\n");
331 if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
332 // else printf(" Centrality by TRK not available!!!\n\n");
333 if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
334 // else printf(" Centrality by TKL not available!!!\n\n");
335 if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
336 // else printf(" Centrality by CL0 not available!!!\n\n");
337 if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(nClusters[1]));
338 /// else printf(" Centrality by CL1 not available!!!\n\n");
340 if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
341 // else printf(" Centrality by V0 vs FMD not available!!!\n\n");
342 if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
343 // else printf(" Centrality by V0 vs TKL not available!!!\n\n");
344 if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
345 // else printf(" Centrality by ZEM vs ZDC not available!!!\n\n");
347 esdCent->SetCentralityV0M(fCentV0M);
348 esdCent->SetCentralityFMD(fCentFMD);
349 esdCent->SetCentralityTRK(fCentTRK);
350 esdCent->SetCentralityTKL(fCentTKL);
351 esdCent->SetCentralityCL0(fCentCL0);
352 esdCent->SetCentralityCL1(fCentCL1);
353 esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
354 esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
355 esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
358 //________________________________________________________________________
359 void AliCentralitySelectionTask::ReadCentralityHistos()
361 fFile = TFile::Open(fCentfilename,"READ");
362 fHtempV0M = (TH1D*) (fFile->Get("hmultV0_percentile"));
363 fHtempFMD = (TH1D*) (fFile->Get("hmultFMD_percentile"));
364 fHtempTRK = (TH1D*) (fFile->Get("hNtracks_percentile"));
365 fHtempTKL = (TH1D*) (fFile->Get("hNtracklets_percentile"));
366 fHtempCL0 = (TH1D*) (fFile->Get("hNclusters0_percentile"));
367 fHtempCL1 = (TH1D*) (fFile->Get("hNclusters1_percentile"));
370 //________________________________________________________________________
371 void AliCentralitySelectionTask::ReadCentralityHistos2()
373 fFile2 = TFile::Open(fCentfilename2);
374 fHtempV0MvsFMD = (TH1D*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
375 fHtempTKLvsV0M = (TH1D*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
376 fHtempZEMvsZDC = (TH1D*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
379 //________________________________________________________________________
380 void AliCentralitySelectionTask::SetPercentileFile(TString filename)
382 fCentfilename = filename;
383 ReadCentralityHistos();
386 //________________________________________________________________________
387 void AliCentralitySelectionTask::SetPercentileFile2(TString filename)
389 fCentfilename2 = filename;
390 ReadCentralityHistos2();
393 //________________________________________________________________________
394 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
396 // Terminate analysis
397 if (fFile && fFile->IsOpen())
399 if (fFile2 && fFile2->IsOpen())
402 //________________________________________________________________________
403 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)
405 // Setup files for run
410 // check if something to be done
411 if (fCurrentRun == esd->GetRunNumber())
414 fCurrentRun = esd->GetRunNumber();
416 AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
418 Int_t runNo = fCurrentRun;
420 // CHANGE HERE FOR RUN RANGES
421 if ( runNo == 137162 ) runNo = 137161;
422 else if ( runNo == 137365 ) runNo = 137366;
423 // CHANGE HERE FOR RUN RANGES
425 TString runName(Form("%d", runNo));
426 TString fileName("");
427 Bool_t isRunKnown = kFALSE;
429 // Check if run is in fileList
430 // if not, take the last name in the list
431 for ( Int_t idx=0 ; idx < fFileList->GetEntries(); ++idx ) {
433 TString str((dynamic_cast<TObjString*>(fFileList->At(idx)))->GetString());
434 if (str.Contains(runName)) {
442 if (fFileList->Last()) {
443 fileName += (dynamic_cast<TObjString*>(fFileList->Last()))->GetString();
444 AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
448 if (fileName.Contains(".root")) {
449 AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
450 SetPercentileFile(fileName.Data());