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