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