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