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