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