AliPhysics  7f2a7c4 (7f2a7c4)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskCRCZDC.cxx
Go to the documentation of this file.
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 thce 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  * analysis task for CRC with ZDC *
18  * *
19  * author: Jacopo Margutti *
20  * (margutti@nikhef.nl) *
21  **********************************/
22 
23 #include "Riostream.h" //needed as include
24 #include "TChain.h"
25 #include "TList.h"
26 #include "TTree.h"
27 #include "TRandom3.h"
28 #include "TTimeStamp.h"
29 #include "TStopwatch.h"
30 #include "TProfile.h"
31 #include <TList.h>
32 #include <TH1F.h>
33 #include <TH2F.h>
34 #include <TH3D.h>
35 #include <TFile.h>
36 #include <TString.h>
37 #include <TCanvas.h>
38 #include <TParticle.h>
39 #include "TProfile2D.h"
40 #include "TProfile3D.h"
41 
42 #include "AliAnalysisManager.h"
43 #include "AliInputEventHandler.h"
44 #include "AliVEvent.h"
45 #include "AliESD.h"
46 #include "AliESDEvent.h"
47 #include "AliESDHeader.h"
48 #include "AliESDInputHandler.h"
49 #include "AliESDZDC.h"
50 #include "AliMultiplicity.h"
51 #include "AliAnalysisUtils.h"
52 #include "AliAODHandler.h"
53 #include "AliAODTrack.h"
54 #include "AliAODEvent.h"
55 #include "AliAODHeader.h"
56 #include "AliAODVertex.h"
57 #include "AliAODVZERO.h"
58 #include "AliAODZDC.h"
59 #include "AliAODMCHeader.h"
60 #include "AliMCEventHandler.h"
61 #include "AliMCEvent.h"
62 #include "AliHeader.h"
63 #include "AliVParticle.h"
64 #include "AliStack.h"
65 #include "AliAODMCParticle.h"
66 #include "AliAnalysisTaskSE.h"
67 #include "AliGenEventHeader.h"
68 #include "AliPhysicsSelectionTask.h"
69 #include "AliPhysicsSelection.h"
70 #include "AliBackgroundSelection.h"
71 #include "AliTriggerAnalysis.h"
72 #include "AliCentrality.h"
73 #include "AliAnalysisTaskCRCZDC.h"
74 #include "AliMultSelection.h"
75 #include "AliLumiTools.h"
76 
77 // ALICE Correction Framework
78 #include "AliCFManager.h"
79 
80 // Interface to Event generators to get Reaction Plane Angle
81 #include "AliGenCocktailEventHeader.h"
82 #include "AliGenPythiaEventHeader.h"
83 #include "AliGenHijingEventHeader.h"
84 #include "AliGenGeVSimEventHeader.h"
85 #include "AliGenEposEventHeader.h"
86 
87 // Interface to Load short life particles
88 #include "TObjArray.h"
89 #include "AliFlowCandidateTrack.h"
90 
91 // Interface to make the Flow Event Simple used in the flow analysis methods
92 #include "AliFlowEvent.h"
93 #include "AliFlowTrackCuts.h"
94 #include "AliFlowEventCuts.h"
95 #include "AliFlowCommonConstants.h"
96 
98 
99 //________________________________________________________________________
102 fAnalysisType("AUTOMATIC"),
103 fRPType(""),
104 fCFManager1(NULL),
105 fCFManager2(NULL),
106 fCutsEvent(NULL),
107 fCutsRP(NULL),
108 fCutsPOI(NULL),
109 fCutContainer(new TList()),
110 fQAList(NULL),
111 fMinMult(0),
112 fMaxMult(10000000),
113 fMinA(-1.0),
114 fMaxA(-0.01),
115 fMinB(0.01),
116 fMaxB(1.0),
117 fGenHeader(NULL),
118 fPythiaGenHeader(NULL),
119 fHijingGenHeader(NULL),
120 fFlowTrack(NULL),
121 fAnalysisUtil(NULL),
122 fQAon(kFALSE),
123 fLoadCandidates(kFALSE),
124 fNbinsMult(10000),
125 fNbinsPt(100),
126 fNbinsPhi(100),
127 fNbinsEta(200),
128 fNbinsQ(500),
129 fNbinsMass(1),
130 fMultMin(0.),
131 fMultMax(10000.),
132 fPtMin(0.),
133 fPtMax(10.),
134 fPhiMin(0.),
135 fPhiMax(TMath::TwoPi()),
136 fEtaMin(-5.),
137 fEtaMax(5.),
138 fQMin(0.),
139 fQMax(3.),
140 fMassMin(-1.),
141 fMassMax(0.),
142 fHistWeightvsPhiMin(0.),
143 fHistWeightvsPhiMax(3.),
144 fExcludedEtaMin(0.),
145 fExcludedEtaMax(0.),
146 fExcludedPhiMin(0.),
147 fExcludedPhiMax(0.),
148 fAfterburnerOn(kFALSE),
149 fNonFlowNumberOfTrackClones(0),
150 fV1(0.),
151 fV2(0.),
152 fV3(0.),
153 fV4(0.),
154 fV5(0.),
155 fDifferentialV2(0),
156 fFlowEvent(NULL),
157 fShuffleTracks(kFALSE),
158 fMyTRandom3(NULL),
159 fAnalysisInput(kAOD),
160 fIsMCInput(kFALSE),
161 fUseMCCen(kTRUE),
162 fRejectPileUp(kTRUE),
163 fRejectPileUpTight(kFALSE),
164 fCentrLowLim(0.),
165 fCentrUpLim(100.),
166 fCentrEstimator(kV0M),
167 fOutput(0x0),
168 fhZNCvsZNA(0x0),
169 fhZDCCvsZDCCA(0x0),
170 fhZNCvsZPC(0x0),
171 fhZNAvsZPA(0x0),
172 fhZNvsZP(0x0),
173 fhZNvsVZERO(0x0),
174 fhZDCvsVZERO(0x0),
175 fhZDCvsTracklets(0x0),
176 fhZDCvsNclu1(0x0),
177 fhDebunch(0x0),
178 fhAsymm(0x0),
179 fhZNAvsAsymm(0x0),
180 fhZNCvsAsymm(0x0),
181 fhZNCvscentrality(0x0),
182 fhZNAvscentrality(0x0),
183 fCRCnRun(0),
184 fZDCGainAlpha(0.395),
185 fDataSet(kAny),
186 fStack(0x0),
187 fCutTPC(kFALSE),
188 fCenDis(0x0),
189 fVZEROMult(0x0),
190 fMultSelection(0x0),
191 fPileUpCount(0x0),
192 fPileUpMultSelCount(0x0),
193 fMultTOFLowCut(0x0),
194 fMultTOFHighCut(0x0),
195 fUseTowerEq(kFALSE),
196 fTowerEqList(NULL),
197 fUseBadTowerCalib(kFALSE),
198 fBadTowerCalibList(NULL),
199 fUseZDCSpectraCorr(kFALSE),
200 fZDCSpectraCorrList(NULL),
201 fSpectraMCList(NULL),
202 fBadTowerStuffList(NULL),
203 fCachedRunNum(0),
204 fhZNSpectra(0x0),
205 fhZNSpectraCor(0x0),
206 fhZNSpectraPow(0x0),
207 fhZNBCCorr(0x0)
208 {
209  for(int i=0; i<5; i++){
210  fhZNCPM[i] = 0x0;
211  fhZNAPM[i] = 0x0;
212  }
213  for(int i=0; i<4; i++){
214  fhZNCPMQiPMC[i] = 0x0;
215  fhZNAPMQiPMC[i] = 0x0;
216  }
217  for(Int_t r=0; r<fCRCMaxnRun; r++) {
218  fRunList[r] = 0;
219  }
220  for(Int_t c=0; c<2; c++) {
221  for(Int_t i=0; i<5; i++) {
222  fTowerGainEq[c][i] = NULL;
223  }
224  }
225  for(Int_t c=0; c<100; c++) {
226  fBadTowerCalibHist[c] = NULL;
227  }
228  for(Int_t i=0; i<8; i++) {
229  SpecCorMu1[i] = NULL;
230  SpecCorMu2[i] = NULL;
231  SpecCorSi[i] = NULL;
232  SpecCorAv[i] = NULL;
233  }
234  this->InitializeRunArrays();
235  fMyTRandom3 = new TRandom3(1);
236  gRandom->SetSeed(fMyTRandom3->Integer(65539));
237  for(Int_t j=0; j<2; j++) {
238  for(Int_t c=0; c<10; c++) {
239  fPtSpecGen[j][c] = NULL;
240  fPtSpecFB32[j][c] = NULL;
241  fPtSpecFB96[j][c] = NULL;
242  fPtSpecFB128[j][c] = NULL;
243  fPtSpecFB768[j][c] = NULL;
244  }
245  }
246  for (Int_t c=0; c<2; c++) {
247  fhZNCenDis[c] = NULL;
248  }
249 }
250 
251 //________________________________________________________________________
252 AliAnalysisTaskCRCZDC::AliAnalysisTaskCRCZDC(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates):
253 AliAnalysisTaskSE(name),
254 fAnalysisType("AUTOMATIC"),
255 fRPType(RPtype),
256 fCFManager1(NULL),
257 fCFManager2(NULL),
258 fCutsEvent(NULL),
259 fCutsRP(NULL),
260 fCutsPOI(NULL),
261 fCutContainer(new TList()),
262 fQAList(NULL),
263 fMinMult(0),
264 fMaxMult(10000000),
265 fMinA(-1.0),
266 fMaxA(-0.01),
267 fMinB(0.01),
268 fMaxB(1.0),
269 fQAon(on),
270 fLoadCandidates(bCandidates),
271 fNbinsMult(10000),
272 fNbinsPt(100),
273 fNbinsPhi(100),
274 fNbinsEta(200),
275 fNbinsQ(500),
276 fNbinsMass(1),
277 fMultMin(0.),
278 fMultMax(10000.),
279 fPtMin(0.),
280 fPtMax(10.),
281 fPhiMin(0.),
282 fPhiMax(TMath::TwoPi()),
283 fEtaMin(-5.),
284 fEtaMax(5.),
285 fQMin(0.),
286 fQMax(3.),
287 fMassMin(-1.),
288 fMassMax(0.),
289 fHistWeightvsPhiMin(0.),
290 fHistWeightvsPhiMax(3.),
291 fExcludedEtaMin(0.),
292 fExcludedEtaMax(0.),
293 fExcludedPhiMin(0.),
294 fExcludedPhiMax(0.),
295 fAfterburnerOn(kFALSE),
296 fNonFlowNumberOfTrackClones(0),
297 fV1(0.),
298 fV2(0.),
299 fV3(0.),
300 fV4(0.),
301 fV5(0.),
302 fDifferentialV2(0),
303 fFlowEvent(NULL),
304 fShuffleTracks(kFALSE),
305 fMyTRandom3(NULL),
306 fAnalysisInput(kAOD),
307 fIsMCInput(kFALSE),
308 fUseMCCen(kTRUE),
309 fRejectPileUp(kTRUE),
310 fRejectPileUpTight(kFALSE),
311 fCentrLowLim(0.),
312 fCentrUpLim(100.),
313 fCentrEstimator(kV0M),
314 fOutput(0x0),
315 fhZNCvsZNA(0x0),
316 fhZDCCvsZDCCA(0x0),
317 fhZNCvsZPC(0x0),
318 fhZNAvsZPA(0x0),
319 fhZNvsZP(0x0),
320 fhZNvsVZERO(0x0),
321 fhZDCvsVZERO(0x0),
322 fhZDCvsTracklets(0x0),
323 fhZDCvsNclu1(0x0),
324 fhDebunch(0x0),
325 fhAsymm(0x0),
326 fhZNAvsAsymm(0x0),
327 fhZNCvsAsymm(0x0),
328 fhZNCvscentrality(0x0),
329 fhZNAvscentrality(0x0),
330 fDataSet(kAny),
331 fCRCnRun(0),
332 fZDCGainAlpha(0.395),
333 fGenHeader(NULL),
334 fPythiaGenHeader(NULL),
335 fHijingGenHeader(NULL),
336 fFlowTrack(NULL),
337 fAnalysisUtil(NULL),
338 fStack(0x0),
339 fCutTPC(kFALSE),
340 fCenDis(0x0),
341 fVZEROMult(0x0),
342 fMultSelection(0x0),
343 fPileUpCount(0x0),
344 fPileUpMultSelCount(0x0),
345 fMultTOFLowCut(0x0),
346 fMultTOFHighCut(0x0),
347 fUseTowerEq(kFALSE),
348 fTowerEqList(NULL),
349 fUseBadTowerCalib(kFALSE),
350 fBadTowerCalibList(NULL),
351 fUseZDCSpectraCorr(kFALSE),
352 fZDCSpectraCorrList(NULL),
353 fSpectraMCList(NULL),
354 fBadTowerStuffList(NULL),
355 fCachedRunNum(0),
356 fhZNSpectra(0x0),
357 fhZNSpectraCor(0x0),
358 fhZNSpectraPow(0x0),
359 fhZNBCCorr(0x0)
360 {
361 
362  for(int i=0; i<5; i++){
363  fhZNCPM[i] = 0x0;
364  fhZNAPM[i] = 0x0;
365  }
366  for(int i=0; i<4; i++){
367  fhZNCPMQiPMC[i] = 0x0;
368  fhZNAPMQiPMC[i] = 0x0;
369  }
370  for(Int_t r=0; r<fCRCMaxnRun; r++) {
371  fRunList[r] = 0;
372  }
373  for(Int_t c=0; c<2; c++) {
374  for(Int_t i=0; i<5; i++) {
375  fTowerGainEq[c][i] = NULL;
376  }
377  }
378  for(Int_t c=0; c<100; c++) {
379  fBadTowerCalibHist[c] = NULL;
380  }
381  for(Int_t i=0; i<8; i++) {
382  SpecCorMu1[i] = NULL;
383  SpecCorMu2[i] = NULL;
384  SpecCorSi[i] = NULL;
385  SpecCorAv[i] = NULL;
386  }
387  this->InitializeRunArrays();
388  fMyTRandom3 = new TRandom3(iseed);
389  gRandom->SetSeed(fMyTRandom3->Integer(65539));
390 
391  DefineInput(0, TChain::Class());
392  // Define output slots here
393  // Define here the flow event output
394  DefineOutput(1, AliFlowEventSimple::Class());
395  DefineOutput(2, TList::Class());
396 
397  for(Int_t j=0; j<2; j++) {
398  for(Int_t c=0; c<10; c++) {
399  fPtSpecGen[j][c] = NULL;
400  fPtSpecFB32[j][c] = NULL;
401  fPtSpecFB96[j][c] = NULL;
402  fPtSpecFB128[j][c] = NULL;
403  fPtSpecFB768[j][c] = NULL;
404  }
405  }
406  for (Int_t c=0; c<2; c++) {
407  fhZNCenDis[c] = NULL;
408  }
409 }
410 
411 //________________________________________________________________________
413 {
414  // Destructor
415  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
416  delete fOutput; fOutput=0;
417  }
418  delete fMyTRandom3;
419  delete fFlowEvent;
420  delete fFlowTrack;
421  delete fCutsEvent;
422  if (fTowerEqList) delete fTowerEqList;
425  if (fAnalysisUtil) delete fAnalysisUtil;
426  if (fQAList) delete fQAList;
427  if (fCutContainer) fCutContainer->Delete(); delete fCutContainer;
428 }
429 
430 //________________________________________________________________________
432 {
433  for(Int_t r=0;r<fCRCMaxnRun;r++) {
434  fCRCQVecListRun[r] = NULL;
435  for(Int_t k=0;k<fCRCnTow;k++) {
436  fZNCTower[r][k] = NULL;
437  fZNATower[r][k] = NULL;
438  }
439 // fhZNSpectraRbR[r] = NULL;
440  }
441  // for(Int_t i=0;i<fnCen;i++) {
442  // fPtPhiEtaRbRFB128[r][i] = NULL;
443  // fPtPhiEtaRbRFB768[r][i] = NULL;
444  // }
445  // }
446 }
447 
448 //________________________________________________________________________
450 {
451  // Create the output containers
452 
453  if (!(fAnalysisType == "AOD" || fAnalysisType == "MCkine" || fAnalysisType == "MCAOD" || fAnalysisType == "AUTOMATIC" || fAnalysisType == "MCESD"))
454  {
455  AliError("WRONG ANALYSIS TYPE! only MCESD, MCkine, MCAOD, AOD and AUTOMATIC are allowed.");
456  exit(1);
457  }
458 
459  //set the common constants
462  cc->SetNbinsPt(fNbinsPt);
463  cc->SetNbinsPhi(fNbinsPhi);
464  cc->SetNbinsEta(fNbinsEta);
465  cc->SetNbinsQ(fNbinsQ);
467  cc->SetMultMin(fMultMin);
468  cc->SetMultMax(fMultMax);
469  cc->SetPtMin(fPtMin);
470  cc->SetPtMax(fPtMax);
471  cc->SetPhiMin(fPhiMin);
472  cc->SetPhiMax(fPhiMax);
473  cc->SetEtaMin(fEtaMin);
474  cc->SetEtaMax(fEtaMax);
475  cc->SetQMin(fQMin);
476  cc->SetQMax(fQMax);
477  cc->SetMassMin(fMassMin);
478  cc->SetMassMax(fMassMax);
481 
482  fFlowEvent = new AliFlowEvent(20000);
483  fFlowTrack = new AliFlowTrack();
484 
485  //printf(" AliAnalysisTaskCRCZDC::UserCreateOutputObjects()\n\n");
486  fOutput = new TList();
487  fOutput->SetOwner(kTRUE);
488  //fOutput->SetName("output");
489 
490  if (fQAon) {
491  fQAList = new TList();
492  fQAList->SetOwner(kTRUE);
493  fQAList->SetName("AliFlowEventCuts QA");
494  if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA()); //0
495  if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA()); //1
496  if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA()); //2
497  fOutput->Add(fQAList);
498  }
499 
500  fCenDis = new TH1F("fCenDis", "fCenDis", 100, 0., 100.);
501  fOutput->Add(fCenDis);
502  fPileUpCount = new TH1F("fPileUpCount", "fPileUpCount", 9, 0., 9.);
503  fPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
504  fPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
505  fPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultiplicityComb08");
506  fPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
507  fPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>7.5");
508  fPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
509  fPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
510  fPileUpCount->GetXaxis()->SetBinLabel(8,"multESDTPCDif");
511  fPileUpCount->GetXaxis()->SetBinLabel(9,"extraPileUpMultSel");
512  fOutput->Add(fPileUpCount);
513  fPileUpMultSelCount = new TH1F("fPileUpMultSelCount", "fPileUpMultSelCount", 8, 0., 8.);
514  fPileUpMultSelCount->GetXaxis()->SetBinLabel(1,"IsNotPileup");
515  fPileUpMultSelCount->GetXaxis()->SetBinLabel(2,"IsNotPileupMV");
516  fPileUpMultSelCount->GetXaxis()->SetBinLabel(3,"IsNotPileupInMultBins");
517  fPileUpMultSelCount->GetXaxis()->SetBinLabel(4,"InconsistentVertices");
518  fPileUpMultSelCount->GetXaxis()->SetBinLabel(5,"TrackletVsCluster");
519  fPileUpMultSelCount->GetXaxis()->SetBinLabel(6,"AsymmetricInVZERO");
520  fPileUpMultSelCount->GetXaxis()->SetBinLabel(7,"IncompleteDAQ");
521  fPileUpMultSelCount->GetXaxis()->SetBinLabel(8,"GoodVertex2016");
523 
524  fMultTOFLowCut = new TF1("fMultTOFLowCut", "[0]+[1]*x+[2]*x*x+[3]*x*x*x - 4.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x+[9]*x*x*x*x*x)", 0, 10000);
525  fMultTOFLowCut->SetParameters(-1.0178, 0.333132, 9.10282e-05, -1.61861e-08, 1.47848, 0.0385923, -5.06153e-05, 4.37641e-08, -1.69082e-11, 2.35085e-15);
526  fOutput->Add(fMultTOFLowCut);
527  fMultTOFHighCut = new TF1("fMultTOFHighCut", "[0]+[1]*x+[2]*x*x+[3]*x*x*x + 4.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x+[9]*x*x*x*x*x)", 0, 10000);
528  fMultTOFHighCut->SetParameters(-1.0178, 0.333132, 9.10282e-05, -1.61861e-08, 1.47848, 0.0385923, -5.06153e-05, 4.37641e-08, -1.69082e-11, 2.35085e-15);
529  fOutput->Add(fMultTOFHighCut);
530 
531  for(Int_t c=0; c<2; c++) {
532  for(Int_t i=0; i<5; i++) {
533  fTowerGainEq[c][i] = new TH1D();
534  fOutput->Add(fTowerGainEq[c][i]);
535  }
536  }
537 
538  for (Int_t c=0; c<2; c++) {
539  fhZNCenDis[c] = new TH3D(Form("fhZNCenDis[%d]",c), Form("fhZNCenDis[%d]",c), 100, 0., 100., 100, -2., 2. , 100., -2., 2.);
540  fOutput->Add(fhZNCenDis[c]);
541  }
542 
543  fBadTowerStuffList = new TList();
544  fBadTowerStuffList->SetOwner(kTRUE);
545  fBadTowerStuffList->SetName("BadTowerCalib");
547 
548  if(fBadTowerCalibList) {
549  for(Int_t c=0; c<100; c++) {
550  fBadTowerCalibHist[c] = (TH2D*)fBadTowerCalibList->FindObject(Form("TH2Resp[%d]",c));
552  }
553  }
554  if(fZDCSpectraCorrList) {
555  for(Int_t i=0; i<8; i++) {
556  SpecCorMu1[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu1[%d]",i));
557  fOutput->Add(SpecCorMu1[i]);
558  SpecCorMu2[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu2[%d]",i));
559  fOutput->Add(SpecCorMu2[i]);
560  SpecCorAv[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorAv[%d]",i));
561  fOutput->Add(SpecCorAv[i]);
562  SpecCorSi[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorSi[%d]",i));
563  fOutput->Add(SpecCorSi[i]);
564  }
565  }
566 
567  fhZNSpectra = new TH3D("fhZNSpectra","fhZNSpectra",100,0.,100.,8,0.,8.,1000,0.,1.E5);
568  fOutput->Add(fhZNSpectra);
569  fhZNSpectraCor = new TH3D("fhZNSpectraCor","fhZNSpectraCor",100,0.,100.,8,0.,8.,1000,0.,1.E5);
570  fOutput->Add(fhZNSpectraCor);
571  fhZNSpectraPow = new TH3D("fhZNSpectraPow","fhZNSpectraPow",100,0.,100.,8,0.,8.,1000,0.,TMath::Power(1.E5,fZDCGainAlpha));
572  fOutput->Add(fhZNSpectraPow);
573  fhZNBCCorr = new TH3D("fhZNBCCorr","fhZNBCCorr",100,0.,100.,500,0.,1.E5,500,0.,1.E5);
574  fOutput->Add(fhZNBCCorr);
575 
576  if(fAnalysisType == "MCAOD") {
577 
578  fSpectraMCList = new TList();
579  fSpectraMCList->SetOwner(kTRUE);
580  fSpectraMCList->SetName("Spectra");
581  fOutput->Add(fSpectraMCList);
582 
583  Double_t xmin[] = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.,1.2,1.4,1.6,1.8,2.,2.33,2.66,3.,3.5,4.,5.,6.,9.,20.};
584  for(Int_t j=0; j<2; j++) {
585  for(Int_t c=0; c<10; c++) {
586  fPtSpecGen[j][c] = new TH1F(Form("fPtSpecGen[%d][%d]",j,c), Form("fPtSpecGen[%d][%d]",j,c), 23, xmin);
587  fSpectraMCList->Add(fPtSpecGen[j][c]);
588  fPtSpecFB32[j][c] = new TH1F(Form("fPtSpecFB32[%d][%d]",j,c), Form("fPtSpecFB32[%d][%d]",j,c), 23, xmin);
589  fSpectraMCList->Add(fPtSpecFB32[j][c]);
590  fPtSpecFB96[j][c] = new TH1F(Form("fPtSpecFB96[%d][%d]",j,c), Form("fPtSpecFB96[%d][%d]",j,c), 23, xmin);
591  fSpectraMCList->Add(fPtSpecFB96[j][c]);
592  fPtSpecFB128[j][c] = new TH1F(Form("fPtSpecFB128[%d][%d]",j,c), Form("fPtSpecFB128[%d][%d]",j,c), 23, xmin);
593  fSpectraMCList->Add(fPtSpecFB128[j][c]);
594  fPtSpecFB768[j][c] = new TH1F(Form("fPtSpecFB768[%d][%d]",j,c), Form("fPtSpecFB768[%d][%d]",j,c), 23, xmin);
595  fSpectraMCList->Add(fPtSpecFB768[j][c]);
596  }
597  }
598  }
599 
600  fAnalysisUtil = new AliAnalysisUtils();
601  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
602  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
603 
604  for(int i=0; i<5; i++){
605  char hname[20];
606  sprintf(hname,"hZNCPM%d",i);
607  fhZNCPM[i] = new TH1F(hname, hname, 200, -50., 140000);
608  fOutput->Add(fhZNCPM[i]);
609  //
610  sprintf(hname,"hZNAPM%d",i);
611  fhZNAPM[i] = new TH1F(hname, hname, 200, -50., 140000);
612  fOutput->Add(fhZNAPM[i]);
613  //
614  if(i<4){
615  //
616  char hnamenc[20];
617  sprintf(hnamenc, "hZNCPMQ%dPMC",i+1);
618  fhZNCPMQiPMC[i] = new TH1F(hnamenc, hnamenc, 100, 0., 1.);
619  fOutput->Add(fhZNCPMQiPMC[i]);
620  //
621  char hnamena[20];
622  sprintf(hnamena, "hZNAPMQ%dPMC",i+1);
623  fhZNAPMQiPMC[i] = new TH1F(hnamena, hnamena, 100, 0., 1.);
624  fOutput->Add(fhZNAPMQiPMC[i]);
625  }
626  }
627 
628  fhZNCvsZNA = new TH2F("hZNCvsZNA","hZNCvsZNA",200,-50.,140000,200,-50.,140000);
629  fOutput->Add(fhZNCvsZNA);
630  fhZDCCvsZDCCA = new TH2F("hZDCCvsZDCCA","hZDCCvsZDCCA",200,0.,180000.,200,0.,200000.);
631  fOutput->Add(fhZDCCvsZDCCA);
632  fhZNCvsZPC = new TH2F("hZNCvsZPC","hZNCvsZPC",200,-50.,50000,200,-50.,140000);
633  fOutput->Add(fhZNCvsZPC);
634  fhZNAvsZPA = new TH2F("hZNAvsZPA","hZNAvsZPA",200,-50.,50000,200,-50.,140000);
635  fOutput->Add(fhZNAvsZPA);
636  fhZNvsZP = new TH2F("hZNvsZP","hZNvsZP",200,-50.,80000,200,-50.,200000);
637  fOutput->Add(fhZNvsZP);
638  fhZNvsVZERO = new TH2F("hZNvsVZERO","hZNvsVZERO",250,0.,25000.,200,0.,200000.);
639  fOutput->Add(fhZNvsVZERO);
640  fhZDCvsVZERO = new TH2F("hZDCvsVZERO","hZDCvsVZERO",250,0.,25000.,250,0.,250000.);
641  fOutput->Add(fhZDCvsVZERO);
642  fhZDCvsTracklets = new TH2F("hZDCvsTracklets","hZDCvsTracklets",200,0.,4000.,250,0.,250000.);
644  fhZDCvsNclu1 = new TH2F("hZDCvsNclu1", "hZDCvsNclu1", 200, 0.,8000.,200,0.,250000.);
645  fOutput->Add(fhZDCvsNclu1);
646  fhDebunch = new TH2F("hDebunch","hDebunch",240,-100.,-40.,240,-30.,30.);
647  fOutput->Add(fhDebunch);
648 
649  fhAsymm = new TH1F("hAsymm" , "Asimmetry ",200,-1.,1.);
650  fOutput->Add(fhAsymm);
651  fhZNAvsAsymm = new TH2F("hZNAvsAsymm","ZNA vs. asymm.",200,-1.,1.,200,0.,80.);
652  fOutput->Add(fhZNAvsAsymm);
653  fhZNCvsAsymm = new TH2F("hZNCvsAsymm","ZNC vs. asymm.",200,-1.,1.,200,0.,80.);
654  fOutput->Add(fhZNCvsAsymm);
655 
656  fhZNCvscentrality = new TH2F("hZNCvscentrality","ZNC vs. centrality",100,0.,100.,100,0.,100.);
658  fhZNAvscentrality = new TH2F("hZNAvscentrality","ZNA vs. centrality",100,0.,100.,100,0.,100.);
660 
661  //********************************************************************
662 
663  Int_t dRun10h[] = {139510, 139507, 139505, 139503, 139465, 139438, 139437, 139360, 139329, 139328, 139314, 139310, 139309, 139173, 139107, 139105, 139038, 139037, 139036, 139029, 139028, 138872, 138871, 138870, 138837, 138732, 138730, 138666, 138662, 138653, 138652, 138638, 138624, 138621, 138583, 138582, 138579, 138578, 138534, 138469, 138442, 138439, 138438, 138396, 138364, 138275, 138225, 138201, 138197, 138192, 138190, 137848, 137844, 137752, 137751, 137724, 137722, 137718, 137704, 137693, 137692, 137691, 137686, 137685, 137639, 137638, 137608, 137595, 137549, 137546, 137544, 137541, 137539, 137531, 137530, 137443, 137441, 137440, 137439, 137434, 137432, 137431, 137430, 137366, 137243, 137236, 137235, 137232, 137231, 137230, 137162, 137161};
664 
665  Int_t dRun11h[] = {167902, 167903, 167915, 167920, 167985, 167987, 167988, 168066, 168068, 168069, 168076, 168104, 168105, 168107, 168108, 168115, 168212, 168310, 168311, 168322, 168325, 168341, 168342, 168361, 168362, 168458, 168460, 168461, 168464, 168467, 168511, 168512, 168514, 168777, 168826, 168984, 168988, 168992, 169035, 169040, 169044, 169045, 169091, 169094, 169099, 169138, 169143, 169144, 169145, 169148, 169156, 169160, 169167, 169238, 169411, 169415, 169417, 169418, 169419, 169420, 169475, 169498, 169504, 169506, 169512, 169515, 169550, 169553, 169554, 169555, 169557, 169586, 169587, 169588, 169590, 169591, 169835, 169837, 169838, 169846, 169855, 169858, 169859, 169923, 169956, 169965, 170027, 170036,170040, 170081, 170083, 170084, 170085, 170088, 170089, 170091, 170155, 170159, 170163, 170193, 170203, 170204, 170207, 170228, 170230, 170268, 170269, 170270, 170306, 170308, 170309, 170311, 170312, 170313, 170315, 170387, 170388, 170572, 170593};
666 
667  // 12 low IR: 244917, 244918, 244975, 244980, 244982, 244983, 245064, 245066, 245068, 246390, 246391, 246392
668  // 80 high IR ("CentralBarrelTracking" good runs): 246994, 246991, 246989, 246984, 246982, 246980, 246948, 246945, 246928, 246871, 246870, 246867, 246865, 246864, 246859, 246858, 246851, 246847, 246846, 246845, 246844, 246810, 246809, 246808, 246807, 246805, 246804, 246766, 246765, 246763, 246760, 246759, 246758, 246757, 246751, 246750, 246676, 246675, 246540, 246495, 246493, 246488, 246487, 246434, 246431, 246428, 246424, 246276, 246275, 246272, 246271, 246225, 246222, 246217, 246185, 246182, 246181, 246180, 246178, 246153, 246152, 246151, 246115, 246113, 246089, 246087, 246053, 246052, 246049, 246048, 246042, 246037, 246036, 246012, 246003, 246001, 245954, 245952, 245949, 245923, 245833, 245831, 245829, 245705, 245702, 245700, 245692, 245683
669 
670  Int_t dRun15o[] = {244917, 244918, 244975, 244980, 244982, 244983, 245064, 245066, 245068, 246390, 246391, 246392, 246994, 246991, 246989, 246984, 246982, 246980, 246948, 246945, 246928, 246851, 246847, 246846, 246845, 246844, 246810, 246809, 246808, 246807, 246805, 246804, 246766, 246765, 246763, 246760, 246759, 246758, 246757, 246751, 246750, 246495, 246493, 246488, 246487, 246434, 246431, 246428, 246424, 246276, 246275, 246272, 246271, 246225, 246222, 246217, 246185, 246182, 246181, 246180, 246178, 246153, 246152, 246151, 246115, 246113, 246089, 246087, 246053, 246052, 246049, 246048, 246042, 246037, 246036, 246012, 246003, 246001, 245954, 245952, 245949, 245923, 245833, 245831, 245829, 245705, 245702, 245700, 245692, 245683};
671 
672  Int_t dRun15ov6[] = {244918, 244975, 244980, 244982, 244983, 245064, 245066, 245068, 246390, 246391, 246392, 246994, 246991, 246989, 246984, 246982, 246980, 246948, 246945, 246928, 246851, 246847, 246846, 246845, 246844, 246810, 246809, 246808, 246807, 246805, 246804, 246766, 246765, 246763, 246760, 246759, 246758, 246757, 246751, 246750, 246495, 246493, 246488, 246487, 246434, 246431, 246428, 246424, 246276, 246275, 246272, 246271, 246225, 246222, 246217, 246185, 246182, 246181, 246180, 246178, 246153, 246152, 246151, 246148, 246115, 246113, 246089, 246087, 246053, 246052, 246049, 246048, 246042, 246037, 246036, 246012, 246003, 246001, 245963, 245954, 245952, 245949, 245923, 245833, 245831, 245829, 245705, 245702, 245700, 245692, 245683};
673 
674  if(fDataSet==k2010) {fCRCnRun=92;}
675  if(fDataSet==k2011) {fCRCnRun=119;}
676  if(fDataSet==k2015) {fCRCnRun=90;}
677  if(fDataSet==k2015v6) {fCRCnRun=91;}
678  if(fDataSet==kAny) {fCRCnRun=1;}
679 
680  Int_t d=0;
681  for(Int_t r=0; r<fCRCnRun; r++) {
682  if(fDataSet==k2010) fRunList[d] = dRun10h[r];
683  if(fDataSet==k2011) fRunList[d] = dRun11h[r];
684  if(fDataSet==k2015) fRunList[d] = dRun15o[r];
685  if(fDataSet==k2015v6) fRunList[d] = dRun15ov6[r];
686  if(fDataSet==kAny) fRunList[d] = 1;
687  d++;
688  }
689 
690  fVZEROMult = new TProfile2D("fVZEROMult","fVZEROMult",fCRCnRun,0.,1.*fCRCnRun,64,0.,64.);
691  for (Int_t i=0; i<fCRCnRun; i++) {
692  fVZEROMult->GetXaxis()->SetBinLabel(i+1,Form("%d",fRunList[i]));
693  }
694  fOutput->Add(fVZEROMult);
695 
696  Double_t ptmin[] = {0.2,0.4,0.6,0.8,1.,1.2,1.4,1.8,2.2,3.,4.,6.,8.,12.,20.};
697  Double_t phimin[] = {0.,TMath::Pi()/8.,2*TMath::Pi()/8.,3*TMath::Pi()/8.,4*TMath::Pi()/8.,5*TMath::Pi()/8.,6*TMath::Pi()/8.,7*TMath::Pi()/8.,8*TMath::Pi()/8.,9*TMath::Pi()/8.,10*TMath::Pi()/8.,11*TMath::Pi()/8.,12*TMath::Pi()/8.,13*TMath::Pi()/8.,14*TMath::Pi()/8.,15*TMath::Pi()/8.,16*TMath::Pi()/8.};
698  Double_t etamin[] = {-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8};
699 
700  for(Int_t r=0;r<fCRCnRun;r++) {
701  fCRCQVecListRun[r] = new TList();
702  fCRCQVecListRun[r]->SetName(Form("Run %d",fRunList[r]));
703  fCRCQVecListRun[r]->SetOwner(kTRUE);
704  fOutput->Add(fCRCQVecListRun[r]);
705 
706  if(!fUseTowerEq) {
707  for(Int_t k=0;k<fCRCnTow;k++) {
708  fZNCTower[r][k] = new TProfile(Form("fZNCTower[%d][%d]",fRunList[r],k),Form("fZNCTower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
709  fZNCTower[r][k]->Sumw2();
710  fCRCQVecListRun[r]->Add(fZNCTower[r][k]);
711  fZNATower[r][k] = new TProfile(Form("fZNATower[%d][%d]",fRunList[r],k),Form("fZNATower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
712  fZNATower[r][k]->Sumw2();
713  fCRCQVecListRun[r]->Add(fZNATower[r][k]);
714  }
715  }
716 
717 // fhZNSpectraRbR[r] = new TH3D(Form("fhZNSpectraRbR[%d]",fRunList[r]),Form("fhZNSpectraRbR[%d]",fRunList[r]),50,0.,100.,8,0.,8.,100,0.,1.E5);
718 // fCRCQVecListRun[r]->Add(fhZNSpectraRbR[r]);
719 
720  // for(Int_t i=0;i<fnCen;i++) {
721  // fPtPhiEtaRbRFB128[r][i] = new TH3F(Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
722  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB128[r][i]);
723  // fPtPhiEtaRbRFB768[r][i] = new TH3F(Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
724  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB768[r][i]);
725  // }
726  }
727 
728  PostData(2, fOutput);
729 }
730 
731 //________________________________________________________________________
733 {
734  // Execute analysis for current event:
735  AliMCEvent* McEvent = MCEvent();
736  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
737  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
738  // AliMultiplicity* myTracklets = NULL;
739  // AliESDPmdTrack* pmdtracks = NULL;
740  // int availableINslot=1;
741 
742  if (!(fCutsRP&&fCutsPOI&&fCutsEvent)) {
743  AliError("cuts not set");
744  return;
745  }
746 
747  Int_t RunBin=-1, bin=0;
748  Int_t RunNum = aod->GetRunNumber();
749  for(Int_t c=0;c<fCRCnRun;c++) {
750  if(fRunList[c]==RunNum) RunBin=bin;
751  else bin++;
752  }
753  if(RunBin==-1) return;
754  if(fDataSet==kAny) RunBin=0;
755 
756  //DEFAULT - automatically takes care of everything
757  if (fAnalysisType == "AUTOMATIC") {
758 
759  // get centrality
760  Double_t centrV0M=300, centrCL1=300, centrCL0=300, centrTRK=300;
761  if(fDataSet!=k2015 && fDataSet!=k2015v6) {
762  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
763  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
764  centrCL0 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL0");
765  centrTRK = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("TRK");
766  } else {
767  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
768  if(!fMultSelection) {
769  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
770  AliWarning("AliMultSelection object not found!");
771  } else {
772  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
773  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
774  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
775  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
776  }
777  }
778 
779  //check event cuts
780  if (InputEvent()) {
781  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
782  if(fRejectPileUp) {
783  if(fDataSet!=k2015 && fDataSet!=k2015v6) {
784 
785  Bool_t BisPileup=kFALSE;
786 
787  // check anyway pileup
788  if (plpMV(aod)) {
789  fPileUpCount->Fill(0.5);
790  BisPileup=kTRUE;
791  }
792 
793  Short_t isPileup = aod->IsPileupFromSPD(3);
794  if (isPileup != 0) {
795  fPileUpCount->Fill(1.5);
796  BisPileup=kTRUE;
797  }
798 
799  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
800  fPileUpCount->Fill(2.5);
801  BisPileup=kTRUE;
802  }
803 
804  if (aod->IsIncompleteDAQ()) {
805  fPileUpCount->Fill(3.5);
806  BisPileup=kTRUE;
807  }
808 
809  // check vertex consistency
810  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
811  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
812 
813  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
814  fPileUpCount->Fill(5.5);
815  BisPileup=kTRUE;
816  }
817 
818  double covTrc[6], covSPD[6];
819  vtTrc->GetCovarianceMatrix(covTrc);
820  vtSPD->GetCovarianceMatrix(covSPD);
821 
822  double dz = vtTrc->GetZ() - vtSPD->GetZ();
823 
824  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
825  double errTrc = TMath::Sqrt(covTrc[5]);
826  double nsigTot = dz/errTot;
827  double nsigTrc = dz/errTrc;
828 
829  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
830  fPileUpCount->Fill(6.5);
831  BisPileup=kTRUE;
832  }
833 
834  if (fAnalysisUtil->IsPileUpEvent(InputEvent())) {
835  fPileUpCount->Fill(7.5);
836  BisPileup=kTRUE;
837  }
838 
839 // if(BisPileup) return;
840  } else {
841  // pileup from AliMultSelection
842  if(!fMultSelection->GetThisEventIsNotPileup()) fPileUpMultSelCount->Fill(0.5);
843  if(!fMultSelection->GetThisEventIsNotPileupMV()) fPileUpMultSelCount->Fill(1.5);
844  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) fPileUpMultSelCount->Fill(2.5);
845  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) fPileUpMultSelCount->Fill(3.5);
846  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) fPileUpMultSelCount->Fill(4.5);
847  if(!fMultSelection->GetThisEventIsNotAsymmetricInVZERO()) fPileUpMultSelCount->Fill(5.5);
848  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) fPileUpMultSelCount->Fill(6.5);
849  if(!fMultSelection->GetThisEventHasGoodVertex2016()) fPileUpMultSelCount->Fill(7.5);
850 
851  Bool_t BisPileup=kFALSE;
852 
853  // pile-up a la Dobrin for LHC15o
854  if (plpMV(aod)) {
855  fPileUpCount->Fill(0.5);
856  BisPileup=kTRUE;
857  }
858 
859  Short_t isPileup = aod->IsPileupFromSPD(3);
860  if (isPileup != 0) {
861  fPileUpCount->Fill(1.5);
862  BisPileup=kTRUE;
863  }
864 
865  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
866  fPileUpCount->Fill(2.5);
867  BisPileup=kTRUE;
868  }
869 
870  if (aod->IsIncompleteDAQ()) {
871  fPileUpCount->Fill(3.5);
872  BisPileup=kTRUE;
873  }
874 
875  if(fabs(centrV0M-centrCL1)>7.5) {
876  fPileUpCount->Fill(4.5);
877  BisPileup=kTRUE;
878  }
879 
880  // check vertex consistency
881  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
882  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
883 
884  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
885  fPileUpCount->Fill(5.5);
886  BisPileup=kTRUE;
887  }
888 
889  double covTrc[6], covSPD[6];
890  vtTrc->GetCovarianceMatrix(covTrc);
891  vtSPD->GetCovarianceMatrix(covSPD);
892 
893  double dz = vtTrc->GetZ() - vtSPD->GetZ();
894 
895  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
896  double errTrc = TMath::Sqrt(covTrc[5]);
897  double nsigTot = dz/errTot;
898  double nsigTrc = dz/errTrc;
899 
900  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
901  fPileUpCount->Fill(6.5);
902  BisPileup=kTRUE;
903  }
904 
905  // cuts on tracks
906  const Int_t nTracks = aod->GetNumberOfTracks();
907  Int_t multEsd = ((AliAODHeader*)aod->GetHeader())->GetNumberOfESDTracks();
908 
909  Int_t multTrk = 0;
910  Int_t multTrkBefC = 0;
911  Int_t multTrkTOFBefC = 0;
912  Int_t multTPC = 0;
913 
914  for (Int_t it = 0; it < nTracks; it++) {
915 
916  AliAODTrack* aodTrk = (AliAODTrack*)aod->GetTrack(it);
917  if (!aodTrk){
918  delete aodTrk;
919  continue;
920  }
921 
922 // if (aodTrk->TestFilterBit(32)){
923 // multTrkBefC++;
924 //
925 // if ( TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
926 // multTrkTOFBefC++;
927 //
928 // if ((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
929 // multTrk++;
930 // }
931 
932  if (aodTrk->TestFilterBit(128))
933  multTPC++;
934 
935  } // end of for (Int_t it = 0; it < nTracks; it++)
936 
937  Double_t multTPCn = multTPC;
938  Double_t multEsdn = multEsd;
939  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
940 
941  if (multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
942  fPileUpCount->Fill(7.5);
943  BisPileup=kTRUE;
944  }
945 
946  if(fRejectPileUpTight) {
947  if(BisPileup==kFALSE) {
948  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
949  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
950  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
951  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
952  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
953  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
954  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
955  if(BisPileup) fPileUpCount->Fill(8.5);
956  }
957  }
958 
959  if(BisPileup) return;
960  }
961  }
962  }
963 
964  //first attach all possible information to the cuts
965  fCutsRP->SetEvent( InputEvent(), MCEvent() ); //attach event
966  fCutsPOI->SetEvent( InputEvent(), MCEvent() );
967 
968  //then make the event
970 
972 
977  fFlowEvent->SetCentralityCL1(centrCL1);
978  fFlowEvent->SetCentralityTRK(centrTRK);
979  // fFlowEvent->SetNITSCL1(((AliVAODHeader*)aod->GetHeader())->GetNumberOfITSClusters(1));
980  Double_t SumV0=0.;
981  for(Int_t i=0; i<64; i++) {
982  if(std::isfinite(aod->GetVZEROEqMultiplicity(i))) SumV0 += aod->GetVZEROEqMultiplicity(i);
983  }
984  fFlowEvent->SetNITSCL1(SumV0);
985 
986  Double_t vtxpos[3]={0.,0.,0.};
987  vtxpos[0] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetX();
988  vtxpos[1] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetY();
989  vtxpos[2] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetZ();
990  fFlowEvent->SetVertexPosition(vtxpos);
991 
992  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
993 
994  // run-by-run QA
995  // for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
996  // AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
997  // if(!track) continue;
998  // // general kinematic & quality cuts
999  // if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70) continue;
1000  // if (track->TestFilterBit(128)) fPtPhiEtaRbRFB128[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1001  // if (track->TestFilterBit(768)) fPtPhiEtaRbRFB768[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1002  // }
1003  fCenDis->Fill(centrV0M);
1004 
1005  }
1006 
1007  if (fAnalysisType == "MCAOD") {
1008 
1009  //check event cuts
1010  if (InputEvent()) {
1011  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
1012  if(fRejectPileUp && fAnalysisUtil->IsPileUpEvent(InputEvent())) return;
1013  }
1014 
1015  fFlowEvent->ClearFast();
1016 
1017  if(!McEvent) {
1018  AliError("ERROR: Could not retrieve MCEvent");
1019  return;
1020  }
1021  fStack = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
1022  if(!fStack){
1023  AliError("ERROR: Could not retrieve MCStack");
1024  return;
1025  }
1026 
1027  // get centrality (from AliMultSelection or AliCentrality)
1028  Double_t centr = 300;
1029  if(fDataSet==k2015 || fDataSet==k2015v6) {
1030  fMultSelection = (AliMultSelection*)aod->FindListObject("MultSelection");
1031  if(!fMultSelection) {
1032  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
1033  AliWarning("AliMultSelection object not found!");
1034  } else {
1035  centr = fMultSelection->GetMultiplicityPercentile("V0M");
1036  }
1037  } else {
1038  centr = (((AliVAODHeader*)aod->GetHeader())->GetCentralityP())->GetCentralityPercentile("V0M");
1039  }
1040  // centrality bin
1041  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1042  Int_t CenBin = -1;
1043  CenBin = GetCenBin(centr);
1044  if(CenBin==-1) return;
1045 
1046  // reconstructed
1047  for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1048 
1049  AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1050  if(!track) continue;
1051 
1052  // to select primaries
1053  Int_t lp = TMath::Abs(track->GetLabel());
1054 
1055  // general kinematic & quality cuts
1056  if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70) continue;
1057 
1058  // Double_t DCAxy = track->DCA();
1059  // Double_t DCAz = track->ZAtDCA();
1060  // if (std::abs((Int_t)DCAxy)==999 || std::abs((Int_t)DCAz)==999) {
1061  // // re-evaluate the dca as it seems to not be natively present
1062  // // allowed only for tracks inside the beam pipe
1063  // Double_t pos[3] = {-99., -99., -99.};
1064  // track->GetPosition(pos);
1065  // if(pos[0]*pos[0]+pos[1]*pos[1] <= 3.*3.) {
1066  // AliAODTrack copy(*track); // stack copy
1067  // Double_t b[2] = {-99., -99.};
1068  // Double_t bCov[3] = {-99., -99., -99.};
1069  // if(copy.PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov)) {
1070  // DCAxy = b[0];
1071  // DCAz = b[1];
1072  // }
1073  // }
1074  // }
1075 
1076  // test filter bits
1077  if (((AliAODMCParticle*)fStack->At(lp))->IsPhysicalPrimary()) {
1078  if (track->TestFilterBit(32)) fPtSpecFB32[0][CenBin]->Fill(track->Pt());
1079  if (track->TestFilterBit(96)) fPtSpecFB96[0][CenBin]->Fill(track->Pt());
1080  if (track->TestFilterBit(128)) fPtSpecFB128[0][CenBin]->Fill(track->Pt());
1081  if (track->TestFilterBit(768)) fPtSpecFB768[0][CenBin]->Fill(track->Pt());
1082  } else {
1083  if (track->TestFilterBit(32)) fPtSpecFB32[1][CenBin]->Fill(track->Pt());
1084  if (track->TestFilterBit(96)) fPtSpecFB96[1][CenBin]->Fill(track->Pt());
1085  if (track->TestFilterBit(128)) fPtSpecFB128[1][CenBin]->Fill(track->Pt());
1086  if (track->TestFilterBit(768)) fPtSpecFB768[1][CenBin]->Fill(track->Pt());
1087  }
1088 
1089  fCenDis->Fill(centr);
1090  }
1091 
1092  // generated (physical primaries)
1093 
1094  for(Int_t jTracks = 0; jTracks<fStack->GetEntriesFast(); jTracks++) {
1095  AliAODMCParticle *MCpart = (AliAODMCParticle*)fStack->At(jTracks);
1096  if (!MCpart) {
1097  printf("ERROR: Could not receive MC track %d\n", jTracks);
1098  continue;
1099  }
1100  // kinematic cuts
1101  if ( MCpart->Pt() < .2 || MCpart->Pt() > 20. || TMath::Abs(MCpart->Eta()) > .8 ) continue;
1102  // select charged primaries
1103  if ( MCpart->Charge() == 0. || !MCpart->IsPhysicalPrimary()) continue;
1104 
1105  fPtSpecGen[0][CenBin]->Fill(MCpart->Pt());
1106  }
1107 
1108  // fGenHeader = McEvent->GenEventHeader();
1109  // if(fGenHeader) fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1110  // printf("#reconstructed : %d (rejected from cuts %d), #MC primaries : %d (rejected from cuts %d) \n",AODPOIs,AODbads,MCPrims,MCSecos);
1111  fFlowEvent->SetReferenceMultiplicity(aod->GetNumberOfTracks());
1112  fFlowEvent->SetCentrality(centr);
1113  // if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1114  fFlowEvent->SetRun(RunNum);
1115  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1116  }
1117 
1118  if(fAnalysisType == "MCESD") {
1119 
1120  fFlowEvent->ClearFast();
1121 
1122  if(!esd) {
1123  AliError("ERROR: Could not retrieve ESDEvent");
1124  return;
1125  }
1126  if(!McEvent) {
1127  AliError("ERROR: Could not retrieve MCEvent");
1128  return;
1129  }
1130  AliStack* fStack = fMCEvent->Stack();
1131  if(!fStack) {
1132  AliError("ERROR: Could not retrieve MCStack");
1133  return;
1134  }
1135 
1136  AliESDVertex *vertex = (AliESDVertex*) esd->GetPrimaryVertex();
1137  if (!vertex) return;
1138  if (TMath::Abs(vertex->GetZ()) > 10. ) return;
1139  if (vertex->GetNContributors() < 1 ) return;
1140  AliCentrality *centrality = esd->GetCentrality();
1141  if (!centrality) return;
1142  Double_t centr = centrality->GetCentralityPercentile("V0M");
1143  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1144  Int_t CenBin = -1;
1145  if (centr>0. && centr<5.) CenBin=0;
1146  if (centr>5. && centr<10.) CenBin=1;
1147  if (centr>10. && centr<20.) CenBin=2;
1148  if (centr>20. && centr<30.) CenBin=3;
1149  if (centr>30. && centr<40.) CenBin=4;
1150  if (centr>40. && centr<50.) CenBin=5;
1151  if (centr>50. && centr<60.) CenBin=6;
1152  if (centr>60. && centr<70.) CenBin=7;
1153  if (centr>70. && centr<80.) CenBin=8;
1154  if (centr>80. && centr<90.) CenBin=9;
1155  if(CenBin==-1) return;
1156 
1157  //Generated
1158  Int_t MCPrims = 0;
1159  for ( Int_t i=0 ; i<fStack->GetNtrack() ; i++ ) {
1160 
1161  //Primaries Selection
1162  TParticle *particle = (TParticle*)fStack->Particle(i);
1163  if (!particle) continue;
1164  if (!fStack->IsPhysicalPrimary(i)) continue;
1165  if ( particle->GetPDG()->Charge() == 0.) continue;
1166 
1167  //Kinematic Cuts
1168  if ( particle->Pt()<0.2 || particle->Pt()>10. ) continue;
1169  if ( TMath::Abs(particle->Eta())>0.8 ) continue;
1170 
1171  fFlowTrack->SetPhi(particle->Phi());
1172  fFlowTrack->SetEta(particle->Eta());
1173  fFlowTrack->SetPt(particle->Pt());
1175  fFlowTrack->SetForRPSelection(kTRUE);
1177  fFlowTrack->SetForPOISelection(kFALSE);
1179  MCPrims++;
1180 
1181  fPtSpecGen[0][CenBin]->Fill(particle->Pt());
1182 
1183  }
1184 
1185  //Reconstructed
1186  Int_t ESDPrims = 0;
1187  for (Int_t i=0 ; i<esd->GetNumberOfTracks() ; i++) {
1188 
1189  //Get reconstructed track
1190  AliVTrack *vtrack = static_cast<AliVTrack*>(esd->GetTrack(i));
1191  AliESDtrack *track = dynamic_cast<AliESDtrack*>(vtrack);
1192  if (!track) continue;
1193 
1194  //Primaries selection
1195  Int_t lp = TMath::Abs(track->GetLabel());
1196  if (!fStack->IsPhysicalPrimary(lp)) continue;
1197  TParticle *particle = (TParticle*)fStack->Particle(lp);
1198  if (!particle) continue;
1199  if (particle->GetPDG()->Charge() == 0.) continue;
1200 
1201  // if(!fCutsPOI->PassesESDcuts(track)) continue;
1202 
1203  Bool_t pass = kTRUE;
1204 
1205  if(fCutTPC) {
1206  // printf("******* cutting TPC ******** \n");
1207  UShort_t ntpccls = track->GetTPCNcls();
1208  Double_t tpcchi2 = track->GetTPCchi2();
1209  if (tpcchi2<0.2 || tpcchi2 >=4.) {
1210  // printf("TPCchi2 : %e %e ",tpcchi2,track->GetTPCchi2Iter1());
1211  pass=kFALSE;
1212  }
1213  if (ntpccls < 70) {
1214  // printf("#TPCcluster : %u %u %u %u ",ntpccls,track->GetTPCNclsF(),track->GetTPCNclsFIter1(),track->GetTPCNclsIter1());
1215  pass=kFALSE;
1216  }
1217  }
1218 
1219  Float_t dcaxy=0.0;
1220  Float_t dcaz=0.0;
1221  track->GetImpactParameters(dcaxy,dcaz);
1222  if (dcaxy > 0.3 || dcaz > 0.3) {
1223  // printf("DCA : %e %e ",dcaxy,dcaz);
1224  pass=kFALSE;
1225  }
1226  if(!pass) continue;
1227 
1228  //Kinematic Cuts
1229  if ( track->Pt()<0.2 || track->Pt()>10. ) continue;
1230  if ( TMath::Abs(track->Eta())>0.8 ) continue;
1231 
1232  fFlowTrack->SetPhi(track->Phi());
1233  fFlowTrack->SetEta(track->Eta());
1234  fFlowTrack->SetPt(track->Pt());
1236  fFlowTrack->SetForRPSelection(kFALSE);
1240  ESDPrims++;
1241 
1242  }
1243 
1244  // printf("#reconstructed : %d , #MC primaries : %d \n",ESDPrims,MCPrims);
1245  fFlowEvent->SetReferenceMultiplicity(esd->GetNumberOfTracks());
1246  fFlowEvent->SetCentrality(centr);
1247  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1248  fFlowEvent->SetRun(esd->GetRunNumber());
1249  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1250 
1251  } // end of if(fAnalysisType == "MCESD")
1252 
1253  if(fAnalysisType == "MCkine") {
1254 
1255  fFlowEvent->ClearFast();
1256 
1257  AliInputEventHandler* McHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1258  if(!McHandler) {
1259  AliError("ERROR: Could not retrieve MCtruthEventHandler");
1260  return;
1261  }
1262  McEvent = McHandler->MCEvent();
1263  if(!McEvent) {
1264  AliError("ERROR: Could not retrieve MC event");
1265  return;
1266  }
1267 
1268  Int_t nTracks = McEvent->GetNumberOfTracks();
1269  // Int_t nPrimTr = McEvent->GetNumberOfPrimaries();
1270 
1271  //loop over tracks
1272  for (Int_t itrkN=0; itrkN<nTracks; itrkN++) {
1273  //get input particle
1274  AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(McEvent->GetTrack(itrkN));
1275  if (!pParticle) continue;
1276 
1277  //check if track passes the cuts
1278  if (McEvent->IsPhysicalPrimary(itrkN) && pParticle->Charge()!=0) {
1279  fFlowTrack->Set(pParticle);
1281  fFlowTrack->SetForRPSelection(kTRUE);
1286  }
1287  }// for all tracks
1288 
1289  // if monte carlo event get reaction plane from monte carlo (depends on generator)
1290  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1291  // set reference multiplicity
1292  fFlowEvent->SetReferenceMultiplicity(McEvent->GetNumberOfTracks());
1293  // tag subevents
1295  // set centrality from impact parameter
1296  Double_t ImpPar=0., CenPer=0.;
1297  fGenHeader = McEvent->GenEventHeader();
1298  if(fGenHeader){
1299  fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1300  if(fPythiaGenHeader) ImpPar = fPythiaGenHeader->GetImpactParameter();
1301  fHijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(fGenHeader);
1302  if(fHijingGenHeader) ImpPar = fHijingGenHeader->ImpactParameter();
1303  if(ImpPar) CenPer = 0.3859796743103508*pow(ImpPar,2.);
1304  if(CenPer>0. && CenPer<100.) fFlowEvent->SetCentrality(CenPer);
1305  else return;
1306  fFlowEvent->SetRun(1);
1307  }
1308 
1309  } // end of if(fAnalysisType == "MCkine")
1310 
1311  if (!fFlowEvent) return; //shuts up coverity
1312 
1313  //check final event cuts
1314  Int_t mult = fFlowEvent->NumberOfTracks();
1315  // AliInfo(Form("FlowEvent has %i tracks",mult));
1316  if (mult<fMinMult || mult>fMaxMult) {
1317  AliWarning("FlowEvent cut on multiplicity"); return;
1318  }
1319 
1320  //define dead zone
1322 
1325  if (fAfterburnerOn)
1326  {
1327  //if reaction plane not set from elsewhere randomize it before adding flow
1329  fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
1330 
1331  if(fDifferentialV2)
1333  else
1334  fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
1335  fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
1336  }
1338 
1339  //tag subEvents
1341 
1342  //do we want to serve shullfed tracks to everybody?
1344 
1345  // associate the mother particles to their daughters in the flow event (if any)
1347 
1348  //fListHistos->Print();
1349  //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
1350 
1351  //********************************************************************************************************************************
1352 
1353  if(fAnalysisType == "AOD" || fAnalysisType == "AUTOMATIC") {
1354 
1355  // PHYSICS SELECTION
1356  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1357  AliInputEventHandler *hdr = (AliInputEventHandler*)am->GetInputEventHandler();
1358 
1359  if(hdr->IsEventSelected() && AliVEvent::kAny) {
1360 
1361  Double_t centrperc = fFlowEvent->GetCentrality();
1362  Int_t cenb = (Int_t)centrperc;
1363 
1364  AliAODTracklets *trackl = aod->GetTracklets();
1365  Int_t nTracklets = trackl->GetNumberOfTracklets();
1366 
1367  AliAODVZERO *vzeroAOD = aod->GetVZEROData();
1368  Double_t multV0A = vzeroAOD->GetMTotV0A();
1369  Double_t multV0C = vzeroAOD->GetMTotV0C();
1370  for(Int_t i=0; i<64; i++) {
1371  Double_t mult = vzeroAOD->GetMultiplicity(i);
1372  fVZEROMult->Fill(RunBin+0.5,i+0.5,mult);
1373  }
1374 
1375 // AliAODForwardMult* aodForward = static_cast<AliAODForwardMult*>(aodEvent->FindListObject("Forward"));
1376 // const TH2D& d2Ndetadphi = aodForward->GetHistogram();
1377 // Int_t nEta = d2Ndetadphi.GetXaxis()->GetNbins();
1378 // Int_t nPhi = d2Ndetadphi.GetYaxis()->GetNbins();
1379 // Double_t ret = 0.;
1380 // // Loop over eta
1381 // for (Int_t iEta = 1; iEta <= nEta; iEta++) {
1382 // Int_t valid = d2Ndetadphi.GetBinContent(iEta, 0);
1383 // if (!valid) continue; // No data expected for this eta
1384 // // Loop over phi
1385 // for (Int_t iPhi = 1; i <= nPhi; i++) {
1386 // ret = d2Ndetadphi.GetBinContent(iEta, iPhi);
1387 // printf("eta %e phi %e : %e \n",d2Ndetadphi.GetXaxis()->GetBinCenter(iEta),d2Ndetadphi.GetYaxis()->GetBinCenter(iPhi),ret);
1388 // }
1389 // }
1390 
1391  AliAODZDC *aodZDC = aod->GetZDCData();
1392 
1393  Double_t energyZNC = (Double_t) (aodZDC->GetZNCEnergy());
1394  Double_t energyZPC = (Double_t) (aodZDC->GetZPCEnergy());
1395  Double_t energyZNA = (Double_t) (aodZDC->GetZNAEnergy());
1396  Double_t energyZPA = (Double_t) (aodZDC->GetZPAEnergy());
1397  Double_t energyZEM1 = (Double_t) (aodZDC->GetZEM1Energy());
1398  Double_t energyZEM2 = (Double_t) (aodZDC->GetZEM2Energy());
1399 
1400  const Double_t * towZNCraw = aodZDC->GetZNCTowerEnergy();
1401  const Double_t * towZNAraw = aodZDC->GetZNATowerEnergy();
1402 
1403  // Get centroid from ZDCs *******************************************************
1404 
1405  Double_t Enucl = (RunNum < 209122 ? 1380. : 2511.);
1406  Double_t xyZNC[2]={0.,0.}, xyZNA[2]={0.,0.};
1407  Double_t towZNC[5]={0.}, towZNA[5]={0.};
1408 
1409  Double_t ZNCcalib=1., ZNAcalib=1.;
1410  if(fUseTowerEq) {
1411  if(RunNum!=fCachedRunNum) {
1412  for(Int_t i=0; i<5; i++) {
1413  fTowerGainEq[0][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNCTower[%d][%d]",RunNum,i)));
1414  fTowerGainEq[1][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNATower[%d][%d]",RunNum,i)));
1415  }
1416  }
1417  for(Int_t i=0; i<5; i++) {
1418  if(fTowerGainEq[0][i]) towZNC[i] = towZNCraw[i]*fTowerGainEq[0][i]->GetBinContent(fTowerGainEq[0][i]->FindBin(centrperc));
1419  if(fTowerGainEq[1][i]) towZNA[i] = towZNAraw[i]*fTowerGainEq[1][i]->GetBinContent(fTowerGainEq[1][i]->FindBin(centrperc));
1420  }
1421  } else {
1422  for(Int_t i=0; i<5; i++) {
1423  towZNC[i] = towZNCraw[i];
1424  towZNA[i] = towZNAraw[i];
1425  fZNCTower[RunBin][i]->Fill(centrperc,towZNC[i]);
1426  fZNATower[RunBin][i]->Fill(centrperc,towZNA[i]);
1427  }
1428  }
1429 
1430  if(RunNum>=245829) towZNA[2] = 0.;
1431  Double_t zncEnergy=0., znaEnergy=0.;
1432  for(Int_t i=0; i<5; i++){
1433  zncEnergy += towZNC[i];
1434  znaEnergy += towZNA[i];
1435  }
1436  if(RunNum>=245829) znaEnergy *= 8./7.;
1437  fFlowEvent->SetZNCEnergy(towZNC[0]);
1438  fFlowEvent->SetZNAEnergy(towZNA[0]);
1439 
1440  const Double_t x[4] = {-1.75, 1.75, -1.75, 1.75};
1441  const Double_t y[4] = {-1.75, -1.75, 1.75, 1.75};
1442  Double_t numXZNC=0., numYZNC=0., denZNC=0., cZNC, wZNC, EZNC, SumEZNC=0.;
1443  Double_t numXZNA=0., numYZNA=0., denZNA=0., cZNA, wZNA, EZNA, SumEZNA=0., BadChOr;
1444  Bool_t fAllChONZNC=kTRUE, fAllChONZNA=kTRUE;
1445 
1446  if (fUseMCCen) {
1447  for(Int_t i=0; i<4; i++){
1448 
1449  // get energy
1450  EZNC = towZNC[i+1];
1451  fhZNSpectra->Fill(centrperc,i+0.5,EZNC);
1452 // fhZNSpectraRbR[RunBin]->Fill(centrperc,i+0.5,EZNC);
1453  if(fUseZDCSpectraCorr && EZNC>0.) {
1454  Double_t mu1 = SpecCorMu1[i]->Interpolate(centrperc);
1455  Double_t mu2 = SpecCorMu2[i]->Interpolate(centrperc);
1456  Double_t av = SpecCorAv[i]->Interpolate(centrperc);
1457  Double_t cor1 = SpecCorSi[i]->Interpolate(centrperc);
1458  EZNC = exp( (log(EZNC) - mu1 + mu2*cor1)/cor1 ) + av;
1459  fhZNSpectraCor->Fill(centrperc,i+0.5,EZNC);
1460  }
1461  if(fUseZDCSpectraCorr && EZNC<=0.) fAllChONZNC=kFALSE;
1462  SumEZNC += EZNC;
1463 
1464  // build centroid
1465  wZNC = TMath::Power(EZNC, fZDCGainAlpha);
1466  numXZNC += x[i]*wZNC;
1467  numYZNC += y[i]*wZNC;
1468  denZNC += wZNC;
1469  fhZNSpectraPow->Fill(centrperc,i+0.5,wZNC);
1470 
1471  // get energy
1472  if(i==1) {
1473  EZNA = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1474  if(fUseBadTowerCalib && fBadTowerCalibHist[cenb]) {
1475  EZNA = GetBadTowerResp(EZNA, fBadTowerCalibHist[cenb]);
1476  }
1477  } else {
1478  EZNA = towZNA[i+1];
1479  }
1480  fhZNSpectra->Fill(centrperc,i+4.5,EZNA);
1481 // fhZNSpectraRbR[RunBin]->Fill(centrperc,i+4.5,EZNA);
1482  if(fUseZDCSpectraCorr && EZNA>0.) {
1483  Double_t mu1 = SpecCorMu1[i+4]->Interpolate(centrperc);
1484  Double_t mu2 = SpecCorMu2[i+4]->Interpolate(centrperc);
1485  Double_t av = SpecCorAv[i+4]->Interpolate(centrperc);
1486  Double_t cor1 = SpecCorSi[i+4]->Interpolate(centrperc);
1487  EZNA = exp( (log(EZNA) - mu1 + mu2*cor1)/cor1 ) + av;
1488  fhZNSpectraCor->Fill(centrperc,i+4.5,EZNA);
1489  }
1490  if(fUseZDCSpectraCorr && EZNA<=0.) fAllChONZNA=kFALSE;
1491  SumEZNA += EZNA;
1492 
1493  // build centroid
1494  wZNA = TMath::Power(EZNA, fZDCGainAlpha);
1495  numXZNA += x[i]*wZNA;
1496  numYZNA += y[i]*wZNA;
1497  denZNA += wZNA;
1498  fhZNSpectraPow->Fill(centrperc,i+4.5,wZNA);
1499  }
1500  // store distribution for unfolding
1501  if(RunNum<245829) {
1502  Double_t recoE = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1503  Double_t trueE = towZNA[2];
1504  fhZNBCCorr->Fill(centrperc,trueE,recoE);
1505  }
1506  if(denZNC>0.){
1507  Double_t nSpecnC = SumEZNC/Enucl;
1508  cZNC = 1.89358-0.71262/(nSpecnC+0.71789);
1509  xyZNC[0] = cZNC*numXZNC/denZNC;
1510  xyZNC[1] = cZNC*numYZNC/denZNC;
1511  denZNC *= cZNC;
1512  }
1513  else{
1514  xyZNC[0] = xyZNC[1] = 0.;
1515  }
1516  if(denZNA>0.){
1517  Double_t nSpecnA = SumEZNA/Enucl;
1518  cZNA = 1.89358-0.71262/(nSpecnA+0.71789);
1519  xyZNA[0] = cZNA*numXZNA/denZNA;
1520  xyZNA[1] = cZNA*numYZNA/denZNA;
1521  denZNA *= cZNA;
1522  }
1523  else{
1524  xyZNA[0] = xyZNA[1] = 0.;
1525  }
1526  } else {
1527  for(Int_t i=0; i<4; i++) {
1528  if(towZNC[i+1]>0.) {
1529  wZNC = TMath::Power(towZNC[i+1], fZDCGainAlpha);
1530  numXZNC += x[i]*wZNC;
1531  numYZNC += y[i]*wZNC;
1532  denZNC += wZNC;
1533  }
1534  if(towZNA[i+1]>0.) {
1535  wZNA = TMath::Power(towZNA[i+1], fZDCGainAlpha);
1536  numXZNA += x[i]*wZNA;
1537  numYZNA += y[i]*wZNA;
1538  denZNA += wZNA;
1539  }
1540  }
1541  if(denZNC!=0) {
1542  xyZNC[0] = numXZNC/denZNC;
1543  xyZNC[1] = numYZNC/denZNC;
1544  }
1545  else{
1546  xyZNC[0] = xyZNC[1] = 999.;
1547  zncEnergy = 0.;
1548  }
1549  if(denZNA!=0) {
1550  xyZNA[0] = numXZNA/denZNA;
1551  xyZNA[1] = numYZNA/denZNA;
1552  }
1553  else{
1554  xyZNA[0] = xyZNA[1] = 999.;
1555  znaEnergy = 0.;
1556  }
1557  }
1558 
1559  if(!fAllChONZNC) denZNC=-1.;
1560  if(!fAllChONZNA) denZNA=-1.;
1561 
1562  if(denZNC>0. && pow(xyZNC[0]*xyZNC[0]+xyZNC[1]*xyZNC[1],0.5)>1.E-6) fhZNCenDis[0]->Fill(centrperc,xyZNC[0],xyZNC[1]);
1563  if(denZNA>0. && pow(xyZNA[0]*xyZNA[0]+xyZNA[1]*xyZNA[1],0.5)>1.E-6) fhZNCenDis[1]->Fill(centrperc,-xyZNA[0], xyZNA[1]);
1564 
1565  fFlowEvent->SetZDC2Qsub(xyZNC,denZNC,xyZNA,denZNA);
1566 
1567  // ******************************************************************************
1568 
1569  Double_t tdcSum = aodZDC->GetZDCTimeSum();
1570  Double_t tdcDiff = aodZDC->GetZDCTimeDiff();
1571  fhDebunch->Fill(tdcDiff, tdcSum);
1572 
1573  for(int i=0; i<5; i++){
1574  fhZNCPM[i]->Fill(towZNC[i]);
1575  if((i<4) && (towZNC[0]>0.)) fhZNCPMQiPMC[i]->Fill(towZNC[i+1]/towZNC[0]);
1576  }
1577  for(int i=0; i<5; i++){
1578  fhZNAPM[i]->Fill(towZNA[i]);
1579  if(((i<4) && towZNA[0]>0.)) fhZNAPMQiPMC[i]->Fill(towZNA[i+1]/towZNA[0]);
1580  }
1581 
1582  fhZNCvsZNA->Fill(energyZNA, energyZNC);
1583  fhZDCCvsZDCCA->Fill(energyZNA+energyZPA, energyZNC+energyZPC);
1584  fhZNCvsZPC->Fill(energyZPC, energyZNC);
1585  fhZNAvsZPA->Fill(energyZPA, energyZNA);
1586  fhZNvsZP->Fill(energyZPA+energyZPC, energyZNA+energyZNC);
1587  fhZNvsVZERO->Fill(multV0A+multV0C, energyZNC+energyZNA);
1588  fhZDCvsVZERO->Fill(multV0A+multV0C, energyZNA+energyZPA+energyZNC+energyZPC);
1589  fhZDCvsTracklets->Fill((Double_t) (nTracklets), energyZNA+energyZPA+energyZNC+energyZPC);
1590 
1591  Double_t asymmetry = -999.;
1592  if((energyZNC+energyZNA)>0.) asymmetry = (energyZNC-energyZNA)/(energyZNC+energyZNA);
1593  fhAsymm->Fill(asymmetry);
1594  fhZNAvsAsymm->Fill(asymmetry, energyZNA/1000.);
1595  fhZNCvsAsymm->Fill(asymmetry, energyZNC/1000.);
1596 
1597  fhZNCvscentrality->Fill(centrperc, energyZNC/1000.);
1598  fhZNAvscentrality->Fill(centrperc, energyZNA/1000.);
1599 
1600  } // PHYSICS SELECTION
1601 
1602  }
1603 
1604  // p) cache run number
1606 
1607  // printf("debug: NoRPs %e, NoPOIs %e, RunNum %d, Cen %e \n",fFlowEvent->GetNumberOfRPs(),fFlowEvent->GetNumberOfPOIs(),fCachedRunNum,fFlowEvent->GetCentrality());
1608 
1609  PostData(1, fFlowEvent);
1610 
1611  PostData(2, fOutput);
1612 }
1613 //________________________________________________________________________
1614 
1616 {
1617  Double_t EtC = BadTowerCalibHist->ProjectionY("",BadTowerCalibHist->GetXaxis()->FindBin(Et),BadTowerCalibHist->GetXaxis()->FindBin(Et))->GetRandom();
1618  return EtC;
1619 }
1620 
1621 //________________________________________________________________________
1622 
1624 {
1625  Int_t CenBin=-1;
1626  if (Centrality>0. && Centrality<5.) CenBin=0;
1627  if (Centrality>5. && Centrality<10.) CenBin=1;
1628  if (Centrality>10. && Centrality<20.) CenBin=2;
1629  if (Centrality>20. && Centrality<30.) CenBin=3;
1630  if (Centrality>30. && Centrality<40.) CenBin=4;
1631  if (Centrality>40. && Centrality<50.) CenBin=5;
1632  if (Centrality>50. && Centrality<60.) CenBin=6;
1633  if (Centrality>60. && Centrality<70.) CenBin=7;
1634  if (Centrality>70. && Centrality<80.) CenBin=8;
1635  if (Centrality>80. && Centrality<90.) CenBin=9;
1636  if (CenBin>=fnCen) CenBin=-1;
1637  if (fnCen==1) CenBin=0;
1638  return CenBin;
1639 } // end of AliFlowAnalysisCRC::GetCRCCenBin(Double_t Centrality)
1640 //_____________________________________________________________________________
1641 
1642 Double_t AliAnalysisTaskCRCZDC::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
1643 {
1644  // calculate sqrt of weighted distance to other vertex
1645  if (!v0 || !v1) {
1646  printf("One of vertices is not valid\n");
1647  return 0;
1648  }
1649  static TMatrixDSym vVb(3);
1650  double dist = -1;
1651  double dx = v0->GetX()-v1->GetX();
1652  double dy = v0->GetY()-v1->GetY();
1653  double dz = v0->GetZ()-v1->GetZ();
1654  double cov0[6],cov1[6];
1655  v0->GetCovarianceMatrix(cov0);
1656  v1->GetCovarianceMatrix(cov1);
1657  vVb(0,0) = cov0[0]+cov1[0];
1658  vVb(1,1) = cov0[2]+cov1[2];
1659  vVb(2,2) = cov0[5]+cov1[5];
1660  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1661  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1662  vVb.InvertFast();
1663  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
1664  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1665  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1666  return dist>0 ? TMath::Sqrt(dist) : -1;
1667 }
1668 //________________________________________________________________________
1669 
1671 {
1672  // check for multi-vertexer pile-up
1673 
1674  const int kMinPlpContrib = 5;
1675  const double kMaxPlpChi2 = 5.0;
1676  const double kMinWDist = 15;
1677 
1678  const AliVVertex* vtPrm = 0;
1679  const AliVVertex* vtPlp = 0;
1680  int nPlp = 0;
1681 
1682  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
1683  vtPrm = aod->GetPrimaryVertex();
1684  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
1685 
1686  //int bcPrim = vtPrm->GetBC();
1687 
1688  for (int ipl=0;ipl<nPlp;ipl++) {
1689  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
1690  //
1691  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1692  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1693  // int bcPlp = vtPlp->GetBC();
1694  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
1695  //
1696  double wDst = GetWDist(vtPrm,vtPlp);
1697  if (wDst<kMinWDist) continue;
1698  //
1699  return kTRUE; // pile-up: well separated vertices
1700  }
1701 
1702  return kFALSE;
1703 }
1704 
1705 //________________________________________________________________________
1707  fCutContainer->AddAt(cutsRP,0); fCutsRP=cutsRP; cutsRP->SetPOItype(0);
1708 }
1709 
1710 //________________________________________________________________________
1712  fCutContainer->AddAt(cutsPOI,1); fCutsPOI=cutsPOI; cutsPOI->SetPOItype(1);
1713 }
1714 
1715 //________________________________________________________________________
1717 {
1718  // Terminate analysis
1719  //
1720  /* if(fDebug > 1) printf(" **** AliAnalysisTaskCRCZDC::Terminate() \n");
1721 
1722  //fOutput = dynamic_cast<TList*> (GetOutputData(1));
1723  //if(!fOutput) printf("ERROR: fOutput not available\n");
1724  */
1725 }
TH2F * fhZDCCvsZDCCA
ZNC vs ZNA;.
TH2F * fhZNvsVZERO
ZNC+ZNA vs ZPC+ZPA;.
TH1F * fPtSpecFB128[2][10]
PtSpecRec FB96.
void FindDaughters(Bool_t keepDaughtersInRPselection=kFALSE)
TH3D * fhZNSpectraCor
ZNA vs. centrality.
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
void SetEta(Double_t eta)
TH1F * fPileUpCount
centrality distribution
Int_t fRunList[fCRCMaxnRun]
Definition: External.C:236
TH2F * fhZDCvsTracklets
ZDC vs VZERO;.
static const Int_t fCRCnTow
TH1F * fhZNCPMQiPMC[4]
ZNA PM high res.
TH2F * fhZNCvsZPC
ZDCC vs ZDCCA.
virtual void SetZDC2Qsub(Double_t *QVC, Double_t MC, Double_t *QVA, Double_t MA)
void AddFlow(Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5)
Int_t GetReferenceMultiplicity(AliVEvent *event, AliMCEvent *mcEvent)
centrality
void Set(const AliVParticle *p)
TH2F * fhZNCvsAsymm
ZNA vs asymmetry.
void SetMCReactionPlaneAngle(const AliMCEvent *mcEvent)
TH3D * fhZNCenDis[2]
Debunch;.
void SetCutsPOI(AliFlowTrackCuts *cutsPOI)
TH1F * fPtSpecFB768[2][10]
PtSpecRec FB128.
TH1F * fPtSpecGen[2][10]
list with pt spectra
virtual void UserExec(Option_t *option)
void SetHistWeightvsPhiMax(Double_t d)
TCanvas * c
Definition: TestFitELoss.C:172
TH2F * fhZNCvscentrality
ZNC vs asymmetry.
void SetZNAEnergy(Double_t const en)
void SetReferenceMultiplicity(Int_t m)
virtual Int_t GetCenBin(Double_t Centrality)
void SetCutsRP(AliFlowTrackCuts *cutsRP)
TRandom * gRandom
TH2F * fhZNCvsZNA
PMQi/PMC for ZNA.
TH2F * fhDebunch
ZDC vs N_cluster layer 1;.
TH1F * fhAsymm
ZN centroid vs centrality.
void SetCentralityCL1(Double_t c)
virtual void SetVertexPosition(Double_t *pos)
void IncrementNumberOfPOIs(Int_t poiType=1)
TList * GetQA() const
TH3D * fhZNBCCorr
ZNA vs. centrality.
void SetForRPSelection(Bool_t b=kTRUE)
const Double_t etamin
void Fill(AliFlowTrackCuts *rpCuts, AliFlowTrackCuts *poiCuts)
TH2F * fhZNAvsZPA
ZNC vs ZPC;.
int Int_t
Definition: External.C:63
TH1F * fhZNAPM[5]
ZNC PM high res.
TF1 * fMultTOFLowCut
centrality distribution
unsigned int UInt_t
Definition: External.C:33
void SetShuffleTracks(Bool_t b)
float Float_t
Definition: External.C:68
AliMultSelection * fMultSelection
void SetNITSCL1(Double_t c)
Definition: External.C:252
Definition: External.C:228
virtual void UserCreateOutputObjects()
Definition: External.C:212
static AliFlowCommonConstants * GetMaster()
TH3D * fhZNSpectra
ZNA vs. centrality.
TH2F * fhZDCvsNclu1
ZDC vs N_tracklets;.
AliFlowEventCuts * fCutsEvent
const Double_t ptmin
TH1F * fPileUpMultSelCount
centrality distribution
Double_t GetCentrality() const
Bool_t fUseTowerEq
MultSelection (RUN2 centrality estimator)
ClassImp(AliAnalysisTaskCRCZDC) AliAnalysisTaskCRCZDC
Bool_t plpMV(const AliAODEvent *aod)
void SetCentrality(Double_t c)
AliFlowTrackCuts * fCutsPOI
void SetRun(Int_t const run)
void SetSource(trackSource s)
Definition: AliFlowTrack.h:37
void SetPhi(Double_t phi)
short Short_t
Definition: External.C:23
void DefineDeadZone(Double_t etaMin, Double_t etaMax, Double_t phiMin, Double_t phiMax)
AliGenEventHeader * fGenHeader
void SetPt(Double_t pt)
TProfile * fZNCTower[fCRCMaxnRun][fCRCnTow]
Q Vectors list per run.
TList * fCRCQVecListRun[fCRCMaxnRun]
Run list.
void SetZNCEnergy(Double_t const en)
virtual void Terminate(Option_t *option)
TProfile * fZNATower[fCRCMaxnRun][fCRCnTow]
ZNC tower spectra.
void SetEvent(AliVEvent *event, AliMCEvent *mcEvent=NULL)
Bool_t IsSetMCReactionPlaneAngle() const
TClonesArray * fStack
ZNA tower spectra.
TList * GetQA() const
static const Int_t fCRCMaxnRun
ZNA vs. centrality.
Double_t GetBadTowerResp(Double_t Et, TH2D *BadTowerCalibHist)
void SetForPOISelection(Bool_t b=kTRUE)
void SetHistWeightvsPhiMin(Double_t d)
TH1F * fhZNAPMQiPMC[4]
PMQi/PMC for ZNC.
void AddV2(Double_t v2)
TH2F * fhZDCvsVZERO
ZN vs VZERO;.
void TagSubeventsInEta(Double_t etaMinA, Double_t etaMaxA, Double_t etaMinB, Double_t etaMaxB)
unsigned short UShort_t
Definition: External.C:28
TH1F * fPtSpecFB96[2][10]
PtSpecRec FB32.
TH3D * fhZNSpectraPow
ZNA vs. centrality.
const char Option_t
Definition: External.C:48
TH2F * fhZNAvsAsymm
ZN asymmetry.
TH2F * fhZNAvscentrality
ZNC vs. centrality.
bool Bool_t
Definition: External.C:53
virtual Bool_t IsSelected(TObject *obj, TObject *objmc)
TH1F * fPtSpecFB32[2][10]
PtSpecGen.
TH1F * fhZNCPM[5]
list send on output slot 0
void SetCentralityTRK(Double_t c)
Double_t GetWDist(const AliVVertex *v0, const AliVVertex *v1)
Bool_t fCutTPC
PtSpecRec FB768.
AliFlowTrackCuts * fCutsRP
void InsertTrack(AliFlowTrack *)
virtual void ClearFast()
Bool_t fUseBadTowerCalib
list for storing calib files
AliAnalysisUtils * fAnalysisUtil
Event selection.
TH2F * fhZNvsZP
ZNA vs ZPA;.
AliGenHijingEventHeader * fHijingGenHeader
const Double_t phimin
AliGenPythiaEventHeader * fPythiaGenHeader
Int_t NumberOfTracks() const