AliPhysics  449db5a (449db5a)
 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(kAUTOMATIC),
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 fCorrectPhiTracklets(kFALSE),
204 fZDCSpectraCorrList(NULL),
205 fSpectraMCList(NULL),
206 fTrackQAList(NULL),
207 fBadTowerStuffList(NULL),
208 fVZEROStuffList(NULL),
209 fVZEROGainEqHist(NULL),
210 fMinRingVZC(1),
211 fMaxRingVZC(4),
212 fMinRingVZA(5),
213 fMaxRingVZA(8),
214 fCachedRunNum(0),
215 fhZNSpectra(0x0),
216 fhZNSpectraCor(0x0),
217 fhZNSpectraPow(0x0),
218 fhZNBCCorr(0x0)
219 {
220  for(int i=0; i<5; i++){
221  fhZNCPM[i] = 0x0;
222  fhZNAPM[i] = 0x0;
223  }
224  for(int i=0; i<4; i++){
225  fhZNCPMQiPMC[i] = 0x0;
226  fhZNAPMQiPMC[i] = 0x0;
227  }
228  for(Int_t r=0; r<fCRCMaxnRun; r++) {
229  fRunList[r] = 0;
230  }
231  for(Int_t c=0; c<2; c++) {
232  for(Int_t i=0; i<5; i++) {
233  fTowerGainEq[c][i] = NULL;
234  }
235  }
236  for(Int_t c=0; c<100; c++) {
237  fBadTowerCalibHist[c] = NULL;
238  }
239  for (Int_t k=0; k<fkVZEROnHar; k++) {
240 // fVZEROQVectorRecQx[k] = NULL;
241 // fVZEROQVectorRecQy[k] = NULL;
242  fVZEROQVectorRecQxStored[k] = NULL;
243  fVZEROQVectorRecQyStored[k] = NULL;
244  for (Int_t t=0; t<fkVZEROnQAplots; t++) {
245  fVZEROQVectorRecFinal[k][t] = NULL;
246  }
247  }
248  for(Int_t i=0; i<8; i++) {
249  SpecCorMu1[i] = NULL;
250  SpecCorMu2[i] = NULL;
251  SpecCorSi[i] = NULL;
252  SpecCorAv[i] = NULL;
253  }
254  this->InitializeRunArrays();
255  fMyTRandom3 = new TRandom3(1);
256  gRandom->SetSeed(fMyTRandom3->Integer(65539));
257  for(Int_t j=0; j<2; j++) {
258  for(Int_t c=0; c<10; c++) {
259  fPtSpecGen[j][c] = NULL;
260  fPtSpecFB32[j][c] = NULL;
261  fPtSpecFB96[j][c] = NULL;
262  fPtSpecFB128[j][c] = NULL;
263  fPtSpecFB768[j][c] = NULL;
264  }
265  }
266  for (Int_t c=0; c<2; c++) {
267  fhZNCenDis[c] = NULL;
268  }
269  for(Int_t fb=0; fb<fKNFBs; fb++){
270  for (Int_t i=0; i<4; i++) {
271  fTrackQADCAxy[fb][i] = NULL;
272  fTrackQADCAz[fb][i] = NULL;
273  fTrackQApT[fb][i] = NULL;
274  fTrackQADphi[fb][i] = NULL;
275  fEbEQRe[fb][i] = NULL;
276  fEbEQIm[fb][i] = NULL;
277  fEbEQMu[fb][i] = NULL;
278  }
279  }
280 }
281 
282 //________________________________________________________________________
283 AliAnalysisTaskCRCZDC::AliAnalysisTaskCRCZDC(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates):
284 AliAnalysisTaskSE(name),
285 fAnalysisType(kAUTOMATIC),
286 fRPType(RPtype),
287 fCFManager1(NULL),
288 fCFManager2(NULL),
289 fCutsEvent(NULL),
290 fCutsRP(NULL),
291 fCutsPOI(NULL),
292 fCutContainer(new TList()),
293 fQAList(NULL),
294 fMinMult(0),
295 fMaxMult(10000000),
296 fMinA(-1.0),
297 fMaxA(-0.01),
298 fMinB(0.01),
299 fMaxB(1.0),
300 fQAon(on),
301 fLoadCandidates(bCandidates),
302 fNbinsMult(10000),
303 fNbinsPt(100),
304 fNbinsPhi(100),
305 fNbinsEta(200),
306 fNbinsQ(500),
307 fNbinsMass(1),
308 fMultMin(0.),
309 fMultMax(10000.),
310 fPtMin(0.),
311 fPtMax(10.),
312 fPhiMin(0.),
313 fPhiMax(TMath::TwoPi()),
314 fEtaMin(-5.),
315 fEtaMax(5.),
316 fQMin(0.),
317 fQMax(3.),
318 fMassMin(-1.),
319 fMassMax(0.),
320 fHistWeightvsPhiMin(0.),
321 fHistWeightvsPhiMax(3.),
322 fExcludedEtaMin(0.),
323 fExcludedEtaMax(0.),
324 fExcludedPhiMin(0.),
325 fExcludedPhiMax(0.),
326 fAfterburnerOn(kFALSE),
327 fNonFlowNumberOfTrackClones(0),
328 fV1(0.),
329 fV2(0.),
330 fV3(0.),
331 fV4(0.),
332 fV5(0.),
333 fDifferentialV2(0),
334 fFlowEvent(NULL),
335 fShuffleTracks(kFALSE),
336 fMyTRandom3(NULL),
337 fAnalysisInput(kAOD),
338 fIsMCInput(kFALSE),
339 fUseMCCen(kTRUE),
340 fRejectPileUp(kTRUE),
341 fRejectPileUpTight(kFALSE),
342 fResetNegativeZDC(kFALSE),
343 fCentrLowLim(0.),
344 fCentrUpLim(100.),
345 fCentrEstimator(kV0M),
346 fOutput(0x0),
347 fhZNCvsZNA(0x0),
348 fhZDCCvsZDCCA(0x0),
349 fhZNCvsZPC(0x0),
350 fhZNAvsZPA(0x0),
351 fhZNvsZP(0x0),
352 fhZNvsVZERO(0x0),
353 fhZDCvsVZERO(0x0),
354 fhZDCvsTracklets(0x0),
355 fhZDCvsNclu1(0x0),
356 fhDebunch(0x0),
357 fhAsymm(0x0),
358 fhZNAvsAsymm(0x0),
359 fhZNCvsAsymm(0x0),
360 fhZNCvscentrality(0x0),
361 fhZNAvscentrality(0x0),
362 fDataSet(kAny),
363 fCRCnRun(0),
364 fZDCGainAlpha(0.395),
365 fGenHeader(NULL),
366 fPythiaGenHeader(NULL),
367 fHijingGenHeader(NULL),
368 fFlowTrack(NULL),
369 fAnalysisUtil(NULL),
370 fStack(0x0),
371 fCutTPC(kFALSE),
372 fCenDis(0x0),
373 fVZEROMult(0x0),
374 fMultSelection(0x0),
375 fPileUpCount(0x0),
376 fPileUpMultSelCount(0x0),
377 fMultTOFLowCut(0x0),
378 fMultTOFHighCut(0x0),
379 fUseTowerEq(kFALSE),
380 fTowerEqList(NULL),
381 fUseBadTowerCalib(kFALSE),
382 fBadTowerCalibList(NULL),
383 fVZEROGainEqList(NULL),
384 fVZEROQVecRecList(NULL),
385 fUseZDCSpectraCorr(kFALSE),
386 fCorrectPhiTracklets(kFALSE),
387 fZDCSpectraCorrList(NULL),
388 fSpectraMCList(NULL),
389 fTrackQAList(NULL),
390 fBadTowerStuffList(NULL),
391 fVZEROStuffList(NULL),
392 fVZEROGainEqHist(NULL),
393 fMinRingVZC(1),
394 fMaxRingVZC(4),
395 fMinRingVZA(5),
396 fMaxRingVZA(8),
397 fCachedRunNum(0),
398 fhZNSpectra(0x0),
399 fhZNSpectraCor(0x0),
400 fhZNSpectraPow(0x0),
401 fhZNBCCorr(0x0)
402 {
403 
404  for(int i=0; i<5; i++){
405  fhZNCPM[i] = 0x0;
406  fhZNAPM[i] = 0x0;
407  }
408  for(int i=0; i<4; i++){
409  fhZNCPMQiPMC[i] = 0x0;
410  fhZNAPMQiPMC[i] = 0x0;
411  }
412  for(Int_t r=0; r<fCRCMaxnRun; r++) {
413  fRunList[r] = 0;
414  }
415  for(Int_t c=0; c<2; c++) {
416  for(Int_t i=0; i<5; i++) {
417  fTowerGainEq[c][i] = NULL;
418  }
419  }
420  for(Int_t c=0; c<100; c++) {
421  fBadTowerCalibHist[c] = NULL;
422  }
423  for (Int_t k=0; k<fkVZEROnHar; k++) {
424  // fVZEROQVectorRecQx[k] = NULL;
425  // fVZEROQVectorRecQy[k] = NULL;
426  fVZEROQVectorRecQxStored[k] = NULL;
427  fVZEROQVectorRecQyStored[k] = NULL;
428  for (Int_t t=0; t<fkVZEROnQAplots; t++) {
429  fVZEROQVectorRecFinal[k][t] = NULL;
430  }
431  }
432  for(Int_t i=0; i<8; i++) {
433  SpecCorMu1[i] = NULL;
434  SpecCorMu2[i] = NULL;
435  SpecCorSi[i] = NULL;
436  SpecCorAv[i] = NULL;
437  }
438  this->InitializeRunArrays();
439  fMyTRandom3 = new TRandom3(iseed);
440  gRandom->SetSeed(fMyTRandom3->Integer(65539));
441 
442  DefineInput(0, TChain::Class());
443  // Define output slots here
444  // Define here the flow event output
445  DefineOutput(1, AliFlowEventSimple::Class());
446  DefineOutput(2, TList::Class());
447 
448  for(Int_t j=0; j<2; j++) {
449  for(Int_t c=0; c<10; c++) {
450  fPtSpecGen[j][c] = NULL;
451  fPtSpecFB32[j][c] = NULL;
452  fPtSpecFB96[j][c] = NULL;
453  fPtSpecFB128[j][c] = NULL;
454  fPtSpecFB768[j][c] = NULL;
455  }
456  }
457  for (Int_t c=0; c<2; c++) {
458  fhZNCenDis[c] = NULL;
459  }
460  for(Int_t fb=0; fb<fKNFBs; fb++){
461  for (Int_t i=0; i<4; i++) {
462  fTrackQADCAxy[fb][i] = NULL;
463  fTrackQADCAz[fb][i] = NULL;
464  fTrackQApT[fb][i] = NULL;
465  fTrackQADphi[fb][i] = NULL;
466  fEbEQRe[fb][i] = NULL;
467  fEbEQIm[fb][i] = NULL;
468  fEbEQMu[fb][i] = NULL;
469  }
470  }
471 }
472 
473 //________________________________________________________________________
475 {
476  // Destructor
477  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
478  delete fOutput; fOutput=0;
479  }
480  delete fMyTRandom3;
481  delete fFlowEvent;
482  delete fFlowTrack;
483  delete fCutsEvent;
484  if (fTowerEqList) delete fTowerEqList;
489  if (fAnalysisUtil) delete fAnalysisUtil;
490  if (fQAList) delete fQAList;
491  if (fCutContainer) fCutContainer->Delete(); delete fCutContainer;
492 }
493 
494 //________________________________________________________________________
496 {
497  for(Int_t r=0;r<fCRCMaxnRun;r++) {
498  fCRCQVecListRun[r] = NULL;
499  for(Int_t k=0;k<fCRCnTow;k++) {
500  fZNCTower[r][k] = NULL;
501  fZNATower[r][k] = NULL;
502  }
503 // fhZNSpectraRbR[r] = NULL;
504  }
505  // for(Int_t i=0;i<fnCen;i++) {
506  // fPtPhiEtaRbRFB128[r][i] = NULL;
507  // fPtPhiEtaRbRFB768[r][i] = NULL;
508  // }
509  // }
510 }
511 
512 //________________________________________________________________________
514 {
515  // Create the output containers
516 
517  //set the common constants
520  cc->SetNbinsPt(fNbinsPt);
521  cc->SetNbinsPhi(fNbinsPhi);
522  cc->SetNbinsEta(fNbinsEta);
523  cc->SetNbinsQ(fNbinsQ);
525  cc->SetMultMin(fMultMin);
526  cc->SetMultMax(fMultMax);
527  cc->SetPtMin(fPtMin);
528  cc->SetPtMax(fPtMax);
529  cc->SetPhiMin(fPhiMin);
530  cc->SetPhiMax(fPhiMax);
531  cc->SetEtaMin(fEtaMin);
532  cc->SetEtaMax(fEtaMax);
533  cc->SetQMin(fQMin);
534  cc->SetQMax(fQMax);
535  cc->SetMassMin(fMassMin);
536  cc->SetMassMax(fMassMax);
539 
540  fFlowEvent = new AliFlowEvent(20000);
541  fFlowTrack = new AliFlowTrack();
542 
543  //printf(" AliAnalysisTaskCRCZDC::UserCreateOutputObjects()\n\n");
544  fOutput = new TList();
545  fOutput->SetOwner(kTRUE);
546  //fOutput->SetName("output");
547 
548  if (fQAon) {
549  fQAList = new TList();
550  fQAList->SetOwner(kTRUE);
551  fQAList->SetName("AliFlowEventCuts QA");
552  if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA()); //0
553  if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA()); //1
554  if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA()); //2
555  fOutput->Add(fQAList);
556  }
557 
558  fVZEROStuffList = new TList();
559  fVZEROStuffList->SetOwner(kTRUE);
560  fVZEROStuffList->SetName("VZERO stuff");
561  fOutput->Add(fVZEROStuffList);
562 
563  fBadTowerStuffList = new TList();
564  fBadTowerStuffList->SetOwner(kTRUE);
565  fBadTowerStuffList->SetName("BadTowerCalib");
567 
568  fCenDis = new TH1F("fCenDis", "fCenDis", 100, 0., 100.);
569  fOutput->Add(fCenDis);
570  fPileUpCount = new TH1F("fPileUpCount", "fPileUpCount", 9, 0., 9.);
571  fPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
572  fPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
573  fPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultiplicityComb08");
574  fPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
575  fPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>7.5");
576  fPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
577  fPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
578  fPileUpCount->GetXaxis()->SetBinLabel(8,"multESDTPCDif");
579  fPileUpCount->GetXaxis()->SetBinLabel(9,"extraPileUpMultSel");
580  fOutput->Add(fPileUpCount);
581  fPileUpMultSelCount = new TH1F("fPileUpMultSelCount", "fPileUpMultSelCount", 8, 0., 8.);
582  fPileUpMultSelCount->GetXaxis()->SetBinLabel(1,"IsNotPileup");
583  fPileUpMultSelCount->GetXaxis()->SetBinLabel(2,"IsNotPileupMV");
584  fPileUpMultSelCount->GetXaxis()->SetBinLabel(3,"IsNotPileupInMultBins");
585  fPileUpMultSelCount->GetXaxis()->SetBinLabel(4,"InconsistentVertices");
586  fPileUpMultSelCount->GetXaxis()->SetBinLabel(5,"TrackletVsCluster");
587  fPileUpMultSelCount->GetXaxis()->SetBinLabel(6,"AsymmetricInVZERO");
588  fPileUpMultSelCount->GetXaxis()->SetBinLabel(7,"IncompleteDAQ");
589  fPileUpMultSelCount->GetXaxis()->SetBinLabel(8,"GoodVertex2016");
591 
592  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);
593  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);
594  fOutput->Add(fMultTOFLowCut);
595  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);
596  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);
597  fOutput->Add(fMultTOFHighCut);
598 
599  for(Int_t c=0; c<2; c++) {
600  for(Int_t i=0; i<5; i++) {
601  fTowerGainEq[c][i] = new TH1D();
602  fOutput->Add(fTowerGainEq[c][i]);
603  }
604  }
605 
606  for (Int_t c=0; c<2; c++) {
607  fhZNCenDis[c] = new TH3D(Form("fhZNCenDis[%d]",c), Form("fhZNCenDis[%d]",c), 100, 0., 100., 100, -2., 2. , 100., -2., 2.);
608  fOutput->Add(fhZNCenDis[c]);
609  }
610 
611  if(fBadTowerCalibList) {
612  for(Int_t c=0; c<100; c++) {
613  fBadTowerCalibHist[c] = (TH2D*)fBadTowerCalibList->FindObject(Form("TH2Resp[%d]",c));
615  }
616  }
617  if(fZDCSpectraCorrList) {
618  for(Int_t i=0; i<8; i++) {
619  SpecCorMu1[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu1[%d]",i));
620  fOutput->Add(SpecCorMu1[i]);
621  SpecCorMu2[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorMu2[%d]",i));
622  fOutput->Add(SpecCorMu2[i]);
623  SpecCorAv[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorAv[%d]",i));
624  fOutput->Add(SpecCorAv[i]);
625  SpecCorSi[i] = (TH1D*)fZDCSpectraCorrList->FindObject(Form("SpecCorSi[%d]",i));
626  fOutput->Add(SpecCorSi[i]);
627  }
628  }
629 
630  fhZNSpectra = new TH3D("fhZNSpectra","fhZNSpectra",100,0.,100.,8,0.,8.,1000,0.,1.E5);
631  fOutput->Add(fhZNSpectra);
632  fhZNSpectraCor = new TH3D("fhZNSpectraCor","fhZNSpectraCor",100,0.,100.,8,0.,8.,1000,0.,1.E5);
633  fOutput->Add(fhZNSpectraCor);
634  fhZNSpectraPow = new TH3D("fhZNSpectraPow","fhZNSpectraPow",100,0.,100.,8,0.,8.,1000,0.,TMath::Power(1.E5,fZDCGainAlpha));
635  fOutput->Add(fhZNSpectraPow);
636  fhZNBCCorr = new TH3D("fhZNBCCorr","fhZNBCCorr",100,0.,100.,500,0.,1.E5,500,0.,1.E5);
637  fOutput->Add(fhZNBCCorr);
638 
639  if(fAnalysisType == kMCAOD) {
640 
641  fSpectraMCList = new TList();
642  fSpectraMCList->SetOwner(kTRUE);
643  fSpectraMCList->SetName("Spectra");
644  fOutput->Add(fSpectraMCList);
645 
646  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.};
647  for(Int_t j=0; j<2; j++) {
648  for(Int_t c=0; c<10; c++) {
649  fPtSpecGen[j][c] = new TH1F(Form("fPtSpecGen[%d][%d]",j,c), Form("fPtSpecGen[%d][%d]",j,c), 23, xmin);
650  fSpectraMCList->Add(fPtSpecGen[j][c]);
651  fPtSpecFB32[j][c] = new TH1F(Form("fPtSpecFB32[%d][%d]",j,c), Form("fPtSpecFB32[%d][%d]",j,c), 23, xmin);
652  fSpectraMCList->Add(fPtSpecFB32[j][c]);
653  fPtSpecFB96[j][c] = new TH1F(Form("fPtSpecFB96[%d][%d]",j,c), Form("fPtSpecFB96[%d][%d]",j,c), 23, xmin);
654  fSpectraMCList->Add(fPtSpecFB96[j][c]);
655  fPtSpecFB128[j][c] = new TH1F(Form("fPtSpecFB128[%d][%d]",j,c), Form("fPtSpecFB128[%d][%d]",j,c), 23, xmin);
656  fSpectraMCList->Add(fPtSpecFB128[j][c]);
657  fPtSpecFB768[j][c] = new TH1F(Form("fPtSpecFB768[%d][%d]",j,c), Form("fPtSpecFB768[%d][%d]",j,c), 23, xmin);
658  fSpectraMCList->Add(fPtSpecFB768[j][c]);
659  }
660  }
661  }
662 
663  fAnalysisUtil = new AliAnalysisUtils();
664  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
665  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
666 
667  for(int i=0; i<5; i++){
668  char hname[20];
669  sprintf(hname,"hZNCPM%d",i);
670  fhZNCPM[i] = new TH1F(hname, hname, 200, -50., 140000);
671  fOutput->Add(fhZNCPM[i]);
672  //
673  sprintf(hname,"hZNAPM%d",i);
674  fhZNAPM[i] = new TH1F(hname, hname, 200, -50., 140000);
675  fOutput->Add(fhZNAPM[i]);
676  //
677  if(i<4){
678  //
679  char hnamenc[20];
680  sprintf(hnamenc, "hZNCPMQ%dPMC",i+1);
681  fhZNCPMQiPMC[i] = new TH1F(hnamenc, hnamenc, 100, 0., 1.);
682  fOutput->Add(fhZNCPMQiPMC[i]);
683  //
684  char hnamena[20];
685  sprintf(hnamena, "hZNAPMQ%dPMC",i+1);
686  fhZNAPMQiPMC[i] = new TH1F(hnamena, hnamena, 100, 0., 1.);
687  fOutput->Add(fhZNAPMQiPMC[i]);
688  }
689  }
690 
691  fhZNCvsZNA = new TH2F("hZNCvsZNA","hZNCvsZNA",200,-50.,140000,200,-50.,140000);
692  fOutput->Add(fhZNCvsZNA);
693  fhZDCCvsZDCCA = new TH2F("hZDCCvsZDCCA","hZDCCvsZDCCA",200,0.,180000.,200,0.,200000.);
694  fOutput->Add(fhZDCCvsZDCCA);
695  fhZNCvsZPC = new TH2F("hZNCvsZPC","hZNCvsZPC",200,-50.,50000,200,-50.,140000);
696  fOutput->Add(fhZNCvsZPC);
697  fhZNAvsZPA = new TH2F("hZNAvsZPA","hZNAvsZPA",200,-50.,50000,200,-50.,140000);
698  fOutput->Add(fhZNAvsZPA);
699  fhZNvsZP = new TH2F("hZNvsZP","hZNvsZP",200,-50.,80000,200,-50.,200000);
700  fOutput->Add(fhZNvsZP);
701  fhZNvsVZERO = new TH2F("hZNvsVZERO","hZNvsVZERO",250,0.,25000.,200,0.,200000.);
702  fOutput->Add(fhZNvsVZERO);
703  fhZDCvsVZERO = new TH2F("hZDCvsVZERO","hZDCvsVZERO",250,0.,25000.,250,0.,250000.);
704  fOutput->Add(fhZDCvsVZERO);
705  fhZDCvsTracklets = new TH2F("hZDCvsTracklets","hZDCvsTracklets",200,0.,4000.,250,0.,250000.);
707  fhZDCvsNclu1 = new TH2F("hZDCvsNclu1", "hZDCvsNclu1", 200, 0.,8000.,200,0.,250000.);
708  fOutput->Add(fhZDCvsNclu1);
709  fhDebunch = new TH2F("hDebunch","hDebunch",240,-100.,-40.,240,-30.,30.);
710  fOutput->Add(fhDebunch);
711 
712  fhAsymm = new TH1F("hAsymm" , "Asimmetry ",200,-1.,1.);
713  fOutput->Add(fhAsymm);
714  fhZNAvsAsymm = new TH2F("hZNAvsAsymm","ZNA vs. asymm.",200,-1.,1.,200,0.,80.);
715  fOutput->Add(fhZNAvsAsymm);
716  fhZNCvsAsymm = new TH2F("hZNCvsAsymm","ZNC vs. asymm.",200,-1.,1.,200,0.,80.);
717  fOutput->Add(fhZNCvsAsymm);
718 
719  fhZNCvscentrality = new TH2F("hZNCvscentrality","ZNC vs. centrality",100,0.,100.,100,0.,100.);
721  fhZNAvscentrality = new TH2F("hZNAvscentrality","ZNA vs. centrality",100,0.,100.,100,0.,100.);
723 
724  //********************************************************************
725 
726  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};
727 
728  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};
729 
730  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};
731 
732  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};
733 
734  Int_t dRun15opidfix[] = {245145, 245146, 245151, 245152, 245231, 245232, 245259, 245343, 245345, 245346, 245347, 245349, 245353, 245396, 245397, 245401, 245407, 245409, 245441, 245446, 245450, 245454, 245496, 245497, 245501, 245504, 245505, 245507, 245535, 245540, 245542, 245543, 245544, 245545, 245554};
735 
736  if(fDataSet==k2010) {fCRCnRun=92;}
737  if(fDataSet==k2011) {fCRCnRun=119;}
738  if(fDataSet==k2015) {fCRCnRun=90;}
739  if(fDataSet==k2015v6) {fCRCnRun=91;}
740  if(fDataSet==k2015pidfix) {fCRCnRun=35;}
741  if(fDataSet==kAny) {fCRCnRun=1;}
742 
743  Int_t d=0;
744  for(Int_t r=0; r<fCRCnRun; r++) {
745  if(fDataSet==k2010) fRunList[d] = dRun10h[r];
746  if(fDataSet==k2011) fRunList[d] = dRun11h[r];
747  if(fDataSet==k2015) fRunList[d] = dRun15o[r];
748  if(fDataSet==k2015v6) fRunList[d] = dRun15ov6[r];
749  if(fDataSet==k2015pidfix) fRunList[d] = dRun15opidfix[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 == kTrackQA) {
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  if(!fUseTowerEq) {
833  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.};
834  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.};
835  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};
836 
837  for(Int_t r=0;r<fCRCnRun;r++) {
838  fCRCQVecListRun[r] = new TList();
839  fCRCQVecListRun[r]->SetName(Form("Run %d",fRunList[r]));
840  fCRCQVecListRun[r]->SetOwner(kTRUE);
841  fOutput->Add(fCRCQVecListRun[r]);
842 
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  // fhZNSpectraRbR[r] = new TH3D(Form("fhZNSpectraRbR[%d]",fRunList[r]),Form("fhZNSpectraRbR[%d]",fRunList[r]),50,0.,100.,8,0.,8.,100,0.,1.E5);
853  // fCRCQVecListRun[r]->Add(fhZNSpectraRbR[r]);
854 
855  // for(Int_t i=0;i<fnCen;i++) {
856  // fPtPhiEtaRbRFB128[r][i] = new TH3F(Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),Form("fPtPhiEtaRbRFB128[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
857  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB128[r][i]);
858  // fPtPhiEtaRbRFB768[r][i] = new TH3F(Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),Form("fPtPhiEtaRbRFB768[%d][%d]",r,i),14, ptmin, 16, phimin, 16, etamin);
859  // fCRCQVecListRun[r]->Add(fPtPhiEtaRbRFB768[r][i]);
860  // }
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
894 
895  // get centrality
896  Double_t centrV0M=300, centrCL1=300, centrCL0=300, centrTRK=300;
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 
931  if(fAnalysisType == kTracklets) {
932  // fill with tracklets
933  AliAODTracklets* anInputTracklets = (AliAODTracklets*)aod->GetTracklets();
934  Int_t multSPD = anInputTracklets->GetNumberOfTracklets();
935  Double_t BField = aod->GetMagneticField();
936  //loop over tracklets
937  for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
938  Float_t thetaTr= anInputTracklets->GetTheta(itracklet);
939  Float_t phiTr= anInputTracklets->GetPhi(itracklet);
940  // calculate eta
941  Float_t etaTr = -TMath::Log(TMath::Tan(thetaTr/2.));
942  //make new AliFLowTrackSimple
943  AliFlowTrack* pTrack = new AliFlowTrack();
944  pTrack->SetPt(0.5);
945  pTrack->SetEta(etaTr);
946  // set charge: "according to Ruben, with positive magnetic field, the positive tracks rotate clockwise. Since the angle stored in AOD is the difference between the hit on the first and second layers of SPD, in positive mag field, positive delta_phi -> positive track charge."
947  Double_t DeltaPhi = anInputTracklets->GetDeltaPhi(itracklet);
948  if(BField>0. && DeltaPhi>0.) pTrack->SetCharge(1);
949  if(BField>0. && DeltaPhi<0.) pTrack->SetCharge(-1);
950  if(BField<0. && DeltaPhi>0.) pTrack->SetCharge(-1);
951  if(BField<0. && DeltaPhi<0.) pTrack->SetCharge(1);
952  // correction of phi
954  phiTr += 39./34.*DeltaPhi;
955  if (phiTr < 0.) phiTr += 2.*TMath::Pi();
956  if (phiTr > 2.*TMath::Pi()) phiTr -= 2.*TMath::Pi();
957  }
958  pTrack->SetPhi(phiTr);
959  //marking the particles as POI type 2:
961  pTrack->SetPOItype(2,kTRUE);
963  //Add the track to the flowevent
964  fFlowEvent->AddTrack(pTrack);
965  }
966  }
967 
969 
974  fFlowEvent->SetCentralityCL1(centrCL1);
975  fFlowEvent->SetCentralityTRK(centrTRK);
976  // fFlowEvent->SetNITSCL1(((AliVAODHeader*)aod->GetHeader())->GetNumberOfITSClusters(1));
977  Double_t SumV0=0.;
978  for(Int_t i=0; i<64; i++) {
979  if(std::isfinite(aod->GetVZEROEqMultiplicity(i))) SumV0 += aod->GetVZEROEqMultiplicity(i);
980  }
981  fFlowEvent->SetNITSCL1(SumV0);
982 
983  // set absolute orbit number
984  UInt_t period = aod->GetPeriodNumber();
985  UInt_t orbit24 = aod->GetOrbitNumber(); // wrapped down to 24 bits
986  if (period > 255) { // 8 bits
987  period = 255;
988  orbit24 = (1<<24) - 1;
989  }
990  if (orbit24 >= (1<<24)) { // 24 bits
991  period = 255;
992  orbit24 = (1<<24) - 1;
993  }
994  UInt_t orbit = period * (1<<24) + orbit24;
995  fFlowEvent->SetAbsOrbit(orbit);
996 
997  Double_t vtxpos[3]={0.,0.,0.};
998  vtxpos[0] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetX();
999  vtxpos[1] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetY();
1000  vtxpos[2] = ((AliAODVertex*)aod->GetPrimaryVertex())->GetZ();
1001  fFlowEvent->SetVertexPosition(vtxpos);
1002 
1003  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1004 
1005  // run-by-run QA
1006  // for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1007  // AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1008  // if(!track) continue;
1009  // // general kinematic & quality cuts
1010  // if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70) continue;
1011  // if (track->TestFilterBit(128)) fPtPhiEtaRbRFB128[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1012  // if (track->TestFilterBit(768)) fPtPhiEtaRbRFB768[RunBin][CenBin]->Fill(track->Pt(),track->Phi(),track->Eta());
1013  // }
1014  fCenDis->Fill(centrV0M);
1015 
1016  }
1017 
1018  if (fAnalysisType == kTrackQA) {
1019 
1020  // empty flow event
1021  fFlowEvent->ClearFast();
1022 
1023  // get centrality
1024  Double_t centrV0M=300;
1026  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
1027  } else {
1028  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
1029  if(!fMultSelection) {
1030  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
1031  AliWarning("AliMultSelection object not found!");
1032  } else {
1033  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
1034  }
1035  }
1036  if(centrV0M<10. || centrV0M>40.) return;
1037 
1038  //check event cuts
1039  if (InputEvent()) {
1040  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
1041  if(fRejectPileUp) {
1042  Bool_t IsPileUp = SelectPileup(aod);
1043  if(IsPileUp) return;
1044  }
1045  }
1046 
1047  Double_t BField = aod->GetMagneticField();
1048 
1049  for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1050 
1051  AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1052  if(!track) continue;
1053 
1054  Double_t dPt = track->Pt();
1055  Double_t dEta = track->Eta();
1056  Double_t dPhi = track->Phi();
1057  Double_t dCharge = track->Charge();
1058  Int_t cw = 0;
1059  if(BField>0. && dCharge>0.) cw=0;
1060  if(BField>0. && dCharge<0.) cw=1;
1061  if(BField<0. && dCharge>0.) cw=2;
1062  if(BField<0. && dCharge<0.) cw=3;
1063 
1064  // general kinematic cuts
1065  if (dPt < .2 || dPt > 50. || TMath::Abs(dEta) > 0.8) continue;
1066 
1067  // cut on DCA
1068  Double_t DCAxy = track->DCA();
1069  Double_t DCAz = track->ZAtDCA();
1070  if (std::abs((Int_t)DCAxy)==999 || std::abs((Int_t)DCAz)==999) {
1071  // re-evaluate the dca as it seems to not be natively present
1072  // allowed only for tracks inside the beam pipe
1073  Double_t pos[3] = {-99., -99., -99.};
1074  track->GetPosition(pos);
1075  if(pos[0]*pos[0]+pos[1]*pos[1] <= 3.*3.) {
1076  AliAODTrack copy(*track); // stack copy
1077  Double_t b[2] = {-99., -99.};
1078  Double_t bCov[3] = {-99., -99., -99.};
1079  if(copy.PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov)) {
1080  DCAxy = b[0];
1081  DCAz = b[1];
1082  }
1083  }
1084  }
1085  if(fabs(DCAxy)>2.4 || fabs(DCAz)>3.2) continue;
1086 
1087  // various cuts on TPC clusters
1088  if (track->GetTPCNcls() < 70) continue;
1089  Double_t chi2_per_tpc = track->Chi2perNDF();
1090  if (chi2_per_tpc < 0.1 || chi2_per_tpc > 4.) continue;
1091  Double_t fraction_shared_tpccls = 1.*track->GetTPCnclsS()/track->GetTPCncls();
1092  if (fraction_shared_tpccls > 0.4) continue;
1093 
1094  Int_t FBarray[4] = {32,96,128,768};
1095 
1096  // test filter bits
1097  for(Int_t fb=0; fb<fKNFBs; fb++){
1098  if (track->TestFilterBit(FBarray[fb])) {
1099  fTrackQADCAxy[fb][cw]->Fill(dEta,dPt,DCAxy);
1100  fTrackQADCAz[fb][cw]->Fill(dEta,dPt,DCAz);
1101  fTrackQApT[fb][cw]->Fill(dEta,dPt);
1102  fEbEQRe[fb][cw]->Fill(dEta,dPt,TMath::Cos(dPhi));
1103  fEbEQIm[fb][cw]->Fill(dEta,dPt,TMath::Sin(dPhi));
1104  fEbEQMu[fb][cw]->Fill(dEta,dPt);
1105  }
1106  }
1107  }
1108 
1109  // compute cos(#Delta#phi)
1110  for(Int_t bx=1; bx<=fEbEQRe[0][0]->GetXaxis()->GetNbins(); bx++) {
1111  for(Int_t by=1; by<=fEbEQRe[0][0]->GetYaxis()->GetNbins(); by++) {
1112 
1113  Double_t dEta = fEbEQRe[0][0]->GetXaxis()->GetBinCenter(bx);
1114  Double_t dPt = fEbEQRe[0][0]->GetYaxis()->GetBinCenter(by);
1115 
1116  for(Int_t fb=0; fb<fKNFBs; fb++){
1117  for(Int_t c=0; c<4; c++){
1118 
1119  Double_t QRe = fEbEQRe[fb][c]->GetBinContent(bx,by);
1120  Double_t QIm = fEbEQIm[fb][c]->GetBinContent(bx,by);
1121  Double_t M = 1.*fEbEQMu[fb][c]->GetBinContent(bx,by);
1122 
1123  if(M>1.) {
1124  Double_t c1 = (QRe*QRe+QIm*QIm-M)/(M*(M-1.));
1125  fTrackQADphi[fb][c]->Fill(dEta,dPt,c1);
1126  }
1127  }
1128  }
1129 
1130  }
1131  }
1132 
1133  // reset EbE histograms
1134  for(Int_t fb=0; fb<fKNFBs; fb++){
1135  for(Int_t c=0; c<4; c++){
1136  fEbEQRe[fb][c]->Reset();
1137  fEbEQIm[fb][c]->Reset();
1138  fEbEQMu[fb][c]->Reset();
1139  }
1140  }
1141 
1142  }
1143 
1144  if (fAnalysisType == kMCAOD) {
1145 
1146  //check event cuts
1147  if (InputEvent()) {
1148  if(!fCutsEvent->IsSelected(InputEvent(),MCEvent())) return;
1149  if(fRejectPileUp && fAnalysisUtil->IsPileUpEvent(InputEvent())) return;
1150  }
1151 
1152  fFlowEvent->ClearFast();
1153 
1154  if(!McEvent) {
1155  AliError("ERROR: Could not retrieve MCEvent");
1156  return;
1157  }
1158  fStack = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
1159  if(!fStack){
1160  AliError("ERROR: Could not retrieve MCStack");
1161  return;
1162  }
1163 
1164  // get centrality (from AliMultSelection or AliCentrality)
1165  Double_t centr = 300;
1167  fMultSelection = (AliMultSelection*)aod->FindListObject("MultSelection");
1168  if(!fMultSelection) {
1169  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
1170  AliWarning("AliMultSelection object not found!");
1171  } else {
1172  centr = fMultSelection->GetMultiplicityPercentile("V0M");
1173  }
1174  } else {
1175  centr = (((AliVAODHeader*)aod->GetHeader())->GetCentralityP())->GetCentralityPercentile("V0M");
1176  }
1177  // centrality bin
1178  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1179  Int_t CenBin = -1;
1180  CenBin = GetCenBin(centr);
1181  if(CenBin==-1) return;
1182  fCenDis->Fill(centr);
1183 
1184  // reconstructed
1185  for(Int_t jTracks = 0; jTracks<aod->GetNumberOfTracks(); jTracks++){
1186 
1187  AliAODTrack* track = (AliAODTrack*)aod->GetTrack(jTracks);
1188  if(!track) continue;
1189 
1190  // to select primaries
1191  Int_t lp = TMath::Abs(track->GetLabel());
1192 
1193  // general kinematic cuts
1194  if (track->Pt() < .2 || track->Pt() > 20. || TMath::Abs(track->Eta()) > 0.8) continue;
1195 
1196  // cut on DCA
1197  Double_t DCAxy = track->DCA();
1198  Double_t DCAz = track->ZAtDCA();
1199  if(fabs(DCAxy)>2.4 || fabs(DCAz)>3.2) continue;
1200 
1201  // various cuts on TPC clusters
1202  if (track->GetTPCNcls() < 70) continue;
1203  Double_t chi2_per_tpc = track->Chi2perNDF();
1204  if (chi2_per_tpc < 0.1 || chi2_per_tpc > 4.) continue;
1205  Double_t fraction_shared_tpccls = 1.*track->GetTPCnclsS()/track->GetTPCncls();
1206  if (fraction_shared_tpccls > 0.4) continue;
1207 
1208  // test filter bits
1209  if (((AliAODMCParticle*)fStack->At(lp))->IsPhysicalPrimary()) {
1210  if (track->TestFilterBit(32)) fPtSpecFB32[0][CenBin]->Fill(track->Pt());
1211  if (track->TestFilterBit(96)) fPtSpecFB96[0][CenBin]->Fill(track->Pt());
1212  if (track->TestFilterBit(128)) fPtSpecFB128[0][CenBin]->Fill(track->Pt());
1213  if (track->TestFilterBit(768)) fPtSpecFB768[0][CenBin]->Fill(track->Pt());
1214  } else {
1215  if (track->TestFilterBit(32)) fPtSpecFB32[1][CenBin]->Fill(track->Pt());
1216  if (track->TestFilterBit(96)) fPtSpecFB96[1][CenBin]->Fill(track->Pt());
1217  if (track->TestFilterBit(128)) fPtSpecFB128[1][CenBin]->Fill(track->Pt());
1218  if (track->TestFilterBit(768)) fPtSpecFB768[1][CenBin]->Fill(track->Pt());
1219  }
1220 
1221  }
1222 
1223  // generated (physical primaries)
1224 
1225  for(Int_t jTracks = 0; jTracks<fStack->GetEntriesFast(); jTracks++) {
1226  AliAODMCParticle *MCpart = (AliAODMCParticle*)fStack->At(jTracks);
1227  if (!MCpart) {
1228  printf("ERROR: Could not receive MC track %d\n", jTracks);
1229  continue;
1230  }
1231 
1232  // kinematic cuts
1233  if ( MCpart->Pt() < .2 || MCpart->Pt() > 20. || TMath::Abs(MCpart->Eta()) > .8 ) continue;
1234  // select charged primaries
1235  if ( MCpart->Charge() == 0. || !MCpart->IsPhysicalPrimary()) continue;
1236 
1237  fPtSpecGen[0][CenBin]->Fill(MCpart->Pt());
1238  }
1239 
1240  // fGenHeader = McEvent->GenEventHeader();
1241  // if(fGenHeader) fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1242  // printf("#reconstructed : %d (rejected from cuts %d), #MC primaries : %d (rejected from cuts %d) \n",AODPOIs,AODbads,MCPrims,MCSecos);
1243  fFlowEvent->SetReferenceMultiplicity(aod->GetNumberOfTracks());
1244  fFlowEvent->SetCentrality(centr);
1245  // if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1246  fFlowEvent->SetRun(RunNum);
1247  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1248  }
1249 
1250  if(fAnalysisType == kMCESD) {
1251 
1252  fFlowEvent->ClearFast();
1253 
1254  if(!esd) {
1255  AliError("ERROR: Could not retrieve ESDEvent");
1256  return;
1257  }
1258  if(!McEvent) {
1259  AliError("ERROR: Could not retrieve MCEvent");
1260  return;
1261  }
1262  AliStack* fStack = fMCEvent->Stack();
1263  if(!fStack) {
1264  AliError("ERROR: Could not retrieve MCStack");
1265  return;
1266  }
1267 
1268  AliESDVertex *vertex = (AliESDVertex*) esd->GetPrimaryVertex();
1269  if (!vertex) return;
1270  if (TMath::Abs(vertex->GetZ()) > 10. ) return;
1271  if (vertex->GetNContributors() < 1 ) return;
1272  AliCentrality *centrality = esd->GetCentrality();
1273  if (!centrality) return;
1274  Double_t centr = centrality->GetCentralityPercentile("V0M");
1275  if (centr<fCentrLowLim || centr>=fCentrUpLim ) return;
1276  Int_t CenBin = -1;
1277  if (centr>0. && centr<5.) CenBin=0;
1278  if (centr>5. && centr<10.) CenBin=1;
1279  if (centr>10. && centr<20.) CenBin=2;
1280  if (centr>20. && centr<30.) CenBin=3;
1281  if (centr>30. && centr<40.) CenBin=4;
1282  if (centr>40. && centr<50.) CenBin=5;
1283  if (centr>50. && centr<60.) CenBin=6;
1284  if (centr>60. && centr<70.) CenBin=7;
1285  if (centr>70. && centr<80.) CenBin=8;
1286  if (centr>80. && centr<90.) CenBin=9;
1287  if(CenBin==-1) return;
1288 
1289  //Generated
1290  Int_t MCPrims = 0;
1291  for ( Int_t i=0 ; i<fStack->GetNtrack() ; i++ ) {
1292 
1293  //Primaries Selection
1294  TParticle *particle = (TParticle*)fStack->Particle(i);
1295  if (!particle) continue;
1296  if (!fStack->IsPhysicalPrimary(i)) continue;
1297  if ( particle->GetPDG()->Charge() == 0.) continue;
1298 
1299  //Kinematic Cuts
1300  if ( particle->Pt()<0.2 || particle->Pt()>10. ) continue;
1301  if ( TMath::Abs(particle->Eta())>0.8 ) continue;
1302 
1303  fFlowTrack->SetPhi(particle->Phi());
1304  fFlowTrack->SetEta(particle->Eta());
1305  fFlowTrack->SetPt(particle->Pt());
1307  fFlowTrack->SetForRPSelection(kTRUE);
1309  fFlowTrack->SetForPOISelection(kFALSE);
1311  MCPrims++;
1312 
1313  fPtSpecGen[0][CenBin]->Fill(particle->Pt());
1314 
1315  }
1316 
1317  //Reconstructed
1318  Int_t ESDPrims = 0;
1319  for (Int_t i=0 ; i<esd->GetNumberOfTracks() ; i++) {
1320 
1321  //Get reconstructed track
1322  AliVTrack *vtrack = static_cast<AliVTrack*>(esd->GetTrack(i));
1323  AliESDtrack *track = dynamic_cast<AliESDtrack*>(vtrack);
1324  if (!track) continue;
1325 
1326  //Primaries selection
1327  Int_t lp = TMath::Abs(track->GetLabel());
1328  if (!fStack->IsPhysicalPrimary(lp)) continue;
1329  TParticle *particle = (TParticle*)fStack->Particle(lp);
1330  if (!particle) continue;
1331  if (particle->GetPDG()->Charge() == 0.) continue;
1332 
1333  // if(!fCutsPOI->PassesESDcuts(track)) continue;
1334 
1335  Bool_t pass = kTRUE;
1336 
1337  if(fCutTPC) {
1338  // printf("******* cutting TPC ******** \n");
1339  UShort_t ntpccls = track->GetTPCNcls();
1340  Double_t tpcchi2 = track->GetTPCchi2();
1341  if (tpcchi2<0.2 || tpcchi2 >=4.) {
1342  // printf("TPCchi2 : %e %e ",tpcchi2,track->GetTPCchi2Iter1());
1343  pass=kFALSE;
1344  }
1345  if (ntpccls < 70) {
1346  // printf("#TPCcluster : %u %u %u %u ",ntpccls,track->GetTPCNclsF(),track->GetTPCNclsFIter1(),track->GetTPCNclsIter1());
1347  pass=kFALSE;
1348  }
1349  }
1350 
1351  Float_t dcaxy=0.0;
1352  Float_t dcaz=0.0;
1353  track->GetImpactParameters(dcaxy,dcaz);
1354  if (dcaxy > 0.3 || dcaz > 0.3) {
1355  // printf("DCA : %e %e ",dcaxy,dcaz);
1356  pass=kFALSE;
1357  }
1358  if(!pass) continue;
1359 
1360  //Kinematic Cuts
1361  if ( track->Pt()<0.2 || track->Pt()>10. ) continue;
1362  if ( TMath::Abs(track->Eta())>0.8 ) continue;
1363 
1364  fFlowTrack->SetPhi(track->Phi());
1365  fFlowTrack->SetEta(track->Eta());
1366  fFlowTrack->SetPt(track->Pt());
1368  fFlowTrack->SetForRPSelection(kFALSE);
1372  ESDPrims++;
1373 
1374  }
1375 
1376  // printf("#reconstructed : %d , #MC primaries : %d \n",ESDPrims,MCPrims);
1377  fFlowEvent->SetReferenceMultiplicity(esd->GetNumberOfTracks());
1378  fFlowEvent->SetCentrality(centr);
1379  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1380  fFlowEvent->SetRun(esd->GetRunNumber());
1381  // printf("Run : %d, RefMult : %d, Cent : %f \n",fFlowEvent->GetRun(),fFlowEvent->GetReferenceMultiplicity(),fFlowEvent->GetCentrality());
1382 
1383  } // end of if(fAnalysisType == kMCESD)
1384 
1385  if(fAnalysisType == kMCkine) {
1386 
1387  fFlowEvent->ClearFast();
1388 
1389  AliInputEventHandler* McHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1390  if(!McHandler) {
1391  AliError("ERROR: Could not retrieve MCtruthEventHandler");
1392  return;
1393  }
1394  McEvent = McHandler->MCEvent();
1395  if(!McEvent) {
1396  AliError("ERROR: Could not retrieve MC event");
1397  return;
1398  }
1399 
1400  Int_t nTracks = McEvent->GetNumberOfTracks();
1401  // Int_t nPrimTr = McEvent->GetNumberOfPrimaries();
1402 
1403  //loop over tracks
1404  for (Int_t itrkN=0; itrkN<nTracks; itrkN++) {
1405  //get input particle
1406  AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(McEvent->GetTrack(itrkN));
1407  if (!pParticle) continue;
1408 
1409  //check if track passes the cuts
1410  if (McEvent->IsPhysicalPrimary(itrkN) && pParticle->Charge()!=0) {
1411  fFlowTrack->Set(pParticle);
1413  fFlowTrack->SetForRPSelection(kTRUE);
1418  }
1419  }// for all tracks
1420 
1421  // if monte carlo event get reaction plane from monte carlo (depends on generator)
1422  if (McEvent && McEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(McEvent);
1423  // set reference multiplicity
1424  fFlowEvent->SetReferenceMultiplicity(McEvent->GetNumberOfTracks());
1425  // tag subevents
1427  // set centrality from impact parameter
1428  Double_t ImpPar=0., CenPer=0.;
1429  fGenHeader = McEvent->GenEventHeader();
1430  if(fGenHeader){
1431  fPythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fGenHeader);
1432  if(fPythiaGenHeader) ImpPar = fPythiaGenHeader->GetImpactParameter();
1433  fHijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(fGenHeader);
1434  if(fHijingGenHeader) ImpPar = fHijingGenHeader->ImpactParameter();
1435  if(ImpPar) CenPer = 0.3859796743103508*pow(ImpPar,2.);
1436  if(CenPer>0. && CenPer<100.) fFlowEvent->SetCentrality(CenPer);
1437  else return;
1438  fFlowEvent->SetRun(1);
1439  }
1440 
1441  } // end of if(fAnalysisType == kMCkine)
1442 
1443  if (!fFlowEvent) return; //shuts up coverity
1444 
1445  //check final event cuts
1446  Int_t mult = fFlowEvent->NumberOfTracks();
1447  // AliInfo(Form("FlowEvent has %i tracks",mult));
1448  if (mult<fMinMult || mult>fMaxMult) {
1449  AliWarning("FlowEvent cut on multiplicity"); return;
1450  }
1451 
1452  //define dead zone
1454 
1457  if (fAfterburnerOn)
1458  {
1459  //if reaction plane not set from elsewhere randomize it before adding flow
1461  fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
1462 
1463  if(fDifferentialV2)
1465  else
1466  fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
1467  fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
1468  }
1470 
1471  //tag subEvents
1473 
1474  //do we want to serve shullfed tracks to everybody?
1476 
1477  // associate the mother particles to their daughters in the flow event (if any)
1479 
1480  //fListHistos->Print();
1481  //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
1482 
1483  //********************************************************************************************************************************
1484 
1486 
1487  // PHYSICS SELECTION
1488  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1489  AliInputEventHandler *hdr = (AliInputEventHandler*)am->GetInputEventHandler();
1490 
1491  if(hdr->IsEventSelected() && AliVEvent::kAny) {
1492 
1493  Double_t centrperc = fFlowEvent->GetCentrality();
1494  Int_t cenb = (Int_t)centrperc;
1495 
1496  AliAODTracklets *trackl = aod->GetTracklets();
1497  Int_t nTracklets = trackl->GetNumberOfTracklets();
1498 
1499  // get VZERO data
1500  AliAODVZERO *vzeroAOD = aod->GetVZEROData();
1501  Double_t multV0A = vzeroAOD->GetMTotV0A();
1502  Double_t multV0C = vzeroAOD->GetMTotV0C();
1503 
1504  // set VZERO Q-vectors
1505  if(fDataSet==k2015 || fDataSet==k2015v6) {
1506  Int_t CachednRing = 1;
1507  Double_t QxTot[fkVZEROnHar] = {0.}, QyTot[fkVZEROnHar] = {0.};
1508  Double_t denom = 0.;
1509  Double_t V0TotQC[fkVZEROnHar][2] = {{0.}}, V0TotQA[fkVZEROnHar][2] = {{0.}};
1510  Double_t MultC[fkVZEROnHar] = {0.}, MultA[fkVZEROnHar] = {0.};
1511 
1512  for(Int_t i=0; i<64; i++) {
1513 
1514  // correct multiplicity per channel
1515  Double_t mult = vzeroAOD->GetMultiplicity(i);
1516  if(fVZEROGainEqHist) {
1517  Double_t EqFactor = fVZEROGainEqHist->GetBinContent(RunBin+1,i+1);
1518  if(EqFactor>0.) mult *= EqFactor;
1519  }
1520  fVZEROMult->Fill(RunBin+0.5,i+0.5,mult);
1521 
1522  // build Q-vector per ring
1523  Int_t nRing = (Int_t)i/8 + 1;
1524  Double_t ChPhi = TMath::PiOver4()*(0.5+i%8);
1525 
1526  if(i == 63) {
1527  for (Int_t k=0; k<fkVZEROnHar; k++) {
1528  QxTot[k] += mult*TMath::Cos((k+1.)*ChPhi);
1529  QyTot[k] += mult*TMath::Sin((k+1.)*ChPhi);
1530  }
1531  denom += mult;
1532  nRing++;
1533  }
1534 
1535  if(nRing!=CachednRing) {
1536  for (Int_t k=0; k<fkVZEROnHar; k++) {
1537  Double_t QxRec = QxTot[k]/denom;
1538  Double_t QyRec = QyTot[k]/denom;
1539  // store values for re-centering
1540  // fVZEROQVectorRecQx[k]->Fill(RunBin+0.5,centrperc,CachednRing-0.5,QxRec);
1541  // fVZEROQVectorRecQy[k]->Fill(RunBin+0.5,centrperc,CachednRing-0.5,QyRec);
1542  // do re-centering
1543  if(fVZEROQVectorRecQxStored[k]) {
1544  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));
1545  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));
1546  }
1547  // sum of Q-vectors over all rings (total V0 Q-vector)
1548  if (CachednRing >= fMinRingVZC && CachednRing <= fMaxRingVZC) {
1549  V0TotQC[k][0] += QxRec*denom;
1550  V0TotQC[k][1] += QyRec*denom;
1551  MultC[k] += denom;
1552  }
1553  if (CachednRing >= fMinRingVZA && CachednRing <= fMaxRingVZA) {
1554  V0TotQA[k][0] += QxRec*denom;
1555  V0TotQA[k][1] += QyRec*denom;
1556  MultA[k] += denom;
1557  }
1558  QxTot[k] = 0.;
1559  QyTot[k] = 0.;
1560  }
1561  denom = 0.;
1562  CachednRing = nRing;
1563  }
1564  for (Int_t k=0; k<fkVZEROnHar; k++) {
1565  QxTot[k] += mult*TMath::Cos((k+1.)*ChPhi);
1566  QyTot[k] += mult*TMath::Sin((k+1.)*ChPhi);
1567  }
1568  denom += mult;
1569  }
1570 
1571  for (Int_t k=0; k<fkVZEROnHar; k++) {
1572  if(MultC[k]>0. && MultA[k]>0.) {
1573  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];
1574  if(!std::isnan(QCx) && !std::isnan(QCy) && !std::isnan(QAx) && !std::isnan(QAy)) {
1575  fFlowEvent->SetV02Qsub(QCx,QCy,MultC[k],QAx,QAy,MultA[k],k+1);
1576  fVZEROQVectorRecFinal[k][0]->Fill(RunBin+0.5,centrperc,QCx);
1577  fVZEROQVectorRecFinal[k][1]->Fill(RunBin+0.5,centrperc,QCy);
1578  fVZEROQVectorRecFinal[k][2]->Fill(RunBin+0.5,centrperc,QAx);
1579  fVZEROQVectorRecFinal[k][3]->Fill(RunBin+0.5,centrperc,QAy);
1580  fVZEROQVectorRecFinal[k][4]->Fill(RunBin+0.5,centrperc,QCx*QAx);
1581  fVZEROQVectorRecFinal[k][5]->Fill(RunBin+0.5,centrperc,QCy*QAy);
1582  fVZEROQVectorRecFinal[k][6]->Fill(RunBin+0.5,centrperc,QCx*QAy);
1583  fVZEROQVectorRecFinal[k][7]->Fill(RunBin+0.5,centrperc,QCy*QAx);
1584  } else {
1585  fFlowEvent->SetV02Qsub(0.,0.,0.,0.,0.,0.,k+1);
1586  }
1587  } else {
1588  fFlowEvent->SetV02Qsub(0.,0.,0.,0.,0.,0.,k+1);
1589  }
1590  }
1591  }
1592 
1593 // AliAODForwardMult* aodForward = static_cast<AliAODForwardMult*>(aodEvent->FindListObject("Forward"));
1594 // const TH2D& d2Ndetadphi = aodForward->GetHistogram();
1595 // Int_t nEta = d2Ndetadphi.GetXaxis()->GetNbins();
1596 // Int_t nPhi = d2Ndetadphi.GetYaxis()->GetNbins();
1597 // Double_t ret = 0.;
1598 // // Loop over eta
1599 // for (Int_t iEta = 1; iEta <= nEta; iEta++) {
1600 // Int_t valid = d2Ndetadphi.GetBinContent(iEta, 0);
1601 // if (!valid) continue; // No data expected for this eta
1602 // // Loop over phi
1603 // for (Int_t iPhi = 1; i <= nPhi; i++) {
1604 // ret = d2Ndetadphi.GetBinContent(iEta, iPhi);
1605 // printf("eta %e phi %e : %e \n",d2Ndetadphi.GetXaxis()->GetBinCenter(iEta),d2Ndetadphi.GetYaxis()->GetBinCenter(iPhi),ret);
1606 // }
1607 // }
1608 
1609  AliAODZDC *aodZDC = aod->GetZDCData();
1610 
1611  const Double_t * towZNCraw = aodZDC->GetZNCTowerEnergy();
1612  const Double_t * towZNAraw = aodZDC->GetZNATowerEnergy();
1613 
1614  // Get centroid from ZDCs *******************************************************
1615 
1616  Double_t Enucl = (RunNum < 209122 ? 1380. : 2511.);
1617  Double_t xyZNC[2]={0.,0.}, xyZNA[2]={0.,0.};
1618  Double_t towZNC[5]={0.}, towZNA[5]={0.};
1619 
1620  Double_t ZNCcalib=1., ZNAcalib=1.;
1621  if(fUseTowerEq) {
1622  if(RunNum!=fCachedRunNum) {
1623  for(Int_t i=0; i<5; i++) {
1624  fTowerGainEq[0][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNCTower[%d][%d]",RunNum,i)));
1625  fTowerGainEq[1][i] = (TH1D*)(fTowerEqList->FindObject(Form("fZNATower[%d][%d]",RunNum,i)));
1626  }
1627  }
1628  for(Int_t i=0; i<5; i++) {
1629  if(fTowerGainEq[0][i]) towZNC[i] = towZNCraw[i]*fTowerGainEq[0][i]->GetBinContent(fTowerGainEq[0][i]->FindBin(centrperc));
1630  if(fTowerGainEq[1][i]) towZNA[i] = towZNAraw[i]*fTowerGainEq[1][i]->GetBinContent(fTowerGainEq[1][i]->FindBin(centrperc));
1631  if(fResetNegativeZDC) {
1632  if(towZNC[i]<0.) towZNC[i] = 0.;
1633  if(towZNA[i]<0.) towZNA[i] = 0.;
1634  }
1635  }
1636  } else {
1637  for(Int_t i=0; i<5; i++) {
1638  towZNC[i] = towZNCraw[i];
1639  towZNA[i] = towZNAraw[i];
1640  if(fResetNegativeZDC) {
1641  if(towZNC[i]<0.) towZNC[i] = 0.;
1642  if(towZNA[i]<0.) towZNA[i] = 0.;
1643  }
1644  fZNCTower[RunBin][i]->Fill(centrperc,towZNC[i]);
1645  fZNATower[RunBin][i]->Fill(centrperc,towZNA[i]);
1646  }
1647  }
1648 
1649  if(RunNum>=245829) towZNA[2] = 0.;
1650  Double_t zncEnergy=0., znaEnergy=0.;
1651  for(Int_t i=0; i<5; i++){
1652  zncEnergy += towZNC[i];
1653  znaEnergy += towZNA[i];
1654  }
1655  if(RunNum>=245829) znaEnergy *= 8./7.;
1656  fFlowEvent->SetZNCQ0(towZNC[0]);
1657  fFlowEvent->SetZNAQ0(towZNA[0]);
1658 
1659  Double_t energyZNC = ((AliVAODHeader*)aod->GetHeader())->GetZDCN1Energy();
1660  Double_t energyZNA = ((AliVAODHeader*)aod->GetHeader())->GetZDCN2Energy();
1661  fFlowEvent->SetZNCEnergy(energyZNC);
1662  fFlowEvent->SetZNAEnergy(energyZNA);
1663 
1664  const Double_t x[4] = {-1.75, 1.75, -1.75, 1.75};
1665  const Double_t y[4] = {-1.75, -1.75, 1.75, 1.75};
1666  Double_t numXZNC=0., numYZNC=0., denZNC=0., cZNC, wZNC, EZNC, SumEZNC=0.;
1667  Double_t numXZNA=0., numYZNA=0., denZNA=0., cZNA, wZNA, EZNA, SumEZNA=0., BadChOr;
1668  Bool_t fAllChONZNC=kTRUE, fAllChONZNA=kTRUE;
1669 
1670  if (fUseMCCen) {
1671  for(Int_t i=0; i<4; i++){
1672 
1673  // get energy
1674  EZNC = towZNC[i+1];
1675  fhZNSpectra->Fill(centrperc,i+0.5,EZNC);
1676 // fhZNSpectraRbR[RunBin]->Fill(centrperc,i+0.5,EZNC);
1677  if(fUseZDCSpectraCorr && EZNC>0.) {
1678  Double_t mu1 = SpecCorMu1[i]->Interpolate(centrperc);
1679  Double_t mu2 = SpecCorMu2[i]->Interpolate(centrperc);
1680  Double_t av = SpecCorAv[i]->Interpolate(centrperc);
1681  Double_t cor1 = SpecCorSi[i]->Interpolate(centrperc);
1682  EZNC = exp( (log(EZNC) - mu1 + mu2*cor1)/cor1 ) + av;
1683  fhZNSpectraCor->Fill(centrperc,i+0.5,EZNC);
1684  }
1685  if(fUseZDCSpectraCorr && EZNC<=0.) fAllChONZNC=kFALSE;
1686 
1687  SumEZNC += EZNC;
1688 
1689  // build centroid
1690  wZNC = TMath::Power(EZNC, fZDCGainAlpha);
1691  numXZNC += x[i]*wZNC;
1692  numYZNC += y[i]*wZNC;
1693  denZNC += wZNC;
1694  fhZNSpectraPow->Fill(centrperc,i+0.5,wZNC);
1695 
1696  // get energy
1697  if(fDataSet==k2015 || fDataSet==k2015v6) {
1698  if(i==1) {
1699  EZNA = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1700  if(fUseBadTowerCalib && fBadTowerCalibHist[cenb]) {
1701  EZNA = GetBadTowerResp(EZNA, fBadTowerCalibHist[cenb]);
1702  }
1703  } else {
1704  EZNA = towZNA[i+1];
1705  }
1706  } else {
1707  EZNA = towZNA[i+1];
1708  }
1709  fhZNSpectra->Fill(centrperc,i+4.5,EZNA);
1710 // fhZNSpectraRbR[RunBin]->Fill(centrperc,i+4.5,EZNA);
1711  if(fUseZDCSpectraCorr && EZNA>0.) {
1712  Double_t mu1 = SpecCorMu1[i+4]->Interpolate(centrperc);
1713  Double_t mu2 = SpecCorMu2[i+4]->Interpolate(centrperc);
1714  Double_t av = SpecCorAv[i+4]->Interpolate(centrperc);
1715  Double_t cor1 = SpecCorSi[i+4]->Interpolate(centrperc);
1716  EZNA = exp( (log(EZNA) - mu1 + mu2*cor1)/cor1 ) + av;
1717  fhZNSpectraCor->Fill(centrperc,i+4.5,EZNA);
1718  }
1719  if(fUseZDCSpectraCorr && EZNA<=0.) fAllChONZNA=kFALSE;
1720  SumEZNA += EZNA;
1721 
1722  // build centroid
1723  wZNA = TMath::Power(EZNA, fZDCGainAlpha);
1724  numXZNA += x[i]*wZNA;
1725  numYZNA += y[i]*wZNA;
1726  denZNA += wZNA;
1727  fhZNSpectraPow->Fill(centrperc,i+4.5,wZNA);
1728  }
1729  // store distribution for unfolding
1730  if(RunNum<245829) {
1731  Double_t recoE = towZNA[0]-towZNA[1]-towZNA[3]-towZNA[4];
1732  Double_t trueE = towZNA[2];
1733  fhZNBCCorr->Fill(centrperc,trueE,recoE);
1734  }
1735  if(denZNC>0.){
1736  Double_t nSpecnC = SumEZNC/Enucl;
1737  cZNC = 1.89358-0.71262/(nSpecnC+0.71789);
1738  xyZNC[0] = cZNC*numXZNC/denZNC;
1739  xyZNC[1] = cZNC*numYZNC/denZNC;
1740  denZNC *= cZNC;
1741  }
1742  else{
1743  xyZNC[0] = xyZNC[1] = 0.;
1744  }
1745  if(denZNA>0.){
1746  Double_t nSpecnA = SumEZNA/Enucl;
1747  cZNA = 1.89358-0.71262/(nSpecnA+0.71789);
1748  xyZNA[0] = cZNA*numXZNA/denZNA;
1749  xyZNA[1] = cZNA*numYZNA/denZNA;
1750  denZNA *= cZNA;
1751  }
1752  else{
1753  xyZNA[0] = xyZNA[1] = 0.;
1754  }
1755  } else {
1756  for(Int_t i=0; i<4; i++) {
1757  if(towZNC[i+1]>0.) {
1758  wZNC = TMath::Power(towZNC[i+1], fZDCGainAlpha);
1759  numXZNC += x[i]*wZNC;
1760  numYZNC += y[i]*wZNC;
1761  denZNC += wZNC;
1762  }
1763  if(towZNA[i+1]>0.) {
1764  wZNA = TMath::Power(towZNA[i+1], fZDCGainAlpha);
1765  numXZNA += x[i]*wZNA;
1766  numYZNA += y[i]*wZNA;
1767  denZNA += wZNA;
1768  }
1769  }
1770  if(denZNC!=0) {
1771  xyZNC[0] = numXZNC/denZNC;
1772  xyZNC[1] = numYZNC/denZNC;
1773  }
1774  else{
1775  xyZNC[0] = xyZNC[1] = 999.;
1776  zncEnergy = 0.;
1777  }
1778  if(denZNA!=0) {
1779  xyZNA[0] = numXZNA/denZNA;
1780  xyZNA[1] = numYZNA/denZNA;
1781  }
1782  else{
1783  xyZNA[0] = xyZNA[1] = 999.;
1784  znaEnergy = 0.;
1785  }
1786  }
1787 
1788  if(!fAllChONZNC) denZNC=-1.;
1789  if(!fAllChONZNA) denZNA=-1.;
1790 
1791  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]);
1792  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]);
1793 
1794  fFlowEvent->SetZDC2Qsub(xyZNC,denZNC,xyZNA,denZNA);
1795 
1796  // ******************************************************************************
1797 
1798  Double_t tdcSum = aodZDC->GetZDCTimeSum();
1799  Double_t tdcDiff = aodZDC->GetZDCTimeDiff();
1800  fhDebunch->Fill(tdcDiff, tdcSum);
1801 
1802  for(int i=0; i<5; i++){
1803  fhZNCPM[i]->Fill(towZNC[i]);
1804  if((i<4) && (towZNC[0]>0.)) fhZNCPMQiPMC[i]->Fill(towZNC[i+1]/towZNC[0]);
1805  }
1806  for(int i=0; i<5; i++){
1807  fhZNAPM[i]->Fill(towZNA[i]);
1808  if(((i<4) && towZNA[0]>0.)) fhZNAPMQiPMC[i]->Fill(towZNA[i+1]/towZNA[0]);
1809  }
1810 
1811  Double_t energyZPC = (Double_t) (aodZDC->GetZPCEnergy());
1812  Double_t energyZPA = (Double_t) (aodZDC->GetZPAEnergy());
1813  Double_t energyZEM1 = (Double_t) (aodZDC->GetZEM1Energy());
1814  Double_t energyZEM2 = (Double_t) (aodZDC->GetZEM2Energy());
1815 
1816  fhZNCvsZNA->Fill(energyZNA, energyZNC);
1817  fhZDCCvsZDCCA->Fill(energyZNA+energyZPA, energyZNC+energyZPC);
1818  fhZNCvsZPC->Fill(energyZPC, energyZNC);
1819  fhZNAvsZPA->Fill(energyZPA, energyZNA);
1820  fhZNvsZP->Fill(energyZPA+energyZPC, energyZNA+energyZNC);
1821  fhZNvsVZERO->Fill(multV0A+multV0C, energyZNC+energyZNA);
1822  fhZDCvsVZERO->Fill(multV0A+multV0C, energyZNA+energyZPA+energyZNC+energyZPC);
1823  fhZDCvsTracklets->Fill((Double_t) (nTracklets), energyZNA+energyZPA+energyZNC+energyZPC);
1824 
1825  Double_t asymmetry = -999.;
1826  if((energyZNC+energyZNA)>0.) asymmetry = (energyZNC-energyZNA)/(energyZNC+energyZNA);
1827  fhAsymm->Fill(asymmetry);
1828  fhZNAvsAsymm->Fill(asymmetry, energyZNA/1000.);
1829  fhZNCvsAsymm->Fill(asymmetry, energyZNC/1000.);
1830 
1831  fhZNCvscentrality->Fill(centrperc, energyZNC/1000.);
1832  fhZNAvscentrality->Fill(centrperc, energyZNA/1000.);
1833 
1834  } // PHYSICS SELECTION
1835 
1836  }
1837 
1838  // p) cache run number
1840 
1841  // printf("debug: NoRPs %e, NoPOIs %e, RunNum %d, Cen %e \n",fFlowEvent->GetNumberOfRPs(),fFlowEvent->GetNumberOfPOIs(),fCachedRunNum,fFlowEvent->GetCentrality());
1842 
1843  PostData(1, fFlowEvent);
1844 
1845  PostData(2, fOutput);
1846 }
1847 //________________________________________________________________________
1848 
1850 {
1851  Bool_t BisPileup=kFALSE;
1852  Double_t centrV0M=300., centrCL1=300.;
1853 
1855 
1856  // pileup for LHC10h and LHC11h
1857 
1858  centrV0M = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
1859  centrCL1 = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP()->GetCentralityPercentile("CL1");
1860 
1861  // check anyway pileup
1862  if (plpMV(aod)) {
1863  fPileUpCount->Fill(0.5);
1864  BisPileup=kTRUE;
1865  }
1866 
1867  Short_t isPileup = aod->IsPileupFromSPD(3);
1868  if (isPileup != 0) {
1869  fPileUpCount->Fill(1.5);
1870  BisPileup=kTRUE;
1871  }
1872 
1873  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
1874  fPileUpCount->Fill(2.5);
1875  BisPileup=kTRUE;
1876  }
1877 
1878  if (aod->IsIncompleteDAQ()) {
1879  fPileUpCount->Fill(3.5);
1880  BisPileup=kTRUE;
1881  }
1882 
1883  // check vertex consistency
1884  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
1885  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
1886 
1887  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
1888  fPileUpCount->Fill(5.5);
1889  BisPileup=kTRUE;
1890  }
1891 
1892  double covTrc[6], covSPD[6];
1893  vtTrc->GetCovarianceMatrix(covTrc);
1894  vtSPD->GetCovarianceMatrix(covSPD);
1895 
1896  double dz = vtTrc->GetZ() - vtSPD->GetZ();
1897 
1898  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
1899  double errTrc = TMath::Sqrt(covTrc[5]);
1900  double nsigTot = dz/errTot;
1901  double nsigTrc = dz/errTrc;
1902 
1903  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
1904  fPileUpCount->Fill(6.5);
1905  BisPileup=kTRUE;
1906  }
1907 
1908  if (fAnalysisUtil->IsPileUpEvent(InputEvent())) {
1909  fPileUpCount->Fill(7.5);
1910  BisPileup=kTRUE;
1911  }
1912 
1913  }
1914  else {
1915 
1916  // pileup for LHC15o, using AliMultSelection
1917 
1918  if(fMultSelection) {
1919  centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
1920  centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
1921  } else {
1922  BisPileup=kTRUE;
1923  }
1924 
1925  // pileup from AliMultSelection
1926  if(!fMultSelection->GetThisEventIsNotPileup()) fPileUpMultSelCount->Fill(0.5);
1927  if(!fMultSelection->GetThisEventIsNotPileupMV()) fPileUpMultSelCount->Fill(1.5);
1928  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) fPileUpMultSelCount->Fill(2.5);
1929  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) fPileUpMultSelCount->Fill(3.5);
1930  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) fPileUpMultSelCount->Fill(4.5);
1931  if(!fMultSelection->GetThisEventIsNotAsymmetricInVZERO()) fPileUpMultSelCount->Fill(5.5);
1932  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) fPileUpMultSelCount->Fill(6.5);
1933  if(!fMultSelection->GetThisEventHasGoodVertex2016()) fPileUpMultSelCount->Fill(7.5);
1934 
1935  // pile-up a la Dobrin for LHC15o
1936  if (plpMV(aod)) {
1937  fPileUpCount->Fill(0.5);
1938  BisPileup=kTRUE;
1939  }
1940 
1941  Short_t isPileup = aod->IsPileupFromSPD(3);
1942  if (isPileup != 0) {
1943  fPileUpCount->Fill(1.5);
1944  BisPileup=kTRUE;
1945  }
1946 
1947  if (((AliAODHeader*)aod->GetHeader())->GetRefMultiplicityComb08() < 0) {
1948  fPileUpCount->Fill(2.5);
1949  BisPileup=kTRUE;
1950  }
1951 
1952  if (aod->IsIncompleteDAQ()) {
1953  fPileUpCount->Fill(3.5);
1954  BisPileup=kTRUE;
1955  }
1956 
1957  if(fabs(centrV0M-centrCL1)>7.5) {
1958  fPileUpCount->Fill(4.5);
1959  BisPileup=kTRUE;
1960  }
1961 
1962  // check vertex consistency
1963  const AliAODVertex* vtTrc = aod->GetPrimaryVertex();
1964  const AliAODVertex* vtSPD = aod->GetPrimaryVertexSPD();
1965 
1966  if (vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
1967  fPileUpCount->Fill(5.5);
1968  BisPileup=kTRUE;
1969  }
1970 
1971  double covTrc[6], covSPD[6];
1972  vtTrc->GetCovarianceMatrix(covTrc);
1973  vtSPD->GetCovarianceMatrix(covSPD);
1974 
1975  double dz = vtTrc->GetZ() - vtSPD->GetZ();
1976 
1977  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
1978  double errTrc = TMath::Sqrt(covTrc[5]);
1979  double nsigTot = dz/errTot;
1980  double nsigTrc = dz/errTrc;
1981 
1982  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
1983  fPileUpCount->Fill(6.5);
1984  BisPileup=kTRUE;
1985  }
1986 
1987  // cuts on tracks
1988  const Int_t nTracks = aod->GetNumberOfTracks();
1989  Int_t multEsd = ((AliAODHeader*)aod->GetHeader())->GetNumberOfESDTracks();
1990 
1991  Int_t multTrk = 0;
1992  Int_t multTrkBefC = 0;
1993  Int_t multTrkTOFBefC = 0;
1994  Int_t multTPC = 0;
1995 
1996  for (Int_t it = 0; it < nTracks; it++) {
1997 
1998  AliAODTrack* aodTrk = (AliAODTrack*)aod->GetTrack(it);
1999  if (!aodTrk){
2000  delete aodTrk;
2001  continue;
2002  }
2003 
2004 // if (aodTrk->TestFilterBit(32)){
2005 // multTrkBefC++;
2006 //
2007 // if ( TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10. && aodTrk->GetTOFsignal() >= 12000. && aodTrk->GetTOFsignal() <= 25000.)
2008 // multTrkTOFBefC++;
2009 //
2010 // if ((TMath::Abs(aodTrk->Eta()) < 0.8) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= 0.2) && (aodTrk->Pt() < 20.))
2011 // multTrk++;
2012 // }
2013 
2014  if (aodTrk->TestFilterBit(128))
2015  multTPC++;
2016 
2017  } // end of for (Int_t it = 0; it < nTracks; it++)
2018 
2019  Double_t multTPCn = multTPC;
2020  Double_t multEsdn = multEsd;
2021  Double_t multESDTPCDif = multEsdn - multTPCn*3.38;
2022 
2023  if (multESDTPCDif > (fRejectPileUpTight?700.:15000.)) {
2024  fPileUpCount->Fill(7.5);
2025  BisPileup=kTRUE;
2026  }
2027 
2028  if(fRejectPileUpTight) {
2029  if(BisPileup==kFALSE) {
2030  if(!fMultSelection->GetThisEventIsNotPileup()) BisPileup=kTRUE;
2031  if(!fMultSelection->GetThisEventIsNotPileupMV()) BisPileup=kTRUE;
2032  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) BisPileup=kTRUE;
2033  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) BisPileup=kTRUE;
2034  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) BisPileup=kTRUE;
2035  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) BisPileup=kTRUE;
2036  if(!fMultSelection->GetThisEventHasGoodVertex2016()) BisPileup=kTRUE;
2037  if(BisPileup) fPileUpCount->Fill(8.5);
2038  }
2039  }
2040  }
2041 
2042  return BisPileup;
2043 }
2044 
2045 //________________________________________________________________________
2046 
2048 {
2049  Double_t EtC = BadTowerCalibHist->ProjectionY("",BadTowerCalibHist->GetXaxis()->FindBin(Et),BadTowerCalibHist->GetXaxis()->FindBin(Et))->GetRandom();
2050  return EtC;
2051 }
2052 
2053 //________________________________________________________________________
2054 
2056 {
2057  Int_t CenBin=-1;
2058  if (Centrality>0. && Centrality<5.) CenBin=0;
2059  if (Centrality>5. && Centrality<10.) CenBin=1;
2060  if (Centrality>10. && Centrality<20.) CenBin=2;
2061  if (Centrality>20. && Centrality<30.) CenBin=3;
2062  if (Centrality>30. && Centrality<40.) CenBin=4;
2063  if (Centrality>40. && Centrality<50.) CenBin=5;
2064  if (Centrality>50. && Centrality<60.) CenBin=6;
2065  if (Centrality>60. && Centrality<70.) CenBin=7;
2066  if (Centrality>70. && Centrality<80.) CenBin=8;
2067  if (Centrality>80. && Centrality<90.) CenBin=9;
2068  if (CenBin>=fnCen) CenBin=-1;
2069  if (fnCen==1) CenBin=0;
2070  return CenBin;
2071 } // end of AliFlowAnalysisCRC::GetCRCCenBin(Double_t Centrality)
2072 //_____________________________________________________________________________
2073 
2074 Double_t AliAnalysisTaskCRCZDC::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
2075 {
2076  // calculate sqrt of weighted distance to other vertex
2077  if (!v0 || !v1) {
2078  printf("One of vertices is not valid\n");
2079  return 0;
2080  }
2081  static TMatrixDSym vVb(3);
2082  double dist = -1;
2083  double dx = v0->GetX()-v1->GetX();
2084  double dy = v0->GetY()-v1->GetY();
2085  double dz = v0->GetZ()-v1->GetZ();
2086  double cov0[6],cov1[6];
2087  v0->GetCovarianceMatrix(cov0);
2088  v1->GetCovarianceMatrix(cov1);
2089  vVb(0,0) = cov0[0]+cov1[0];
2090  vVb(1,1) = cov0[2]+cov1[2];
2091  vVb(2,2) = cov0[5]+cov1[5];
2092  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
2093  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
2094  vVb.InvertFast();
2095  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
2096  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
2097  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
2098  return dist>0 ? TMath::Sqrt(dist) : -1;
2099 }
2100 //________________________________________________________________________
2101 
2103 {
2104  // check for multi-vertexer pile-up
2105 
2106  const int kMinPlpContrib = 5;
2107  const double kMaxPlpChi2 = 5.0;
2108  const double kMinWDist = 15;
2109 
2110  const AliVVertex* vtPrm = 0;
2111  const AliVVertex* vtPlp = 0;
2112  int nPlp = 0;
2113 
2114  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
2115  vtPrm = aod->GetPrimaryVertex();
2116  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
2117 
2118  //int bcPrim = vtPrm->GetBC();
2119 
2120  for (int ipl=0;ipl<nPlp;ipl++) {
2121  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
2122  //
2123  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
2124  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
2125  // int bcPlp = vtPlp->GetBC();
2126  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
2127  //
2128  double wDst = GetWDist(vtPrm,vtPlp);
2129  if (wDst<kMinWDist) continue;
2130  //
2131  return kTRUE; // pile-up: well separated vertices
2132  }
2133 
2134  return kFALSE;
2135 }
2136 
2137 //________________________________________________________________________
2139  fCutContainer->AddAt(cutsRP,0); fCutsRP=cutsRP; cutsRP->SetPOItype(0);
2140 }
2141 
2142 //________________________________________________________________________
2144  fCutContainer->AddAt(cutsPOI,1); fCutsPOI=cutsPOI; cutsPOI->SetPOItype(1);
2145 }
2146 
2147 //________________________________________________________________________
2149 {
2150  // Terminate analysis
2151  //
2152  /* if(fDebug > 1) printf(" **** AliAnalysisTaskCRCZDC::Terminate() \n");
2153 
2154  //fOutput = dynamic_cast<TList*> (GetOutputData(1));
2155  //if(!fOutput) printf("ERROR: fOutput not available\n");
2156  */
2157 }
void SetCharge(Int_t charge)
TH2F * fhZDCCvsZDCCA
ZNC vs ZNA;.
static const Int_t fkVZEROnQAplots
TH2F * fhZNvsVZERO
ZNC+ZNA vs ZPC+ZPA;.
TH1F * fPtSpecFB128[2][10]
PtSpecRec FB96.
void SetPOItype(Int_t poiType, Bool_t b=kTRUE)
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.
void SetZNCQ0(Double_t const en)
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)
void SetZNAQ0(Double_t const en)
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 AddTrack(AliFlowTrackSimple *track)
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
void SetAbsOrbit(UInt_t const en)
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