AliPhysics  2e7186c (2e7186c)
 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 "TProfile3D.h"
40 
41 #include "AliAnalysisManager.h"
42 #include "AliInputEventHandler.h"
43 #include "AliVEvent.h"
44 #include "AliESD.h"
45 #include "AliESDEvent.h"
46 #include "AliESDHeader.h"
47 #include "AliESDInputHandler.h"
48 #include "AliESDZDC.h"
49 #include "AliMultiplicity.h"
50 #include "AliAnalysisUtils.h"
51 #include "AliAODHandler.h"
52 #include "AliAODTrack.h"
53 #include "AliAODEvent.h"
54 #include "AliAODHeader.h"
55 #include "AliAODVertex.h"
56 #include "AliAODVZERO.h"
57 #include "AliAODZDC.h"
58 #include "AliAODMCHeader.h"
59 #include "AliMCEventHandler.h"
60 #include "AliMCEvent.h"
61 #include "AliHeader.h"
62 #include "AliVParticle.h"
63 #include "AliStack.h"
64 #include "AliAODMCParticle.h"
65 #include "AliAnalysisTaskSE.h"
66 #include "AliGenEventHeader.h"
67 #include "AliPhysicsSelectionTask.h"
68 #include "AliPhysicsSelection.h"
69 #include "AliBackgroundSelection.h"
70 #include "AliTriggerAnalysis.h"
71 #include "AliCentrality.h"
72 #include "AliAnalysisTaskCRCZDC.h"
73 #include "AliMultSelection.h"
74 #include "AliLumiTools.h"
75 
76 // ALICE Correction Framework
77 #include "AliCFManager.h"
78 
79 // Interface to Event generators to get Reaction Plane Angle
80 #include "AliGenCocktailEventHeader.h"
81 #include "AliGenPythiaEventHeader.h"
82 #include "AliGenHijingEventHeader.h"
83 #include "AliGenGeVSimEventHeader.h"
84 #include "AliGenEposEventHeader.h"
85 
86 // Interface to Load short life particles
87 #include "TObjArray.h"
88 #include "AliFlowCandidateTrack.h"
89 
90 // Interface to make the Flow Event Simple used in the flow analysis methods
91 #include "AliFlowEvent.h"
92 #include "AliFlowTrackCuts.h"
93 #include "AliFlowEventCuts.h"
94 #include "AliFlowCommonConstants.h"
95 
97 
98 //________________________________________________________________________
101 fAnalysisType("AUTOMATIC"),
102 fRPType(""),
103 fCFManager1(NULL),
104 fCFManager2(NULL),
105 fCutsEvent(NULL),
106 fCutsRP(NULL),
107 fCutsPOI(NULL),
108 fCutContainer(new TList()),
109 fQAList(NULL),
110 fMinMult(0),
111 fMaxMult(10000000),
112 fMinA(-1.0),
113 fMaxA(-0.01),
114 fMinB(0.01),
115 fMaxB(1.0),
116 fGenHeader(NULL),
117 fPythiaGenHeader(NULL),
118 fHijingGenHeader(NULL),
119 fFlowTrack(NULL),
120 fAnalysisUtil(NULL),
121 fQAon(kFALSE),
122 fLoadCandidates(kFALSE),
123 fNbinsMult(10000),
124 fNbinsPt(100),
125 fNbinsPhi(100),
126 fNbinsEta(200),
127 fNbinsQ(500),
128 fNbinsMass(1),
129 fMultMin(0.),
130 fMultMax(10000.),
131 fPtMin(0.),
132 fPtMax(10.),
133 fPhiMin(0.),
134 fPhiMax(TMath::TwoPi()),
135 fEtaMin(-5.),
136 fEtaMax(5.),
137 fQMin(0.),
138 fQMax(3.),
139 fMassMin(-1.),
140 fMassMax(0.),
141 fHistWeightvsPhiMin(0.),
142 fHistWeightvsPhiMax(3.),
143 fExcludedEtaMin(0.),
144 fExcludedEtaMax(0.),
145 fExcludedPhiMin(0.),
146 fExcludedPhiMax(0.),
147 fAfterburnerOn(kFALSE),
148 fNonFlowNumberOfTrackClones(0),
149 fV1(0.),
150 fV2(0.),
151 fV3(0.),
152 fV4(0.),
153 fV5(0.),
154 fDifferentialV2(0),
155 fFlowEvent(NULL),
156 fShuffleTracks(kFALSE),
157 fMyTRandom3(NULL),
158 fAnalysisInput(kAOD),
159 fIsMCInput(kFALSE),
160 fUseMCCen(kTRUE),
161 fRejectPileUp(kTRUE),
162 fRejectPileUpTight(kFALSE),
163 fCentrLowLim(0.),
164 fCentrUpLim(100.),
165 fCentrEstimator(kV0M),
166 fOutput(0x0),
167 fhZNCvsZNA(0x0),
168 fhZDCCvsZDCCA(0x0),
169 fhZNCvsZPC(0x0),
170 fhZNAvsZPA(0x0),
171 fhZNvsZP(0x0),
172 fhZNvsVZERO(0x0),
173 fhZDCvsVZERO(0x0),
174 fhZDCvsTracklets(0x0),
175 fhZDCvsNclu1(0x0),
176 fhDebunch(0x0),
177 fhZNCcentroid(0x0),
178 fhZNAcentroid(0x0),
179 fhAsymm(0x0),
180 fhZNAvsAsymm(0x0),
181 fhZNCvsAsymm(0x0),
182 fhZNCvscentrality(0x0),
183 fhZNAvscentrality(0x0),
184 fCRCnRun(0),
185 fZDCGainAlpha(0.395),
186 fDataSet(kAny),
187 fStack(0x0),
188 fCutTPC(kFALSE),
189 fCenDis(0x0),
190 fMultSelection(0x0),
191 fPileUpCount(0x0),
192 fPileUpMultSelCount(0x0),
193 fMultTOFLowCut(0x0),
194 fMultTOFHighCut(0x0),
195 fUseTowerEq(kFALSE),
196 fTowerEqList(NULL),
197 fBadTowerCalibList(NULL),
198 fSpectraMCList(NULL),
199 fBadTowerStuffList(NULL),
200 fCachedRunNum(0),
201 fhZNSpectra(0x0),
202 fhZNBCCorr(0x0)
203 {
204  for(int i=0; i<5; i++){
205  fhZNCPM[i] = 0x0;
206  fhZNAPM[i] = 0x0;
207  }
208  for(int i=0; i<4; i++){
209  fhZNCPMQiPMC[i] = 0x0;
210  fhZNAPMQiPMC[i] = 0x0;
211  }
212  for(Int_t r=0; r<fCRCMaxnRun; r++) {
213  fRunList[r] = 0;
214  }
215  for(Int_t c=0; c<2; c++) {
216  for(Int_t i=0; i<5; i++) {
217  fTowerGainEq[c][i] = NULL;
218  }
219  }
220  for(Int_t c=0; c<100; c++) {
221  fBadTowerCalibHist[c] = NULL;
222  }
223  this->InitializeRunArrays();
224  fMyTRandom3 = new TRandom3(1);
225  gRandom->SetSeed(fMyTRandom3->Integer(65539));
226  for(Int_t j=0; j<2; j++) {
227  for(Int_t c=0; c<10; c++) {
228  fPtSpecGen[j][c] = NULL;
229  fPtSpecFB32[j][c] = NULL;
230  fPtSpecFB96[j][c] = NULL;
231  fPtSpecFB128[j][c] = NULL;
232  fPtSpecFB768[j][c] = NULL;
233  }
234  }
235 }
236 
237 //________________________________________________________________________
238 AliAnalysisTaskCRCZDC::AliAnalysisTaskCRCZDC(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates):
239 AliAnalysisTaskSE(name),
240 fAnalysisType("AUTOMATIC"),
241 fRPType(RPtype),
242 fCFManager1(NULL),
243 fCFManager2(NULL),
244 fCutsEvent(NULL),
245 fCutsRP(NULL),
246 fCutsPOI(NULL),
247 fCutContainer(new TList()),
248 fQAList(NULL),
249 fMinMult(0),
250 fMaxMult(10000000),
251 fMinA(-1.0),
252 fMaxA(-0.01),
253 fMinB(0.01),
254 fMaxB(1.0),
255 fQAon(on),
256 fLoadCandidates(bCandidates),
257 fNbinsMult(10000),
258 fNbinsPt(100),
259 fNbinsPhi(100),
260 fNbinsEta(200),
261 fNbinsQ(500),
262 fNbinsMass(1),
263 fMultMin(0.),
264 fMultMax(10000.),
265 fPtMin(0.),
266 fPtMax(10.),
267 fPhiMin(0.),
268 fPhiMax(TMath::TwoPi()),
269 fEtaMin(-5.),
270 fEtaMax(5.),
271 fQMin(0.),
272 fQMax(3.),
273 fMassMin(-1.),
274 fMassMax(0.),
275 fHistWeightvsPhiMin(0.),
276 fHistWeightvsPhiMax(3.),
277 fExcludedEtaMin(0.),
278 fExcludedEtaMax(0.),
279 fExcludedPhiMin(0.),
280 fExcludedPhiMax(0.),
281 fAfterburnerOn(kFALSE),
282 fNonFlowNumberOfTrackClones(0),
283 fV1(0.),
284 fV2(0.),
285 fV3(0.),
286 fV4(0.),
287 fV5(0.),
288 fDifferentialV2(0),
289 fFlowEvent(NULL),
290 fShuffleTracks(kFALSE),
291 fMyTRandom3(NULL),
292 fAnalysisInput(kAOD),
293 fIsMCInput(kFALSE),
294 fUseMCCen(kTRUE),
295 fRejectPileUp(kTRUE),
296 fRejectPileUpTight(kFALSE),
297 fCentrLowLim(0.),
298 fCentrUpLim(100.),
299 fCentrEstimator(kV0M),
300 fOutput(0x0),
301 fhZNCvsZNA(0x0),
302 fhZDCCvsZDCCA(0x0),
303 fhZNCvsZPC(0x0),
304 fhZNAvsZPA(0x0),
305 fhZNvsZP(0x0),
306 fhZNvsVZERO(0x0),
307 fhZDCvsVZERO(0x0),
308 fhZDCvsTracklets(0x0),
309 fhZDCvsNclu1(0x0),
310 fhDebunch(0x0),
311 fhZNCcentroid(0x0),
312 fhZNAcentroid(0x0),
313 fhAsymm(0x0),
314 fhZNAvsAsymm(0x0),
315 fhZNCvsAsymm(0x0),
316 fhZNCvscentrality(0x0),
317 fhZNAvscentrality(0x0),
318 fDataSet(kAny),
319 fCRCnRun(0),
320 fZDCGainAlpha(0.395),
321 fGenHeader(NULL),
322 fPythiaGenHeader(NULL),
323 fHijingGenHeader(NULL),
324 fFlowTrack(NULL),
325 fAnalysisUtil(NULL),
326 fStack(0x0),
327 fCutTPC(kFALSE),
328 fCenDis(0x0),
329 fMultSelection(0x0),
330 fPileUpCount(0x0),
331 fPileUpMultSelCount(0x0),
332 fMultTOFLowCut(0x0),
333 fMultTOFHighCut(0x0),
334 fUseTowerEq(kFALSE),
335 fTowerEqList(NULL),
336 fBadTowerCalibList(NULL),
337 fCachedRunNum(0),
338 fhZNSpectra(0x0),
339 fhZNBCCorr(0x0)
340 {
341 
342  for(int i=0; i<5; i++){
343  fhZNCPM[i] = 0x0;
344  fhZNAPM[i] = 0x0;
345  }
346  for(int i=0; i<4; i++){
347  fhZNCPMQiPMC[i] = 0x0;
348  fhZNAPMQiPMC[i] = 0x0;
349  }
350  for(Int_t r=0; r<fCRCMaxnRun; r++) {
351  fRunList[r] = 0;
352  }
353  for(Int_t c=0; c<2; c++) {
354  for(Int_t i=0; i<5; i++) {
355  fTowerGainEq[c][i] = NULL;
356  }
357  }
358  for(Int_t c=0; c<100; c++) {
359  fBadTowerCalibHist[c] = NULL;
360  }
361  this->InitializeRunArrays();
362  fMyTRandom3 = new TRandom3(iseed);
363  gRandom->SetSeed(fMyTRandom3->Integer(65539));
364 
365  DefineInput(0, TChain::Class());
366  // Define output slots here
367  // Define here the flow event output
368  DefineOutput(1, AliFlowEventSimple::Class());
369  DefineOutput(2, TList::Class());
370 
371  for(Int_t j=0; j<2; j++) {
372  for(Int_t c=0; c<10; c++) {
373  fPtSpecGen[j][c] = NULL;
374  fPtSpecFB32[j][c] = NULL;
375  fPtSpecFB96[j][c] = NULL;
376  fPtSpecFB128[j][c] = NULL;
377  fPtSpecFB768[j][c] = NULL;
378  }
379  }
380 
381 }
382 
383 //________________________________________________________________________
385 {
386  // Destructor
387  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
388  delete fOutput; fOutput=0;
389  }
390  delete fMyTRandom3;
391  delete fFlowEvent;
392  delete fFlowTrack;
393  delete fCutsEvent;
394  if (fTowerEqList) delete fTowerEqList;
396  if (fAnalysisUtil) delete fAnalysisUtil;
397  if (fQAList) delete fQAList;
398  if (fCutContainer) fCutContainer->Delete(); delete fCutContainer;
399 }
400 
401 //________________________________________________________________________
403 {
404  for(Int_t r=0;r<fCRCMaxnRun;r++) {
405  fCRCQVecListRun[r] = NULL;
406  for(Int_t k=0;k<fCRCnTow;k++) {
407  fZNCTower[r][k] = NULL;
408  fZNATower[r][k] = NULL;
409  }
410  }
411  // for(Int_t i=0;i<fnCen;i++) {
412  // fPtPhiEtaRbRFB128[r][i] = NULL;
413  // fPtPhiEtaRbRFB768[r][i] = NULL;
414  // }
415  // }
416 }
417 
418 //________________________________________________________________________
420 {
421  // Create the output containers
422 
423  if (!(fAnalysisType == "AOD" || fAnalysisType == "MCkine" || fAnalysisType == "MCAOD" || fAnalysisType == "AUTOMATIC" || fAnalysisType == "MCESD"))
424  {
425  AliError("WRONG ANALYSIS TYPE! only MCESD, MCkine, MCAOD, AOD and AUTOMATIC are allowed.");
426  exit(1);
427  }
428 
429  //set the common constants
432  cc->SetNbinsPt(fNbinsPt);
433  cc->SetNbinsPhi(fNbinsPhi);
434  cc->SetNbinsEta(fNbinsEta);
435  cc->SetNbinsQ(fNbinsQ);
437  cc->SetMultMin(fMultMin);
438  cc->SetMultMax(fMultMax);
439  cc->SetPtMin(fPtMin);
440  cc->SetPtMax(fPtMax);
441  cc->SetPhiMin(fPhiMin);
442  cc->SetPhiMax(fPhiMax);
443  cc->SetEtaMin(fEtaMin);
444  cc->SetEtaMax(fEtaMax);
445  cc->SetQMin(fQMin);
446  cc->SetQMax(fQMax);
447  cc->SetMassMin(fMassMin);
448  cc->SetMassMax(fMassMax);
451 
452  fFlowEvent = new AliFlowEvent(20000);
453  fFlowTrack = new AliFlowTrack();
454 
455  //printf(" AliAnalysisTaskCRCZDC::UserCreateOutputObjects()\n\n");
456  fOutput = new TList();
457  fOutput->SetOwner(kTRUE);
458  //fOutput->SetName("output");
459 
460  if (fQAon) {
461  fQAList = new TList();
462  fQAList->SetOwner(kTRUE);
463  fQAList->SetName("AliFlowEventCuts QA");
464  if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA()); //0
465  if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA()); //1
466  if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA()); //2
467  fOutput->Add(fQAList);
468  }
469 
470  fCenDis = new TH1F("fCenDis", "fCenDis", 100, 0., 100.);
471  fOutput->Add(fCenDis);
472  fPileUpCount = new TH1F("fPileUpCount", "fPileUpCount", 9, 0., 9.);
473  fPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
474  fPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
475  fPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultiplicityComb08");
476  fPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
477  fPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>7.5");
478  fPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
479  fPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
480  fPileUpCount->GetXaxis()->SetBinLabel(8,"multESDTPCDif");
481  fPileUpCount->GetXaxis()->SetBinLabel(9,"extraPileUpMultSel");
482  fOutput->Add(fPileUpCount);
483  fPileUpMultSelCount = new TH1F("fPileUpMultSelCount", "fPileUpMultSelCount", 8, 0., 8.);
484  fPileUpMultSelCount->GetXaxis()->SetBinLabel(1,"IsNotPileup");
485  fPileUpMultSelCount->GetXaxis()->SetBinLabel(2,"IsNotPileupMV");
486  fPileUpMultSelCount->GetXaxis()->SetBinLabel(3,"IsNotPileupInMultBins");
487  fPileUpMultSelCount->GetXaxis()->SetBinLabel(4,"InconsistentVertices");
488  fPileUpMultSelCount->GetXaxis()->SetBinLabel(5,"TrackletVsCluster");
489  fPileUpMultSelCount->GetXaxis()->SetBinLabel(6,"AsymmetricInVZERO");
490  fPileUpMultSelCount->GetXaxis()->SetBinLabel(7,"IncompleteDAQ");
491  fPileUpMultSelCount->GetXaxis()->SetBinLabel(8,"GoodVertex2016");
493 
494  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);
495  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);
496  fOutput->Add(fMultTOFLowCut);
497  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);
498  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);
499  fOutput->Add(fMultTOFHighCut);
500 
501  for(Int_t c=0; c<2; c++) {
502  for(Int_t i=0; i<5; i++) {
503  fTowerGainEq[c][i] = new TH1D();
504  fOutput->Add(fTowerGainEq[c][i]);
505  }
506  }
507 
508  fBadTowerStuffList = new TList();
509  fBadTowerStuffList->SetOwner(kTRUE);
510  fBadTowerStuffList->SetName("BadTowerCalib");
512 
513  if(fBadTowerCalibList) {
514  for(Int_t c=0; c<100; c++) {
515  fBadTowerCalibHist[c] = (TH2D*)fBadTowerCalibList->FindObject(Form("TH2Resp[%d]",c));
517  }
518  }
519 
520  fhZNSpectra = new TH3D("fhZNSpectra","fhZNSpectra",100,0.,100.,8,0.,8.,1000,0.,1.E5);
521  fOutput->Add(fhZNSpectra);
522  fhZNBCCorr = new TH3D("fhZNBCCorr","fhZNBCCorr",100,0.,100.,500,0.,1.E5,500,0.,1.E5);
523  fOutput->Add(fhZNBCCorr);
524 
525  if(fAnalysisType == "MCAOD") {
526 
527  fSpectraMCList = new TList();
528  fSpectraMCList->SetOwner(kTRUE);
529  fSpectraMCList->SetName("Spectra");
530  fOutput->Add(fSpectraMCList);
531 
532  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.};
533  for(Int_t j=0; j<2; j++) {
534  for(Int_t c=0; c<10; c++) {
535  fPtSpecGen[j][c] = new TH1F(Form("fPtSpecGen[%d][%d]",j,c), Form("fPtSpecGen[%d][%d]",j,c), 23, xmin);
536  fSpectraMCList->Add(fPtSpecGen[j][c]);
537  fPtSpecFB32[j][c] = new TH1F(Form("fPtSpecFB32[%d][%d]",j,c), Form("fPtSpecFB32[%d][%d]",j,c), 23, xmin);
538  fSpectraMCList->Add(fPtSpecFB32[j][c]);
539  fPtSpecFB96[j][c] = new TH1F(Form("fPtSpecFB96[%d][%d]",j,c), Form("fPtSpecFB96[%d][%d]",j,c), 23, xmin);
540  fSpectraMCList->Add(fPtSpecFB96[j][c]);
541  fPtSpecFB128[j][c] = new TH1F(Form("fPtSpecFB128[%d][%d]",j,c), Form("fPtSpecFB128[%d][%d]",j,c), 23, xmin);
542  fSpectraMCList->Add(fPtSpecFB128[j][c]);
543  fPtSpecFB768[j][c] = new TH1F(Form("fPtSpecFB768[%d][%d]",j,c), Form("fPtSpecFB768[%d][%d]",j,c), 23, xmin);
544  fSpectraMCList->Add(fPtSpecFB768[j][c]);
545  }
546  }
547  }
548 
549  fAnalysisUtil = new AliAnalysisUtils();
550  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
551  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
552 
553  for(int i=0; i<5; i++){
554  char hname[20];
555  sprintf(hname,"hZNCPM%d",i);
556  fhZNCPM[i] = new TH1F(hname, hname, 200, -50., 140000);
557  fOutput->Add(fhZNCPM[i]);
558  //
559  sprintf(hname,"hZNAPM%d",i);
560  fhZNAPM[i] = new TH1F(hname, hname, 200, -50., 140000);
561  fOutput->Add(fhZNAPM[i]);
562  //
563  if(i<4){
564  //
565  char hnamenc[20];
566  sprintf(hnamenc, "hZNCPMQ%dPMC",i+1);
567  fhZNCPMQiPMC[i] = new TH1F(hnamenc, hnamenc, 100, 0., 1.);
568  fOutput->Add(fhZNCPMQiPMC[i]);
569  //
570  char hnamena[20];
571  sprintf(hnamena, "hZNAPMQ%dPMC",i+1);
572  fhZNAPMQiPMC[i] = new TH1F(hnamena, hnamena, 100, 0., 1.);
573  fOutput->Add(fhZNAPMQiPMC[i]);
574  }
575  }
576 
577  fhZNCvsZNA = new TH2F("hZNCvsZNA","hZNCvsZNA",200,-50.,140000,200,-50.,140000);
578  fOutput->Add(fhZNCvsZNA);
579  fhZDCCvsZDCCA = new TH2F("hZDCCvsZDCCA","hZDCCvsZDCCA",200,0.,180000.,200,0.,200000.);
580  fOutput->Add(fhZDCCvsZDCCA);
581  fhZNCvsZPC = new TH2F("hZNCvsZPC","hZNCvsZPC",200,-50.,50000,200,-50.,140000);
582  fOutput->Add(fhZNCvsZPC);
583  fhZNAvsZPA = new TH2F("hZNAvsZPA","hZNAvsZPA",200,-50.,50000,200,-50.,140000);
584  fOutput->Add(fhZNAvsZPA);
585  fhZNvsZP = new TH2F("hZNvsZP","hZNvsZP",200,-50.,80000,200,-50.,200000);
586  fOutput->Add(fhZNvsZP);
587  fhZNvsVZERO = new TH2F("hZNvsVZERO","hZNvsVZERO",250,0.,25000.,200,0.,200000.);
588  fOutput->Add(fhZNvsVZERO);
589  fhZDCvsVZERO = new TH2F("hZDCvsVZERO","hZDCvsVZERO",250,0.,25000.,250,0.,250000.);
590  fOutput->Add(fhZDCvsVZERO);
591  fhZDCvsTracklets = new TH2F("hZDCvsTracklets","hZDCvsTracklets",200,0.,4000.,250,0.,250000.);
593  fhZDCvsNclu1 = new TH2F("hZDCvsNclu1", "hZDCvsNclu1", 200, 0.,8000.,200,0.,250000.);
594  fOutput->Add(fhZDCvsNclu1);
595  fhDebunch = new TH2F("hDebunch","hDebunch",240,-100.,-40.,240,-30.,30.);
596  fOutput->Add(fhDebunch);
597  fhZNCcentroid = new TH2F("hZNCcentroid","hZNCcentroid",100,-3.5,3.5,100,-3.5,3.5);
598  fOutput->Add(fhZNCcentroid);
599  fhZNAcentroid = new TH2F("hZNAcentroid","hZNAcentroid",100,-3.5,3.5,100,-3.5,3.5);
600  fOutput->Add(fhZNAcentroid);
601 
602  fhAsymm = new TH1F("hAsymm" , "Asimmetry ",200,-1.,1.);
603  fOutput->Add(fhAsymm);
604  fhZNAvsAsymm = new TH2F("hZNAvsAsymm","ZNA vs. asymm.",200,-1.,1.,200,0.,80.);
605  fOutput->Add(fhZNAvsAsymm);
606  fhZNCvsAsymm = new TH2F("hZNCvsAsymm","ZNC vs. asymm.",200,-1.,1.,200,0.,80.);
607  fOutput->Add(fhZNCvsAsymm);
608 
609  fhZNCvscentrality = new TH2F("hZNCvscentrality","ZNC vs. centrality",100,0.,100.,100,0.,100.);
611  fhZNAvscentrality = new TH2F("hZNAvscentrality","ZNA vs. centrality",100,0.,100.,100,0.,100.);
613 
614  //********************************************************************
615 
616  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};
617 
618  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};
619 
620  // 12 low IR: 244917, 244918, 244975, 244980, 244982, 244983, 245064, 245066, 245068, 246390, 246391, 246392
621  // 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
622 
623  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};
624 
625  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};
626 
627  if(fDataSet==k2010) {fCRCnRun=92;}
628  if(fDataSet==k2011) {fCRCnRun=119;}
629  if(fDataSet==k2015) {fCRCnRun=90;}
630  if(fDataSet==k2015v6) {fCRCnRun=91;}
631  if(fDataSet==kAny) {fCRCnRun=1;}
632 
633  Int_t d=0;
634  for(Int_t r=0; r<fCRCnRun; r++) {
635  if(fDataSet==k2010) fRunList[d] = dRun10h[r];
636  if(fDataSet==k2011) fRunList[d] = dRun11h[r];
637  if(fDataSet==k2015) fRunList[d] = dRun15o[r];
638  if(fDataSet==k2015v6) fRunList[d] = dRun15ov6[r];
639  if(fDataSet==kAny) fRunList[d] = 1;
640  d++;
641  }
642 
643  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.};
644  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.};
645  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};
646 
647  for(Int_t r=0;r<fCRCnRun;r++) {
648  fCRCQVecListRun[r] = new TList();
649  fCRCQVecListRun[r]->SetName(Form("Run %d",fRunList[r]));
650  fCRCQVecListRun[r]->SetOwner(kTRUE);
651  fOutput->Add(fCRCQVecListRun[r]);
652 
653  for(Int_t k=0;k<fCRCnTow;k++) {
654  fZNCTower[r][k] = new TProfile(Form("fZNCTower[%d][%d]",fRunList[r],k),Form("fZNCTower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
655  fZNCTower[r][k]->Sumw2();
656  fCRCQVecListRun[r]->Add(fZNCTower[r][k]);
657  fZNATower[r][k] = new TProfile(Form("fZNATower[%d][%d]",fRunList[r],k),Form("fZNATower[%d][%d]",fRunList[r],k),100,0.,100.,"s");
658  fZNATower[r][k]->Sumw2();
659  fCRCQVecListRun[r]->Add(fZNATower[r][k]);
660  }
661 
662  // for(Int_t i=0;i<fnCen;i++) {
663  // fPtPhiEtaRbRFB128[r][i] = new TH3F(Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
664  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB128[r][i]);
665  // fPtPhiEtaRbRFB768[r][i] = new TH3F(Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
666  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB768[r][i]);
667  // }
668  }
669 
670  PostData(2, fOutput);
671 }
672 
673 //________________________________________________________________________
675 {
676  // Execute analysis for current event:
677  AliMCEvent* McEvent = MCEvent();
678  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
679  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
680  // AliMultiplicity* myTracklets = NULL;
681  // AliESDPmdTrack* pmdtracks = NULL;
682  // int availableINslot=1;
683 
684  if (!(fCutsRP&&fCutsPOI&&fCutsEvent)) {
685  AliError("cuts not set");
686  return;
687  }
688 
689  Int_t RunBin=-1, bin=0;
690  Int_t RunNum = aod->GetRunNumber();
691  for(Int_t c=0;c<fCRCnRun;c++) {
692  if(fRunList[c]==RunNum) RunBin=bin;
693  else bin++;
694  }
695  if(RunBin==-1) return;
696  if(fDataSet==kAny) RunBin=0;
697 
698  //DEFAULT - automatically takes care of everything
699  if (fAnalysisType == "AUTOMATIC") {
700 
701  // get centrality
702  Double_t centrV0M=300, centrCL1=300, centrCL0=300, centrTRK=300;
703  if(fDataSet!=k2015 && fDataSet!=k2015v6) {
704  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
705  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
706  centrCL0 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL0");
707  centrTRK = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("TRK");
708  } else {
709  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
710  if(!fMultSelection) {
711  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
712  AliWarning("AliMultSelection object not found!");
713  }else{
714  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
715  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
716  centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
717  centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
718  }
719  }
720 
721  //check event cuts
722  if (InputEvent()) {
723  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
724  if(fRejectPileUp) {
725  if(fDataSet!=k2015 && fDataSet!=k2015v6) {
726 
727  Bool_t BisPileup=kFALSE;
728 
729  // check anyway pileup
730  if (plpMV(aod)) {
731  fPileUpCount->Fill(0.5);
732  BisPileup=kTRUE;
733  }
734 
735  Short_t isPileup = aod->IsPileupFromSPD(3);
736  if (isPileup != 0) {
737  fPileUpCount->Fill(1.5);
738  BisPileup=kTRUE;
739  }
740 
741  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
742  fPileUpCount->Fill(2.5);
743  BisPileup=kTRUE;
744  }
745 
746  if (aod->IsIncompleteDAQ()) {
747  fPileUpCount->Fill(3.5);
748  BisPileup=kTRUE;
749  }
750 
751  // check vertex consistency
752  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
753  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
754 
755  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
756  fPileUpCount->Fill(5.5);
757  BisPileup=kTRUE;
758  }
759 
760  double covTrc[6], covSPD[6];
761  vtTrc->GetCovarianceMatrix(covTrc);
762  vtSPD->GetCovarianceMatrix(covSPD);
763 
764  double dz = vtTrc->GetZ() - vtSPD->GetZ();
765 
766  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
767  double errTrc = TMath::Sqrt(covTrc[5]);
768  double nsigTot = dz/errTot;
769  double nsigTrc = dz/errTrc;
770 
771  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
772  fPileUpCount->Fill(6.5);
773  BisPileup=kTRUE;
774  }
775 
776  if (fAnalysisUtil->IsPileUpEvent(InputEvent())) {
777  fPileUpCount->Fill(7.5);
778  BisPileup=kTRUE;
779  }
780 
781 // if(BisPileup) return;
782  } else {
783  // pileup from AliMultSelection
784  if(!fMultSelection->GetThisEventIsNotPileup()) fPileUpMultSelCount->Fill(0.5);
785  if(!fMultSelection->GetThisEventIsNotPileupMV()) fPileUpMultSelCount->Fill(1.5);
786  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) fPileUpMultSelCount->Fill(2.5);
787  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) fPileUpMultSelCount->Fill(3.5);
788  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) fPileUpMultSelCount->Fill(4.5);
789  if(!fMultSelection->GetThisEventIsNotAsymmetricInVZERO()) fPileUpMultSelCount->Fill(5.5);
790  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) fPileUpMultSelCount->Fill(6.5);
791  if(!fMultSelection->GetThisEventHasGoodVertex2016()) fPileUpMultSelCount->Fill(7.5);
792 
793  Bool_t BisPileup=kFALSE;
794 
795  // pile-up a la Dobrin for LHC15o
796  if (plpMV(aod)) {
797  fPileUpCount->Fill(0.5);
798  BisPileup=kTRUE;
799  }
800 
801  Short_t isPileup = aod->IsPileupFromSPD(3);
802  if (isPileup != 0) {
803  fPileUpCount->Fill(1.5);
804  BisPileup=kTRUE;
805  }
806 
807  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
808  fPileUpCount->Fill(2.5);
809  BisPileup=kTRUE;
810  }
811 
812  if (aod->IsIncompleteDAQ()) {
813  fPileUpCount->Fill(3.5);
814  BisPileup=kTRUE;
815  }
816 
817  if(fabs(centrV0M-centrCL1)>7.5) {
818  fPileUpCount->Fill(4.5);
819  BisPileup=kTRUE;
820  }
821 
822  // check vertex consistency
823  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
824  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
825 
826  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
827  fPileUpCount->Fill(5.5);
828  BisPileup=kTRUE;
829  }
830 
831  double covTrc[6], covSPD[6];
832  vtTrc->GetCovarianceMatrix(covTrc);
833  vtSPD->GetCovarianceMatrix(covSPD);
834 
835  double dz = vtTrc->GetZ() - vtSPD->GetZ();
836 
837  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
838  double errTrc = TMath::Sqrt(covTrc[5]);
839  double nsigTot = dz/errTot;
840  double nsigTrc = dz/errTrc;
841 
842  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
843  fPileUpCount->Fill(6.5);
844  BisPileup=kTRUE;
845  }
846 
847  // cuts on tracks
848  const Int_t nTracks = aod->GetNumberOfTracks();
849  Int_t multEsd = ((AliAODHeader*)aod->GetHeader())->GetNumberOfESDTracks();
850 
851  Int_t multTrk = 0;
852  Int_t multTrkBefC = 0;
853  Int_t multTrkTOFBefC = 0;
854  Int_t multTPC = 0;
855 
856  for (Int_t it = 0; it < nTracks; it++) {
857 
858  AliAODTrack* aodTrk = (AliAODTrack*)aod->GetTrack(it);
859  if (!aodTrk){
860  delete aodTrk;
861  continue;
862  }
863 
864 // if (aodTrk->TestFilterBit(32)){
865 // multTrkBefC++;
866 //
867 // if ( TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
868 // multTrkTOFBefC++;
869 //
870 // if ((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
871 // multTrk++;
872 // }
873 
874  if (aodTrk->TestFilterBit(128))
875  multTPC++;
876 
877  } // end of for (Int_t it = 0; it < nTracks; it++)
878 
879  Double_t multTPCn = multTPC;
880  Double_t multEsdn = multEsd;
881  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
882 
883  if (multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
884  fPileUpCount->Fill(7.5);
885  BisPileup=kTRUE;
886  }
887 
888  if(fRejectPileUpTight) {
889  if(BisPileup==kFALSE) {
890  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
891  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
892  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
893  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
894  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
895  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
896  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
897  if(BisPileup) fPileUpCount->Fill(8.5);
898  }
899  }
900 
901  if(BisPileup) return;
902  }
903  }
904  }
905 
906  //first attach all possible information to the cuts
907  fCutsRP->SetEvent( InputEvent(), MCEvent() ); //attach event
908  fCutsPOI->SetEvent( InputEvent(), MCEvent() );
909 
910  //then make the event
912 
914 
919  fFlowEvent->SetCentralityCL1(centrCL1);
920  fFlowEvent->SetCentralityTRK(centrTRK);
921  // fFlowEvent->SetNITSCL1(((AliVAODHeader*)aod->GetHeader())->GetNumberOfITSClusters(1));
922  Double_t SumV0=0.;
923  for(Int_t i=0; i<64; i++) {
924  if(std::isfinite(aod->GetVZEROEqMultiplicity(i))) SumV0 += aod->GetVZEROEqMultiplicity(i);
925  }
926  fFlowEvent->SetNITSCL1(SumV0);
927 
928  Double_t vtxpos[3]={0.,0.,0.};
929  vtxpos[0] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetX();
930  vtxpos[1] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetY();
931  vtxpos[2] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetZ();
932  fFlowEvent->SetVertexPosition(vtxpos);
933 
934  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
935 
936  // run-by-run QA
937  // for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
938  // AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
939  // if(!track) continue;
940  // // general kinematic & quality cuts
941  // if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70) continue;
942  // if (track->TestFilterBit(128)) fPtPhiEtaRbRFB128[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
943  // if (track->TestFilterBit(768)) fPtPhiEtaRbRFB768[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
944  // }
945  fCenDis->Fill(centrV0M);
946 
947  }
948 
949  if (fAnalysisType == "MCAOD") {
950 
951  //check event cuts
952  if (InputEvent()) {
953  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
954  if(fRejectPileUp && fAnalysisUtil->IsPileUpEvent(InputEvent())) return;
955  }
956 
958 
959  if(!McEvent) {
960  AliError("ERROR: Could not retrieve MCEvent");
961  return;
962  }
963  fStack = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
964  if(!fStack){
965  AliError("ERROR: Could not retrieve MCStack");
966  return;
967  }
968 
969  // get centrality (from AliMultSelection or AliCentrality)
970  Double_t centr = 300;
971  if(fDataSet==k2015 || fDataSet==k2015v6) {
972  fMultSelection = (AliMultSelection*)aod->FindListObject("MultSelection");
973  if(!fMultSelection) {
974  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
975  AliWarning("AliMultSelection object not found!");
976  } else {
977  centr = fMultSelection->GetMultiplicityPercentile("V0M");
978  }
979  } else {
980  centr = (((AliVAODHeader*)aod->GetHeader())->GetCentralityP())->GetCentralityPercentile("V0M");
981  }
982  // centrality bin
983  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
984  Int_t CenBin = -1;
985  CenBin = GetCenBin(centr);
986  if(CenBin==-1) return;
987 
988  // reconstructed
989  for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
990 
991  AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
992  if(!track) continue;
993 
994  // to select primaries
995  Int_t lp = TMath::Abs(track->GetLabel());
996 
997  // general kinematic & quality cuts
998  if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70) continue;
999 
1000  // Double_t DCAxy = track->DCA();
1001  // Double_t DCAz = track->ZAtDCA();
1002  // if (std::abs((Int_t)DCAxy)==999 || std::abs((Int_t)DCAz)==999) {
1003  // // re-evaluate the dca as it seems to not be natively present
1004  // // allowed only for tracks inside the beam pipe
1005  // Double_t pos[3] = {-99., -99., -99.};
1006  // track->GetPosition(pos);
1007  // if(pos[0]*pos[0]+pos[1]*pos[1] <= 3.*3.) {
1008  // AliAODTrack copy(*track); // stack copy
1009  // Double_t b[2] = {-99., -99.};
1010  // Double_t bCov[3] = {-99., -99., -99.};
1011  // if(copy.PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov)) {
1012  // DCAxy = b[0];
1013  // DCAz = b[1];
1014  // }
1015  // }
1016  // }
1017 
1018  // test filter bits
1019  if (((AliAODMCParticle*)fStack->At(lp))->IsPhysicalPrimary()) {
1020  if (track->TestFilterBit(32)) fPtSpecFB32[0][CenBin]->Fill(track->Pt());
1021  if (track->TestFilterBit(96)) fPtSpecFB96[0][CenBin]->Fill(track->Pt());
1022  if (track->TestFilterBit(128)) fPtSpecFB128[0][CenBin]->Fill(track->Pt());
1023  if (track->TestFilterBit(768)) fPtSpecFB768[0][CenBin]->Fill(track->Pt());
1024  } else {
1025  if (track->TestFilterBit(32)) fPtSpecFB32[1][CenBin]->Fill(track->Pt());
1026  if (track->TestFilterBit(96)) fPtSpecFB96[1][CenBin]->Fill(track->Pt());
1027  if (track->TestFilterBit(128)) fPtSpecFB128[1][CenBin]->Fill(track->Pt());
1028  if (track->TestFilterBit(768)) fPtSpecFB768[1][CenBin]->Fill(track->Pt());
1029  }
1030 
1031  fCenDis->Fill(centr);
1032  }
1033 
1034  // generated (physical primaries)
1035 
1036  for(Int_t jTracks = 0; jTracks<fStack->GetEntriesFast(); jTracks++) {
1037  AliAODMCParticle *MCpart = (AliAODMCParticle*)fStack->At(jTracks);
1038  if (!MCpart) {
1039  printf("ERROR: Could not receive MC track %d\n", jTracks);
1040  continue;
1041  }
1042  // kinematic cuts
1043  if ( MCpart->Pt() < .2 || MCpart->Pt() > 20. || TMath::Abs(MCpart->Eta()) > .8 ) continue;
1044  // select charged primaries
1045  if ( MCpart->Charge() == 0. || !MCpart->IsPhysicalPrimary()) continue;
1046 
1047  fPtSpecGen[0][CenBin]->Fill(MCpart->Pt());
1048  }
1049 
1050  // fGenHeader = McEvent->GenEventHeader();
1051  // if(fGenHeader) fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1052  // printf("#reconstructed : %d (rejected from cuts %d), #MC primaries : %d (rejected from cuts %d) \n",AODPOIs,AODbads,MCPrims,MCSecos);
1053  fFlowEvent->SetReferenceMultiplicity(aod->GetNumberOfTracks());
1054  fFlowEvent->SetCentrality(centr);
1055  // if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1056  fFlowEvent->SetRun(RunNum);
1057  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1058  }
1059 
1060  if(fAnalysisType == "MCESD") {
1061 
1062  fFlowEvent->ClearFast();
1063 
1064  if(!esd) {
1065  AliError("ERROR: Could not retrieve ESDEvent");
1066  return;
1067  }
1068  if(!McEvent) {
1069  AliError("ERROR: Could not retrieve MCEvent");
1070  return;
1071  }
1072  AliStack* fStack = fMCEvent->Stack();
1073  if(!fStack) {
1074  AliError("ERROR: Could not retrieve MCStack");
1075  return;
1076  }
1077 
1078  AliESDVertex *vertex = (AliESDVertex*) esd->GetPrimaryVertex();
1079  if (!vertex) return;
1080  if (TMath::Abs(vertex->GetZ()) > 10. ) return;
1081  if (vertex->GetNContributors() < 1 ) return;
1082  AliCentrality *centrality = esd->GetCentrality();
1083  if (!centrality) return;
1084  Double_t centr = centrality->GetCentralityPercentile("V0M");
1085  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1086  Int_t CenBin = -1;
1087  if (centr>0. && centr<5.) CenBin=0;
1088  if (centr>5. && centr<10.) CenBin=1;
1089  if (centr>10. && centr<20.) CenBin=2;
1090  if (centr>20. && centr<30.) CenBin=3;
1091  if (centr>30. && centr<40.) CenBin=4;
1092  if (centr>40. && centr<50.) CenBin=5;
1093  if (centr>50. && centr<60.) CenBin=6;
1094  if (centr>60. && centr<70.) CenBin=7;
1095  if (centr>70. && centr<80.) CenBin=8;
1096  if (centr>80. && centr<90.) CenBin=9;
1097  if(CenBin==-1) return;
1098 
1099  //Generated
1100  Int_t MCPrims = 0;
1101  for ( Int_t i=0 ; i<fStack->GetNtrack() ; i++ ) {
1102 
1103  //Primaries Selection
1104  TParticle *particle = (TParticle*)fStack->Particle(i);
1105  if (!particle) continue;
1106  if (!fStack->IsPhysicalPrimary(i)) continue;
1107  if ( particle->GetPDG()->Charge() == 0.) continue;
1108 
1109  //Kinematic Cuts
1110  if ( particle->Pt()<0.2 || particle->Pt()>10. ) continue;
1111  if ( TMath::Abs(particle->Eta())>0.8 ) continue;
1112 
1113  fFlowTrack->SetPhi(particle->Phi());
1114  fFlowTrack->SetEta(particle->Eta());
1115  fFlowTrack->SetPt(particle->Pt());
1117  fFlowTrack->SetForRPSelection(kTRUE);
1119  fFlowTrack->SetForPOISelection(kFALSE);
1121  MCPrims++;
1122 
1123  fPtSpecGen[0][CenBin]->Fill(particle->Pt());
1124 
1125  }
1126 
1127  //Reconstructed
1128  Int_t ESDPrims = 0;
1129  for (Int_t i=0 ; i<esd->GetNumberOfTracks() ; i++) {
1130 
1131  //Get reconstructed track
1132  AliVTrack *vtrack = static_cast<AliVTrack*>(esd->GetTrack(i));
1133  AliESDtrack *track = dynamic_cast<AliESDtrack*>(vtrack);
1134  if (!track) continue;
1135 
1136  //Primaries selection
1137  Int_t lp = TMath::Abs(track->GetLabel());
1138  if (!fStack->IsPhysicalPrimary(lp)) continue;
1139  TParticle *particle = (TParticle*)fStack->Particle(lp);
1140  if (!particle) continue;
1141  if (particle->GetPDG()->Charge() == 0.) continue;
1142 
1143  // if(!fCutsPOI->PassesESDcuts(track)) continue;
1144 
1145  Bool_t pass = kTRUE;
1146 
1147  if(fCutTPC) {
1148  // printf("******* cutting TPC ******** \n");
1149  UShort_t ntpccls = track->GetTPCNcls();
1150  Double_t tpcchi2 = track->GetTPCchi2();
1151  if (tpcchi2<0.2 || tpcchi2 >=4.) {
1152  // printf("TPCchi2 : %e %e ",tpcchi2,track->GetTPCchi2Iter1());
1153  pass=kFALSE;
1154  }
1155  if (ntpccls < 70) {
1156  // printf("#TPCcluster : %u %u %u %u ",ntpccls,track->GetTPCNclsF(),track->GetTPCNclsFIter1(),track->GetTPCNclsIter1());
1157  pass=kFALSE;
1158  }
1159  }
1160 
1161  Float_t dcaxy=0.0;
1162  Float_t dcaz=0.0;
1163  track->GetImpactParameters(dcaxy,dcaz);
1164  if (dcaxy > 0.3 || dcaz > 0.3) {
1165  // printf("DCA : %e %e ",dcaxy,dcaz);
1166  pass=kFALSE;
1167  }
1168  if(!pass) continue;
1169 
1170  //Kinematic Cuts
1171  if ( track->Pt()<0.2 || track->Pt()>10. ) continue;
1172  if ( TMath::Abs(track->Eta())>0.8 ) continue;
1173 
1174  fFlowTrack->SetPhi(track->Phi());
1175  fFlowTrack->SetEta(track->Eta());
1176  fFlowTrack->SetPt(track->Pt());
1178  fFlowTrack->SetForRPSelection(kFALSE);
1182  ESDPrims++;
1183 
1184  }
1185 
1186  // printf("#reconstructed : %d , #MC primaries : %d \n",ESDPrims,MCPrims);
1187  fFlowEvent->SetReferenceMultiplicity(esd->GetNumberOfTracks());
1188  fFlowEvent->SetCentrality(centr);
1189  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1190  fFlowEvent->SetRun(esd->GetRunNumber());
1191  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1192 
1193  } // end of if(fAnalysisType == "MCESD")
1194 
1195  if(fAnalysisType == "MCkine") {
1196 
1197  fFlowEvent->ClearFast();
1198 
1199  AliInputEventHandler* McHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1200  if(!McHandler) {
1201  AliError("ERROR: Could not retrieve MCtruthEventHandler");
1202  return;
1203  }
1204  McEvent = McHandler->MCEvent();
1205  if(!McEvent) {
1206  AliError("ERROR: Could not retrieve MC event");
1207  return;
1208  }
1209 
1210  Int_t nTracks = McEvent->GetNumberOfTracks();
1211  // Int_t nPrimTr = McEvent->GetNumberOfPrimaries();
1212 
1213  //loop over tracks
1214  for (Int_t itrkN=0; itrkN<nTracks; itrkN++) {
1215  //get input particle
1216  AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(McEvent->GetTrack(itrkN));
1217  if (!pParticle) continue;
1218 
1219  //check if track passes the cuts
1220  if (McEvent->IsPhysicalPrimary(itrkN) && pParticle->Charge()!=0) {
1221  fFlowTrack->Set(pParticle);
1223  fFlowTrack->SetForRPSelection(kTRUE);
1228  }
1229  }// for all tracks
1230 
1231  // if monte carlo event get reaction plane from monte carlo (depends on generator)
1232  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1233  // set reference multiplicity
1234  fFlowEvent->SetReferenceMultiplicity(McEvent->GetNumberOfTracks());
1235  // tag subevents
1237  // set centrality from impact parameter
1238  Double_t ImpPar=0., CenPer=0.;
1239  fGenHeader = McEvent->GenEventHeader();
1240  if(fGenHeader){
1241  fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1242  if(fPythiaGenHeader) ImpPar = fPythiaGenHeader->GetImpactParameter();
1243  fHijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(fGenHeader);
1244  if(fHijingGenHeader) ImpPar = fHijingGenHeader->ImpactParameter();
1245  if(ImpPar) CenPer = 0.3859796743103508*pow(ImpPar,2.);
1246  if(CenPer>0. && CenPer<100.) fFlowEvent->SetCentrality(CenPer);
1247  else return;
1248  fFlowEvent->SetRun(1);
1249  }
1250 
1251  } // end of if(fAnalysisType == "MCkine")
1252 
1253  if (!fFlowEvent) return; //shuts up coverity
1254 
1255  //check final event cuts
1256  Int_t mult = fFlowEvent->NumberOfTracks();
1257  // AliInfo(Form("FlowEvent has %i tracks",mult));
1258  if (mult<fMinMult || mult>fMaxMult) {
1259  AliWarning("FlowEvent cut on multiplicity"); return;
1260  }
1261 
1262  //define dead zone
1264 
1267  if (fAfterburnerOn)
1268  {
1269  //if reaction plane not set from elsewhere randomize it before adding flow
1271  fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
1272 
1273  if(fDifferentialV2)
1275  else
1276  fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
1277  fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
1278  }
1280 
1281  //tag subEvents
1283 
1284  //do we want to serve shullfed tracks to everybody?
1286 
1287  // associate the mother particles to their daughters in the flow event (if any)
1289 
1290  //fListHistos->Print();
1291  //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
1292 
1293  //********************************************************************************************************************************
1294 
1295  if(fAnalysisType == "AOD" || fAnalysisType == "AUTOMATIC") {
1296 
1297  // PHYSICS SELECTION
1298  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1299  AliInputEventHandler *hdr = (AliInputEventHandler*)am->GetInputEventHandler();
1300 
1301  if(hdr->IsEventSelected() && AliVEvent::kAny) {
1302 
1303  Double_t centrperc = fFlowEvent->GetCentrality();
1304  Int_t cenb = (Int_t)centrperc;
1305 
1306  AliAODTracklets *trackl = aod->GetTracklets();
1307  Int_t nTracklets = trackl->GetNumberOfTracklets();
1308 
1309  AliAODVZERO *vzeroAOD = aod->GetVZEROData();
1310  Double_t multV0A = vzeroAOD->GetMTotV0A();
1311  Double_t multV0C = vzeroAOD->GetMTotV0C();
1312 
1313  AliAODZDC *aodZDC = aod->GetZDCData();
1314 
1315  Double_t energyZNC = (Double_t) (aodZDC->GetZNCEnergy());
1316  Double_t energyZPC = (Double_t) (aodZDC->GetZPCEnergy());
1317  Double_t energyZNA = (Double_t) (aodZDC->GetZNAEnergy());
1318  Double_t energyZPA = (Double_t) (aodZDC->GetZPAEnergy());
1319  Double_t energyZEM1 = (Double_t) (aodZDC->GetZEM1Energy());
1320  Double_t energyZEM2 = (Double_t) (aodZDC->GetZEM2Energy());
1321 
1322  const Double_t * towZNCraw = aodZDC->GetZNCTowerEnergy();
1323  const Double_t * towZNAraw = aodZDC->GetZNATowerEnergy();
1324 
1325  // Get centroid from ZDCs *******************************************************
1326 
1327  Double_t Enucl = (RunNum < 209122 ? 1380. : 2511.);
1328  Double_t xyZNC[2]={999.,999.}, xyZNA[2]={999.,999.};
1329  Double_t towZNC[5]={0.}, towZNA[5]={0.};
1330 
1331  Double_t ZNCcalib=1., ZNAcalib=1.;
1332  if(fUseTowerEq) {
1333  if(RunNum!=fCachedRunNum) {
1334  for(Int_t i=0; i<5; i++) {
1335  fTowerGainEq[0][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNCTower[%d][%d]",RunNum,i)));
1336  fTowerGainEq[1][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNATower[%d][%d]",RunNum,i)));
1337  }
1338  }
1339  for(Int_t i=0; i<5; i++) {
1340  if(fTowerGainEq[0][i]) towZNC[i] = towZNCraw[i]*fTowerGainEq[0][i]->GetBinContent(fTowerGainEq[0][i]->FindBin(centrperc));
1341  if(fTowerGainEq[1][i]) towZNA[i] = towZNAraw[i]*fTowerGainEq[1][i]->GetBinContent(fTowerGainEq[1][i]->FindBin(centrperc));
1342  }
1343  } else {
1344  for(Int_t i=0; i<5; i++) {
1345  towZNC[i] = towZNCraw[i];
1346  towZNA[i] = towZNAraw[i];
1347  fZNCTower[RunBin][i]->Fill(centrperc,towZNC[i]);
1348  fZNATower[RunBin][i]->Fill(centrperc,towZNA[i]);
1349  }
1350  }
1351 
1352  if(RunNum>=245829) towZNA[2] = 0.;
1353  Double_t zncEnergy=0., znaEnergy=0.;
1354  for(Int_t i=0; i<5; i++){
1355  zncEnergy += towZNC[i];
1356  znaEnergy += towZNA[i];
1357  }
1358  if(RunNum>=245829) znaEnergy *= 8./7.;
1359  fFlowEvent->SetZNCEnergy(towZNC[0]);
1360  fFlowEvent->SetZNAEnergy(towZNA[0]);
1361 
1362  const Double_t x[4] = {-1.75, 1.75, -1.75, 1.75};
1363  const Double_t y[4] = {-1.75, -1.75, 1.75, 1.75};
1364  Double_t numXZNC=0., numYZNC=0., denZNC=0., cZNC, wZNC;
1365  Double_t numXZNA=0., numYZNA=0., denZNA=0., cZNA, wZNA, BadChOr;
1366 
1367  if (fUseMCCen) {
1368  for(Int_t i=0; i<4; i++){
1369  wZNC = TMath::Power(towZNC[i+1], fZDCGainAlpha);
1370  numXZNC += x[i]*wZNC;
1371  numYZNC += y[i]*wZNC;
1372  denZNC += wZNC;
1373  fhZNSpectra->Fill(centrperc,i+0.5,towZNC[i+1]);
1374 
1375  if(i==1) {
1376  wZNA = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1377  if(fBadTowerCalibHist[cenb]) {
1378  wZNA = GetBadTowerResp(wZNA, fBadTowerCalibHist[cenb]);
1379  }
1380  BadChOr = wZNA;
1381  wZNA = TMath::Power(wZNA, fZDCGainAlpha);
1382  }
1383  else wZNA = TMath::Power(towZNA[i+1], fZDCGainAlpha);
1384  numXZNA += x[i]*wZNA;
1385  numYZNA += y[i]*wZNA;
1386  denZNA += wZNA;
1387  fhZNSpectra->Fill(centrperc,i+4.5,(i!=1?towZNA[i+1]:BadChOr));
1388  }
1389  // store distribution for unfolding
1390  if(RunNum<245829) {
1391  Double_t recoE = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1392  Double_t trueE = towZNA[2];
1393  fhZNBCCorr->Fill(centrperc,trueE,recoE);
1394  }
1395  if(denZNC!=0){
1396  Double_t nSpecnC = zncEnergy/Enucl;
1397  cZNC = 1.89358-0.71262/(nSpecnC+0.71789);
1398  xyZNC[0] = cZNC*numXZNC/denZNC;
1399  xyZNC[1] = cZNC*numYZNC/denZNC;
1400  }
1401  else{
1402  xyZNC[0] = xyZNC[1] = 999.;
1403  }
1404  if(denZNA!=0){
1405  Double_t nSpecnA = znaEnergy/Enucl;
1406  cZNA = 1.89358-0.71262/(nSpecnA+0.71789);
1407  xyZNA[0] = cZNA*numXZNA/denZNA;
1408  xyZNA[1] = cZNA*numYZNA/denZNA;
1409  }
1410  else{
1411  xyZNA[0] = xyZNA[1] = 999.;
1412  }
1413  } else {
1414  for(Int_t i=0; i<4; i++) {
1415  if(towZNC[i+1]>0.) {
1416  wZNC = TMath::Power(towZNC[i+1], fZDCGainAlpha);
1417  numXZNC += x[i]*wZNC;
1418  numYZNC += y[i]*wZNC;
1419  denZNC += wZNC;
1420  }
1421  if(towZNA[i+1]>0.) {
1422  wZNA = TMath::Power(towZNA[i+1], fZDCGainAlpha);
1423  numXZNA += x[i]*wZNA;
1424  numYZNA += y[i]*wZNA;
1425  denZNA += wZNA;
1426  }
1427  }
1428  if(denZNC!=0) {
1429  xyZNC[0] = numXZNC/denZNC;
1430  xyZNC[1] = numYZNC/denZNC;
1431  }
1432  else{
1433  xyZNC[0] = xyZNC[1] = 999.;
1434  zncEnergy = 0.;
1435  }
1436  if(denZNA!=0) {
1437  xyZNA[0] = numXZNA/denZNA;
1438  xyZNA[1] = numYZNA/denZNA;
1439  }
1440  else{
1441  xyZNA[0] = xyZNA[1] = 999.;
1442  znaEnergy = 0.;
1443  }
1444  }
1445 
1446  fhZNCcentroid->Fill(xyZNC[0], xyZNC[1]);
1447  fhZNAcentroid->Fill(xyZNA[0], xyZNA[1]);
1448  fFlowEvent->SetZDC2Qsub(xyZNC,denZNC,xyZNA,denZNA);
1449 
1450  // ******************************************************************************
1451 
1452  Double_t tdcSum = aodZDC->GetZDCTimeSum();
1453  Double_t tdcDiff = aodZDC->GetZDCTimeDiff();
1454  fhDebunch->Fill(tdcDiff, tdcSum);
1455 
1456  for(int i=0; i<5; i++){
1457  fhZNCPM[i]->Fill(towZNC[i]);
1458  if((i<4) && (towZNC[0]>0.)) fhZNCPMQiPMC[i]->Fill(towZNC[i+1]/towZNC[0]);
1459  }
1460  for(int i=0; i<5; i++){
1461  fhZNAPM[i]->Fill(towZNA[i]);
1462  if(((i<4) && towZNA[0]>0.)) fhZNAPMQiPMC[i]->Fill(towZNA[i+1]/towZNA[0]);
1463  }
1464 
1465  fhZNCvsZNA->Fill(energyZNA, energyZNC);
1466  fhZDCCvsZDCCA->Fill(energyZNA+energyZPA, energyZNC+energyZPC);
1467  fhZNCvsZPC->Fill(energyZPC, energyZNC);
1468  fhZNAvsZPA->Fill(energyZPA, energyZNA);
1469  fhZNvsZP->Fill(energyZPA+energyZPC, energyZNA+energyZNC);
1470  fhZNvsVZERO->Fill(multV0A+multV0C, energyZNC+energyZNA);
1471  fhZDCvsVZERO->Fill(multV0A+multV0C, energyZNA+energyZPA+energyZNC+energyZPC);
1472  fhZDCvsTracklets->Fill((Double_t) (nTracklets), energyZNA+energyZPA+energyZNC+energyZPC);
1473 
1474  Double_t asymmetry = -999.;
1475  if((energyZNC+energyZNA)>0.) asymmetry = (energyZNC-energyZNA)/(energyZNC+energyZNA);
1476  fhAsymm->Fill(asymmetry);
1477  fhZNAvsAsymm->Fill(asymmetry, energyZNA/1000.);
1478  fhZNCvsAsymm->Fill(asymmetry, energyZNC/1000.);
1479 
1480  fhZNCvscentrality->Fill(centrperc, energyZNC/1000.);
1481  fhZNAvscentrality->Fill(centrperc, energyZNA/1000.);
1482 
1483  } // PHYSICS SELECTION
1484 
1485  }
1486 
1487  // p) cache run number
1489 
1490  // printf("debug: NoRPs %e, NoPOIs %e, RunNum %d, Cen %e \n",fFlowEvent->GetNumberOfRPs(),fFlowEvent->GetNumberOfPOIs(),fCachedRunNum,fFlowEvent->GetCentrality());
1491 
1492  PostData(1, fFlowEvent);
1493 
1494  PostData(2, fOutput);
1495 }
1496 //________________________________________________________________________
1497 
1499 {
1500  Double_t EtC = BadTowerCalibHist->ProjectionY("",BadTowerCalibHist->GetXaxis()->FindBin(Et),BadTowerCalibHist->GetXaxis()->FindBin(Et))->GetRandom();
1501  return EtC;
1502 }
1503 
1504 //________________________________________________________________________
1505 
1507 {
1508  Int_t CenBin=-1;
1509  if (Centrality>0. && Centrality<5.) CenBin=0;
1510  if (Centrality>5. && Centrality<10.) CenBin=1;
1511  if (Centrality>10. && Centrality<20.) CenBin=2;
1512  if (Centrality>20. && Centrality<30.) CenBin=3;
1513  if (Centrality>30. && Centrality<40.) CenBin=4;
1514  if (Centrality>40. && Centrality<50.) CenBin=5;
1515  if (Centrality>50. && Centrality<60.) CenBin=6;
1516  if (Centrality>60. && Centrality<70.) CenBin=7;
1517  if (Centrality>70. && Centrality<80.) CenBin=8;
1518  if (Centrality>80. && Centrality<90.) CenBin=9;
1519  if (CenBin>=fnCen) CenBin=-1;
1520  if (fnCen==1) CenBin=0;
1521  return CenBin;
1522 } // end of AliFlowAnalysisCRC::GetCRCCenBin(Double_t Centrality)
1523 //_____________________________________________________________________________
1524 
1525 Double_t AliAnalysisTaskCRCZDC::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
1526 {
1527  // calculate sqrt of weighted distance to other vertex
1528  if (!v0 || !v1) {
1529  printf("One of vertices is not valid\n");
1530  return 0;
1531  }
1532  static TMatrixDSym vVb(3);
1533  double dist = -1;
1534  double dx = v0->GetX()-v1->GetX();
1535  double dy = v0->GetY()-v1->GetY();
1536  double dz = v0->GetZ()-v1->GetZ();
1537  double cov0[6],cov1[6];
1538  v0->GetCovarianceMatrix(cov0);
1539  v1->GetCovarianceMatrix(cov1);
1540  vVb(0,0) = cov0[0]+cov1[0];
1541  vVb(1,1) = cov0[2]+cov1[2];
1542  vVb(2,2) = cov0[5]+cov1[5];
1543  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1544  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1545  vVb.InvertFast();
1546  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
1547  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1548  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1549  return dist>0 ? TMath::Sqrt(dist) : -1;
1550 }
1551 //________________________________________________________________________
1552 
1554 {
1555  // check for multi-vertexer pile-up
1556 
1557  const int kMinPlpContrib = 5;
1558  const double kMaxPlpChi2 = 5.0;
1559  const double kMinWDist = 15;
1560 
1561  const AliVVertex* vtPrm = 0;
1562  const AliVVertex* vtPlp = 0;
1563  int nPlp = 0;
1564 
1565  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
1566  vtPrm = aod->GetPrimaryVertex();
1567  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
1568 
1569  //int bcPrim = vtPrm->GetBC();
1570 
1571  for (int ipl=0;ipl<nPlp;ipl++) {
1572  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
1573  //
1574  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1575  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1576  // int bcPlp = vtPlp->GetBC();
1577  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
1578  //
1579  double wDst = GetWDist(vtPrm,vtPlp);
1580  if (wDst<kMinWDist) continue;
1581  //
1582  return kTRUE; // pile-up: well separated vertices
1583  }
1584 
1585  return kFALSE;
1586 }
1587 
1588 //________________________________________________________________________
1590  fCutContainer->AddAt(cutsRP,0); fCutsRP=cutsRP; cutsRP->SetPOItype(0);
1591 }
1592 
1593 //________________________________________________________________________
1595  fCutContainer->AddAt(cutsPOI,1); fCutsPOI=cutsPOI; cutsPOI->SetPOItype(1);
1596 }
1597 
1598 //________________________________________________________________________
1600 {
1601  // Terminate analysis
1602  //
1603  /* if(fDebug > 1) printf(" **** AliAnalysisTaskCRCZDC::Terminate() \n");
1604 
1605  //fOutput = dynamic_cast<TList*> (GetOutputData(1));
1606  //if(!fOutput) printf("ERROR: fOutput not available\n");
1607  */
1608 }
TH2F * fhZDCCvsZDCCA
ZNC vs ZNA;.
TH2F * fhZNvsVZERO
ZNC+ZNA vs ZPC+ZPA;.
TH1F * fPtSpecFB128[2][10]
PtSpecRec FB96.
void FindDaughters(Bool_t keepDaughtersInRPselection=kFALSE)
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)
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
ZNA centroid.
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)
TH2F * fhZNAcentroid
ZNC centroid.
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
TList * fBadTowerCalibList
list for storing calib files
TH1F * fPtSpecFB96[2][10]
PtSpecRec FB32.
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()
AliAnalysisUtils * fAnalysisUtil
Event selection.
TH2F * fhZNCcentroid
Debunch;.
TH2F * fhZNvsZP
ZNA vs ZPA;.
AliGenHijingEventHeader * fHijingGenHeader
const Double_t phimin
AliGenPythiaEventHeader * fPythiaGenHeader
Int_t NumberOfTracks() const