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