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