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