AliPhysics  3bba2fe (3bba2fe)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskJetV2.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 the 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  * Jet V2 task
18  *
19  * this task is part of the emcal jet framework and should be run in the emcaljet train
20  * the following extensions to an accepted AliVEvent are expected:
21  * - (anti-kt) jets
22  * - background estimate rho
23  * - pico tracks
24  * aod's and esd's are handled transparently
25  * the task will attempt to estimate a phi-dependent background density rho
26  * by fitting vn harmonics to the dpt/dphi distribution
27  *
28  * author: Redmer Alexander Bertens, Utrecht Univeristy, Utrecht, Netherlands
29  * rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
30  */
31 
32 // root includes
33 #include <TStyle.h>
34 #include <TRandom3.h>
35 #include <TChain.h>
36 #include <TMath.h>
37 #include <TF1.h>
38 #include <TF2.h>
39 #include <TH1F.h>
40 #include <TH2F.h>
41 #include <TH3F.h>
42 #include <TProfile.h>
43 #include <TFile.h>
44 // aliroot includes
45 #include <AliAnalysisTask.h>
46 #include <AliAnalysisManager.h>
47 #include <AliCentrality.h>
48 #include <AliVVertex.h>
49 #include <AliVTrack.h>
50 #include <AliVVZERO.h>
51 #include <AliESDEvent.h>
52 #include <AliAODEvent.h>
53 #include <AliAODTrack.h>
54 #include <AliOADBContainer.h>
55 //#include <AliMultSelection.h>
56 #include <AliInputEventHandler.h>
57 // emcal jet framework includes
58 #include <AliPicoTrack.h>
59 #include <AliEmcalJet.h>
60 #include <AliRhoParameter.h>
61 #include <AliLocalRhoParameter.h>
62 #include <AliAnalysisTaskJetV2.h>
63 #include <AliClusterContainer.h>
64 
66 using namespace std;
67 
69 
71  fRunToyMC(kFALSE), fLocalInit(0), fAttachToEvent(kTRUE), fFillHistograms(kTRUE), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(-1.), fReduceBinsYByFactor(-1.), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fExpectedRuns(0), fExpectedSemiGoodRuns(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fAcceptanceWeights(kFALSE), fEventPlaneWeight(1.), fTracksCont(0), fClusterCont(0), fJetsCont(0), fLeadingJet(0), fLeadingJetAfterSub(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fFitGoodnessTest(kChi2Poisson), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fUsePtWeightErrorPropagation(kTRUE), fUse2DIntegration(kFALSE), fDetectorType(kVZEROComb), fAnalysisType(kCharged), fFitModulationOptions("QWLI"), fRunModeType(kGrid), fDataType(kESD), fCollisionType(kPbPb), fRandom(0), fRunNumber(-1), fRunNumberCaliInfo(-1), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fFitControl(0), fMinPvalue(0.01), fMaxPvalue(1), fNameSmallRho(""), fCachedRho(0), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.), fSemiGoodJetMinPhi(0.), fSemiGoodJetMaxPhi(4.), fSemiGoodTrackMinPhi(0.), fSemiGoodTrackMaxPhi(4.), fHistCentrality(0), fHistCentralityPercIn(0), fHistCentralityPercOut(0), fHistCentralityPercLost(0), fHistVertexz(0), fHistMultCorAfterCuts(0), fHistMultvsCentr(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistRunnumbersCaliInfo(0), fHistPvalueCDFROOT(0), fHistPvalueCDFROOTCent(0), fHistChi2ROOTCent(0), fHistPChi2Root(0), fHistPvalueCDF(0), fHistPvalueCDFCent(0), fHistChi2Cent(0), fHistPChi2(0), fHistKolmogorovTest(0), fHistKolmogorovTestCent(0), fHistPKolmogorov(0), fHistRhoStatusCent(0), fHistUndeterminedRunQA(0), fMinDisanceRCtoLJ(0), fMaxCones(-1), fExcludeLeadingJetsFromFit(1.), fExcludeJetsWithTrackPt(9999.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiVZEROAV0M(0), fHistPsiVZEROCV0M(0), fHistPsiVZEROVV0M(0), fHistPsiTPCV0M(0), fHistPsiVZEROATRK(0), fHistPsiVZEROCTRK(0), fHistPsiVZEROTRK(0), fHistPsiTPCTRK(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0), fVZEROgainEqualization(0x0), fVZEROApol(0), fVZEROCpol(0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0), fSigma2A(0x0), fSigma2C(0x0), fSigma3A(0x0), fSigma3C(0x0), fWeightForVZERO(kChi), fOADB(0x0), fHistQxV0aBC(0x0), fHistQyV0aBC(0x0), fHistQxV0cBC(0x0), fHistQyV0cBC(0x0), fHistQxV0a(0x0), fHistQyV0a(0x0), fHistQxV0c(0x0), fHistQyV0c(0x0), fHistMultVsCellBC(0x0), fHistMultVsCell(0x0), fHistEPBC(0x0), fHistEP(0x0)
72 {
73  for(Int_t i(0); i < 10; i++) {
74  fEventPlaneWeights[i] = 0;
75  fProfV2Resolution[i] = 0;
76  fProfV3Resolution[i] = 0;
77  fHistPicoTrackPt[i] = 0;
78  fHistPicoTrackMult[i] = 0;
79  fHistPicoCat1[i] = 0;
80  fHistPicoCat2[i] = 0;
81  fHistPicoCat3[i] = 0;
82  fHistClusterPt[i] = 0;
83  fHistClusterEtaPhi[i] = 0;
85  fHistTriggerQAIn[i] = 0;
86  fHistTriggerQAOut[i] = 0;
87  fHistEPCorrelations[i] = 0;
88  fHistEPCorrAvChi[i] = 0;
89  fHistEPCorrAvSigma[i] = 0;
90  fHistEPCorrChiSigma[i] = 0;
93  fHistPsiTPCLeadingJet[i] = 0;
97  fHistPsi2Correlation[i] = 0;
99  fHistRhoPackage[i] = 0;
100  fHistRho[i] = 0;
101  fHistRhoEtaBC[i] = 0;
102  fHistRCPhiEta[i] = 0;
103  fHistRhoVsRCPt[i] = 0;
104  fHistRCPt[i] = 0;
105  fHistDeltaPtDeltaPhi2[i] = 0;
107  fHistRCPhiEtaExLJ[i] = 0;
108  fHistRhoVsRCPtExLJ[i] = 0;
109  fHistRCPtExLJ[i] = 0;
112  fHistJetPtRaw[i] = 0;
113  fHistJetPt[i] = 0;
114  fHistJetPtBC[i] = 0;
115  fHistJetEtaPhi[i] = 0;
116  fHistJetEtaPhiBC[i] = 0;
117  fHistJetPtArea[i] = 0;
118  fHistJetPtAreaBC[i] = 0;
119  fHistJetPtEta[i] = 0;
120  fHistJetPtConstituents[i] = 0;
121  fHistJetEtaRho[i] = 0;
122  fHistJetPsi2Pt[i] = 0;
123  fHistJetLJPsi2Pt[i] = 0;
124  fHistJetLJPsi2PtRatio[i] = 0;
125  fHistJetPsi2PtRho0[i] = 0;
126  }
127  for(Int_t i(0); i < 9; i++) {
128  for(Int_t j(0); j < 2; j++) {
129  for(Int_t k(0); k < 2; k++) {
130  fMeanQ[i][j][k] = 0.;
131  fWidthQ[i][j][k] = 0.;
132  fMeanQv3[i][j][k] = 0.;
133  fWidthQv3[i][j][k] = 0.;
134  }
135  }
136  }
137  // default constructor
138 }
139 //_____________________________________________________________________________
140 AliAnalysisTaskJetV2::AliAnalysisTaskJetV2(const char* name, runModeType type, Bool_t baseClassHistos) : AliAnalysisTaskEmcalJet(name, baseClassHistos),
141  fRunToyMC(kFALSE), fLocalInit(0), fAttachToEvent(kTRUE), fFillHistograms(kTRUE), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(-1.), fReduceBinsYByFactor(-1.), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fExpectedRuns(0), fExpectedSemiGoodRuns(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fAcceptanceWeights(kFALSE), fEventPlaneWeight(1.), fTracksCont(0), fClusterCont(0), fJetsCont(0), fLeadingJet(0), fLeadingJetAfterSub(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fFitGoodnessTest(kChi2Poisson), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fUsePtWeightErrorPropagation(kTRUE), fUse2DIntegration(kFALSE), fDetectorType(kVZEROComb), fAnalysisType(kCharged), fFitModulationOptions("QWLI"), fRunModeType(type), fDataType(kESD), fCollisionType(kPbPb), fRandom(0), fRunNumber(-1), fRunNumberCaliInfo(-1), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fFitControl(0), fMinPvalue(0.01), fMaxPvalue(1), fNameSmallRho(""), fCachedRho(0), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.), fSemiGoodJetMinPhi(0.), fSemiGoodJetMaxPhi(4.), fSemiGoodTrackMinPhi(0.), fSemiGoodTrackMaxPhi(4.), fHistCentrality(0), fHistCentralityPercIn(0), fHistCentralityPercOut(0), fHistCentralityPercLost(0), fHistVertexz(0), fHistMultCorAfterCuts(0), fHistMultvsCentr(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistRunnumbersCaliInfo(0), fHistPvalueCDFROOT(0), fHistPvalueCDFROOTCent(0), fHistChi2ROOTCent(0), fHistPChi2Root(0), fHistPvalueCDF(0), fHistPvalueCDFCent(0), fHistChi2Cent(0), fHistPChi2(0), fHistKolmogorovTest(0), fHistKolmogorovTestCent(0), fHistPKolmogorov(0), fHistRhoStatusCent(0), fHistUndeterminedRunQA(0), fMinDisanceRCtoLJ(0), fMaxCones(-1), fExcludeLeadingJetsFromFit(1.), fExcludeJetsWithTrackPt(9999), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiVZEROAV0M(0), fHistPsiVZEROCV0M(0), fHistPsiVZEROVV0M(0), fHistPsiTPCV0M(0), fHistPsiVZEROATRK(0), fHistPsiVZEROCTRK(0), fHistPsiVZEROTRK(0), fHistPsiTPCTRK(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0), fVZEROgainEqualization(0x0), fVZEROApol(0), fVZEROCpol(0), fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0), fSigma2A(0x0), fSigma2C(0x0), fSigma3A(0x0), fSigma3C(0x0), fWeightForVZERO(kChi), fOADB(0x0), fHistQxV0aBC(0x0), fHistQyV0aBC(0x0), fHistQxV0cBC(0x0), fHistQyV0cBC(0x0), fHistQxV0a(0x0), fHistQyV0a(0x0), fHistQxV0c(0x0), fHistQyV0c(0x0), fHistMultVsCellBC(0x0), fHistMultVsCell(0x0), fHistEPBC(0x0), fHistEP(0x0)
142 {
143  for(Int_t i(0); i < 10; i++) {
144  fEventPlaneWeights[i] = 0;
145  fProfV2Resolution[i] = 0;
146  fProfV3Resolution[i] = 0;
147  fHistPicoTrackPt[i] = 0;
148  fHistPicoTrackMult[i] = 0;
149  fHistPicoCat1[i] = 0;
150  fHistPicoCat2[i] = 0;
151  fHistPicoCat3[i] = 0;
152  fHistClusterPt[i] = 0;
153  fHistClusterEtaPhi[i] = 0;
155  fHistTriggerQAIn[i] = 0;
156  fHistTriggerQAOut[i] = 0;
157  fHistEPCorrelations[i] = 0;
158  fHistEPCorrAvChi[i] = 0;
159  fHistEPCorrAvSigma[i] = 0;
160  fHistEPCorrChiSigma[i] = 0;
163  fHistPsiTPCLeadingJet[i] = 0;
164  fHistPsiVZEROALeadingJet[i] = 0;
167  fHistPsi2Correlation[i] = 0;
169  fHistRhoPackage[i] = 0;
170  fHistRho[i] = 0;
171  fHistRhoEtaBC[i] = 0;
172  fHistRCPhiEta[i] = 0;
173  fHistRhoVsRCPt[i] = 0;
174  fHistRCPt[i] = 0;
175  fHistDeltaPtDeltaPhi2[i] = 0;
177  fHistRCPhiEtaExLJ[i] = 0;
178  fHistRhoVsRCPtExLJ[i] = 0;
179  fHistRCPtExLJ[i] = 0;
182  fHistJetPtRaw[i] = 0;
183  fHistJetPt[i] = 0;
184  fHistJetPtBC[i] = 0;
185  fHistJetEtaPhi[i] = 0;
186  fHistJetEtaPhiBC[i] = 0;
187  fHistJetPtArea[i] = 0;
188  fHistJetPtAreaBC[i] = 0;
189  fHistJetPtEta[i] = 0;
190  fHistJetPtConstituents[i] = 0;
191  fHistJetEtaRho[i] = 0;
192  fHistJetPsi2Pt[i] = 0;
193  fHistJetLJPsi2Pt[i] = 0;
194  fHistJetLJPsi2PtRatio[i] = 0;
195  fHistJetPsi2PtRho0[i] = 0;
196  }
197  for(Int_t i(0); i < 9; i++) {
198  for(Int_t j(0); j < 2; j++) {
199  for(Int_t k(0); k < 2; k++) {
200  fMeanQ[i][j][k] = 0.;
201  fWidthQ[i][j][k] = 0.;
202  fMeanQv3[i][j][k] = 0.;
203  fWidthQv3[i][j][k] = 0.;
204  }
205  }
206  }
207 
208  // constructor
209  DefineInput(0, TChain::Class());
210  Int_t startAt(1);
211  if(fCreateHisto) startAt++;
212  DefineOutput(startAt, TList::Class());
213  switch (fRunModeType) {
214  case kLocal : {
215  gStyle->SetOptFit(1);
216  DefineOutput(startAt+1, TList::Class());
217  DefineOutput(startAt+2, TList::Class());
218  } break;
219  default: break;
220  }
221  switch (fCollisionType) {
222  case kPythia : {
224  } break;
225  default : break;
226  }
227  if(fLocalRhoName=="") fLocalRhoName = Form("LocalRhoFrom_%s", GetName());
228  SetMakeGeneralHistograms(baseClassHistos);
229 }
230 //_____________________________________________________________________________
232 {
233  // destructor
234  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
235  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
236  #endif
237 
238  if(fOutputList) {delete fOutputList; fOutputList = 0x0;}
239  if(fOutputListGood) {delete fOutputListGood; fOutputListGood = 0x0;}
240  if(fOutputListBad) {delete fOutputListBad; fOutputListBad = 0x0;}
241  if(fFitModulation) {delete fFitModulation; fFitModulation = 0x0;}
242  if(fHistSwap) {delete fHistSwap; fHistSwap = 0x0;}
244  if(fExpectedRuns) {delete fExpectedRuns; fExpectedRuns = 0x0;}
246  if(fFitControl) {delete fFitControl; fFitControl = 0x0;}
248  if(fChi2A) {delete fChi2A; fChi2A = 0x0;}
249  if(fChi2C) {delete fChi2C; fChi2C = 0x0;}
250  if(fChi3A) {delete fChi3A; fChi3A = 0x0;}
251  if(fChi3C) {delete fChi3C; fChi3C = 0x0;}
252  if(fSigma2A) {delete fSigma2A; fSigma2A = 0x0;}
253  if(fSigma2C) {delete fSigma2C; fSigma2C = 0x0;}
254  if(fSigma3A) {delete fSigma3A; fSigma3A = 0x0;}
255  if(fSigma3C) {delete fSigma3C; fSigma3C = 0x0;}
256  if(fOADB && !fOADB->IsZombie()) {
257  fOADB->Close(); fOADB = 0x0;
258  } else if (fOADB) fOADB = 0x0;
259 }
260 //_____________________________________________________________________________
262 {
263  // Init the analysis
264  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
265  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
266  #endif
267  if(!fLocalRho) {
268  fLocalRho = new AliLocalRhoParameter(fLocalRhoName.Data(), 0);
269  if(fAttachToEvent) {
270  if(!(InputEvent()->FindListObject(fLocalRho->GetName()))) {
271  InputEvent()->AddObject(fLocalRho);
272  } else {
273  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), fLocalRho->GetName()));
274  }
275  }
276  }
277  AliAnalysisTaskEmcalJet::ExecOnce(); // init the base class
278  if(!GetJetContainer()) AliFatal(Form("%s: Couldn't find jet container. Aborting !", GetName()));
279 }
280 //_____________________________________________________________________________
282 {
283  // determine the run number to see if the track and jet cuts should be refreshed for semi-good TPC runs
284  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
285  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
286  #endif
287  if(fRunNumber != InputEvent()->GetRunNumber()) {
288  fRunNumber = InputEvent()->GetRunNumber(); // set the current run number
289  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
290  printf("__FUNC__ %s > NEW RUNNUMBER DETECTED \n ", __func__);
291  #endif
292  // check if this is 10h or 11h data
293  switch (fCollisionType) {
294  case kPbPb10h : {
295  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
296  printf(" LHC10h data, assuming full acceptance, reading VZERO calibration DB \n ");
297  #endif
298  // for 10h data the vzero event plane calibration needs to be cached
300  // no need to change rho or acceptance for 10h, so we're done
301  return kTRUE;
302  } break;
303  case kJetFlowMC : {
304  return kTRUE;
305  } break;
306  default : {
307  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
308  printf(" checking runnumber to adjust acceptance on the fly \n");
309  #endif
311  } break;
312  }
313  // reset the cuts. should be a pointless operation except for the case where the run number changes
314  // from semi-good back to good on one node, which is not a likely scenario (unless trains will
315  // run as one masterjob)
316  switch (fAnalysisType) {
317  case kCharged: {
319  } break;
320  case kFull: {
322  } break;
323  default: {
325  } break;
326  }
327  if(fCachedRho) { // if there's a cached rho, it's the default, so switch back
328  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
329  printf("__FUNC__ %s > replacing rho with cached rho \n ", __func__);
330  #endif
331  fRho = fCachedRho; // reset rho back to cached value. again, should be pointless
332  }
333  Bool_t flaggedAsSemiGood(kFALSE); // not flagged as anything
334  for(Int_t i(0); i < fExpectedSemiGoodRuns->GetSize(); i++) {
335  if(fExpectedSemiGoodRuns->At(i) == fRunNumber) { // run is semi-good
336  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
337  printf("__FUNC__ %s > semi-good tpc run detected, adjusting acceptance \n ", __func__);
338  #endif
339  flaggedAsSemiGood = kTRUE;
340  switch (fAnalysisType) {
341  // for full jets the jet acceptance does not have to be changed as emcal does not
342  // cover the tpc low voltage readout strips
343  case kCharged: {
344  AliAnalysisTaskEmcalJet::SetJetPhiLimits(fSemiGoodJetMinPhi, fSemiGoodJetMaxPhi); // just an acceptance cut, jets are obtained from full azimuth, so no edge effects
345  } break;
346  default: break;
347  }
348  AliAnalysisTaskEmcal::SetTrackPhiLimits(fSemiGoodTrackMinPhi, fSemiGoodTrackMaxPhi); // only affects vn extraction, NOT jet finding
349  // for semi-good runs, also try to get the 'small rho' estimate, if it is available
350  AliRhoParameter* tempRho(dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fNameSmallRho.Data())));
351  if(tempRho) {
352  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
353  printf("__FUNC__ %s > switching to small rho, caching normal rho \n ", __func__);
354  #endif
355  fHistAnalysisSummary->SetBinContent(54, 1.); // bookkeep the fact that small rho is used
356  fCachedRho = fRho; // cache the original rho ...
357  fRho = tempRho; // ... and use the small rho
358  }
359  }
360  }
361  if(!flaggedAsSemiGood) {
362  // in case the run is not a semi-good run, check if it is recognized as another run
363  // only done to catch unexpected runs
364  for(Int_t i(0); i < fExpectedRuns->GetSize(); i++) {
365  if(fExpectedRuns->At(i) == fRunNumber) break; // run is known, break the loop else store the number in a random bin
366  fHistUndeterminedRunQA->SetBinContent(TMath::Nint(10.*gRandom->Uniform(0.,.9))+1, fRunNumber);
367  }
368  fHistAnalysisSummary->SetBinContent(53, 1.); // bookkeep which rho estimate is used
369  }
370  }
371  return kTRUE;
372 }
373 //_____________________________________________________________________________
375 {
376  // initialize the anaysis
377  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
378  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
379  #endif
380  // if not set, estimate the number of cones that would fit into the selected acceptance
381  if(fMaxCones <= 0) fMaxCones = TMath::CeilNint((TMath::Abs(GetJetContainer()->GetJetEtaMax()-GetJetContainer()->GetJetEtaMin())*TMath::Abs(GetJetContainer()->GetJetPhiMax()-GetJetContainer()->GetJetPhiMin()))/(TMath::Pi()*GetJetRadius()*GetJetRadius()));
382  // manually 'override' the default acceptance cuts of the emcal framework (use with caution)
384  if(dynamic_cast<AliAODEvent*>(InputEvent())) fDataType = kAOD; // determine the datatype
385  else if(dynamic_cast<AliESDEvent*>(InputEvent())) fDataType = kESD;
386  fHistAnalysisSummary->SetBinContent(36, (int)fDataType);
387  if(!fRandom) fRandom = new TRandom3(0); // set randomizer and random seed
388  switch (fFitModulationType) {
389  case kNoFit : { SetModulationFit(new TF1("fix_kNoFit", "[0]", 0, TMath::TwoPi())); } break;
390  case kV2 : {
391  SetModulationFit(new TF1("fit_kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
392  fFitModulation->SetParameter(0, 0.); // normalization
393  fFitModulation->SetParameter(3, 0.2); // v2
394  fFitModulation->FixParameter(1, 1.); // constant
395  fFitModulation->FixParameter(2, 2.); // constant
396  } break;
397  case kV3: {
398  SetModulationFit(new TF1("fit_kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
399  fFitModulation->SetParameter(0, 0.); // normalization
400  fFitModulation->SetParameter(3, 0.2); // v3
401  fFitModulation->FixParameter(1, 1.); // constant
402  fFitModulation->FixParameter(2, 3.); // constant
403  } break;
404  default : { // for the combined fit, the 'direct fourier series' or the user supplied vn values we use v2 and v3
405  SetModulationFit(new TF1("fit_kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi()));
406  fFitModulation->SetParameter(0, 0.); // normalization
407  fFitModulation->SetParameter(3, 0.2); // v2
408  fFitModulation->FixParameter(1, 1.); // constant
409  fFitModulation->FixParameter(2, 2.); // constant
410  fFitModulation->FixParameter(5, 3.); // constant
411  fFitModulation->SetParameter(7, 0.2); // v3
412  } break;
413  }
414  switch (fRunModeType) {
415  case kGrid : { fFitModulationOptions += "N0"; } break;
416  default : break;
417  }
419  return kTRUE;
420 }
421 //_____________________________________________________________________________
422 TH1F* AliAnalysisTaskJetV2::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)
423 {
424  // book a TH1F and connect it to the output container
425  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
426  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
427  #endif
428  if(fReduceBinsXByFactor > 0 ) bins = TMath::Nint(bins/fReduceBinsXByFactor);
429  if(append && !fOutputList) return 0x0;
430  TString title(name);
431  if(c!=-1) { // format centrality dependent histograms accordingly
432  name = Form("%s_%i", name, c);
433  title += Form("_%i-%i", (int)(fCentralityClasses->At(c)), (int)(fCentralityClasses->At((1+c))));
434  }
435  title += Form(";%s;[counts]", x);
436  TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
437  histogram->Sumw2();
438  if(append) fOutputList->Add(histogram);
439  return histogram;
440 }
441 //_____________________________________________________________________________
442 TH2F* AliAnalysisTaskJetV2::BookTH2F(const char* name, const char* x, const char* y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t c, Bool_t append)
443 {
444  // book a TH2F and connect it to the output container
445  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
446  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
447  #endif
448  if(fReduceBinsXByFactor > 0 ) binsx = TMath::Nint(binsx/fReduceBinsXByFactor);
449  if(fReduceBinsYByFactor > 0 ) binsy = TMath::Nint(binsy/fReduceBinsYByFactor);
450  if(append && !fOutputList) return 0x0;
451  TString title(name);
452  if(c!=-1) { // format centrality dependent histograms accordingly
453  name = Form("%s_%i", name, c);
454  title += Form("_%i-%i", (int)fCentralityClasses->At(c), (int)(fCentralityClasses->At((1+c))));
455  }
456  title += Form(";%s;%s", x, y);
457  TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
458  histogram->Sumw2();
459  if(append) fOutputList->Add(histogram);
460  return histogram;
461 }
462 //_____________________________________________________________________________
463 TH3F* AliAnalysisTaskJetV2::BookTH3F(const char* name, const char* x, const char* y, const char* z, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t binsz, Double_t minz, Double_t maxz, Int_t c, Bool_t append)
464 {
465  // book a TH2F and connect it to the output container
466  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
467  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
468  #endif
469  if(fReduceBinsXByFactor > 0 ) {
470  binsx = TMath::Nint(binsx/fReduceBinsXByFactor);
471  binsy = TMath::Nint(binsy/fReduceBinsXByFactor);
472  binsz = TMath::Nint(binsz/fReduceBinsXByFactor);
473  }
474  if(append && !fOutputList) return 0x0;
475  TString title(name);
476  if(c!=-1) { // format centrality dependent histograms accordingly
477  name = Form("%s_%i", name, c);
478  title += Form("_%i-%i", (int)fCentralityClasses->At(c), (int)(fCentralityClasses->At((1+c))));
479  }
480  title += Form(";%s;%s;%s", x, y, z);
481  TH3F* histogram = new TH3F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy, binsz, minz, maxz);
482  histogram->Sumw2();
483  if(append) fOutputList->Add(histogram);
484  return histogram;
485 }
486 //_____________________________________________________________________________
488 {
489  // create output objects. also initializes some default values in case they aren't
490  // loaded via the AddTask macro
491  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
492  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
493  #endif
495  fOutputList = new TList();
496  fOutputList->SetOwner(kTRUE);
497  if(!fCentralityClasses) { // classes must be defined at this point
498  Double_t c[] = {0., 20., 40., 60., 80., 100.};
499  fCentralityClasses = new TArrayD(sizeof(c)/sizeof(c[0]), c);
500  }
501  if(!fExpectedRuns) { // expected runs must be defined at this point
502  Int_t r[] = {167813, 167988, 168066, 168068, 168069, 168076, 168104, 168212, 168311, 168322, 168325, 168341, 168361, 168362, 168458, 168460, 168461, 168992, 169091, 169094, 169138, 169143, 169167, 169417, 169835, 169837, 169838, 169846, 169855, 169858, 169859, 169923, 169956, 170027, 170036, 170081, /* up till here original good TPC list */169975, 169981, 170038, 170040, 170083, 170084, 170085, 170088, 170089, 170091, 170152, 170155, 170159, 170163, 170193, 170195, 170203, 170204, 170205, 170228, 170230, 170264, 170268, 170269, 170270, 170306, 170308, 170309, /* original semi-good tpc list */169415, 169411, 169035, 168988, 168984, 168826, 168777, 168512, 168511, 168467, 168464, 168342, 168310, 168115, 168108, 168107, 167987, 167915, 167903, /*new runs, good according to RCT */ 169238, 169160, 169156, 169148, 169145, 169144 /* run swith missing OROC 8 but seem ok in QA */};
503  fExpectedRuns = new TArrayI(sizeof(r)/sizeof(r[0]), r);
504  }
505  // set default semi-good runs only for 11h data
506  switch (fCollisionType) {
507  case kPbPb10h : {
508  fHistMultCorAfterCuts = new TH2F("fHistMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
510  fHistMultvsCentr = new TH2F("fHistMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 100.5, 101, 0, 3000);
512  } break;
513  default : {
514  if(!fExpectedSemiGoodRuns) {
515  Int_t r[] = {169975, 169981, 170038, 170040, 170083, 170084, 170085, 170088, 170089, 170091, 170152, 170155, 170159, 170163, 170193, 170195, 170203, 170204, 170205, 170228, 170230, 170264, 170268, 170269, 170270, 170306, 170308, 170309};
516  fExpectedSemiGoodRuns = new TArrayI(sizeof(r)/sizeof(r[0]), r);
517  }
518  }
519  }
520 
521  // global QA
522  fHistCentrality = BookTH1F("fHistCentrality", "centrality", 102, -2, 100);
523  fHistVertexz = BookTH1F("fHistVertexz", "vertex z (cm)", 100, -12, 12);
524  if(fAcceptanceWeights) {
525  fHistCentralityPercIn = new TProfile("fHistCentralityPercIn", "fHistCentralityPercIn", 102, -2, 100);
526  fHistCentralityPercOut = new TProfile("fHistCentralityPercOut", "fHistCentralityPercOut", 102, -2, 100);
527  fHistCentralityPercLost = new TProfile("fHistCentralityPercLost", "fHistCentralityPercLost", 102, -2, 100);
528  }
529  // for some histograms adjust the bounds according to analysis acceptance
530  Double_t etaMin(-1.), etaMax(1.), phiMin(0.), phiMax(TMath::TwoPi());
531  switch (fAnalysisType) {
532  case kFull : {
533  etaMin = -.7;
534  etaMax = .7;
535  phiMin = 1.405;
536  phiMax = 3.135;
537  } break;
538  default : break;
539  }
540 
541  // pico track and emcal cluster kinematics, trigger qa
542  for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
543  fHistPicoTrackPt[i] = BookTH1F("fHistPicoTrackPt", "p_{t} [GeV/c]", 100, 0, 100, i);
544  fHistPicoTrackMult[i] = BookTH1F("fHistPicoTrackMult", "multiplicity", 100, 0, 5000, i);
545  if(fFillQAHistograms) {
546  fHistPicoCat1[i] = BookTH2F("fHistPicoCat1", "#eta", "#phi", 50, -1, 1, 50, 0, TMath::TwoPi(), i);
547  fHistPicoCat2[i] = BookTH2F("fHistPicoCat2", "#eta", "#phi", 50, -1, 1, 50, 0, TMath::TwoPi(), i);
548  fHistPicoCat3[i] = BookTH2F("fHistPicoCat3", "#eta", "#phi", 50, -1, 1, 50, 0, TMath::TwoPi(), i);
550  fHistClusterPt[i] = BookTH1F("fHistClusterPt", "p_{t} [GeV/c]", 100, 0, 100, i);
551  fHistClusterEtaPhi[i] = BookTH2F("fHistClusterEtaPhi", "#eta", "#phi", 100, etaMax, etaMax, 100, phiMin, phiMax, i);
552  fHistClusterEtaPhiWeighted[i] = BookTH2F("fHistClusterEtaPhiWeighted", "#eta", "#phi", 100, etaMin, etaMax, 100, phiMin, phiMax, i);
553  }
554  fHistPsiTPCLeadingJet[i] = BookTH3F("fHistPsiTPCLeadingJet", "p_{t} [GeV/c]", "#Psi_{TPC}", "#varphi_{jet}", 70, 0, 210, 50, -1.*TMath::Pi()/2., TMath::Pi()/2., 50, phiMin, phiMax, i);
555  fHistEPCorrelations[i] = BookTH3F("fHistEPCorrelations", "EP_V0 average", "EP_V0 #chi", "EP_V0 #sigma", 50, -TMath::Pi()/2., TMath::Pi()/2., 50, -TMath::Pi()/2., TMath::Pi()/2., 50, -TMath::Pi()/2., TMath::Pi()/2.);
556  fHistEPCorrAvChi[i] = BookTH2F("fHistEPCorrAvChi", "EP_V0 average", "EP_V0 #chi", 50, -TMath::Pi()/2., TMath::Pi()/2., 50, -TMath::Pi()/2., TMath::Pi()/2., i);
557  fHistEPCorrAvSigma[i] = BookTH2F("fHistEPCorrAvSigma", "EP_V0 average", "EP_V0 #sigma", 50, -TMath::Pi()/2., TMath::Pi()/2., 50, -TMath::Pi()/2., TMath::Pi()/2., i);
558  fHistEPCorrChiSigma[i] = BookTH2F("fHistEPCorrChiSigma", "EP_V0 #chi", "EP_V0 #sigma", 50, -TMath::Pi()/2., TMath::Pi()/2., 50, -TMath::Pi()/2., TMath::Pi()/2., i);
559  fHistIntegralCorrelations[i] = BookTH2F("fHistIntegralCorrelations", "square [GeV/c/A]", "circle [GeVc/A]", 100, 0, 100, 100, 0, 100);
560  fProfIntegralCorrelations[i] = new TProfile(Form("fProfIntegralCorrelations_%i", i), Form("fProfIntegralCorrelations_%i", i), 100, 0, 100);
561  fProfIntegralCorrelations[i]->GetXaxis()->SetTitle("RC energy, #eta #varphi scale");
562  fProfIntegralCorrelations[i]->GetYaxis()->SetTitle("#phi / #eta, #varphi");
564  fHistPsiVZEROALeadingJet[i] = BookTH3F("fHistPsiVZEROALeadingJet", "p_{t} [GeV/c]", "#Psi_{VZEROA}", "#varphi_{jet}", 70, 0, 210, 50, -1.*TMath::Pi()/2., TMath::Pi()/2., 50, phiMin, phiMax, i);
565  fHistPsiVZEROCLeadingJet[i] = BookTH3F("fHistPsiVZEROCLeadingJet", "p_{t} [GeV/c]", "#Psi_{VZEROC}", "#varphi_{jet}", 70, 0, 210, 50, -1.*TMath::Pi()/2., TMath::Pi()/2., 50, phiMin, phiMax, i);
566  fHistPsiVZEROCombLeadingJet[i] = BookTH3F("fHistPsiVZEROCombLeadingJet", "p_{t} [GeV/c]", "#Psi_{VZEROComb}", "#varphi_{jet}", 70, 0, 210, 50, -1.*TMath::Pi()/2., TMath::Pi()/2., 50, phiMin, phiMax, i);
567  fHistPsi2Correlation[i] = BookTH3F("fHistPsi2Correlation", "#Psi_{TPC}", "#Psi_{VZEROA}", "#Psi_{VZEROC}", 20, -1.*TMath::Pi()/2., TMath::Pi()/2., 20, -1.*TMath::Pi()/2., TMath::Pi()/2., 20, -1.*TMath::Pi()/2., TMath::Pi()/2., i);
568  fHistLeadingJetBackground[i] = BookTH2F("fHistLeadingJetBackground", "#Delta #eta (leading jet with, without sub)", "Delta #varphi (leading jet with, without sub)", 50, 0., 2, 50, 0., TMath::TwoPi(), i);
569  // trigger qa
570  fHistTriggerQAIn[i] = BookTH2F("fHistTriggerQAIn", "trigger configuration", "p_{T}^{jet} (GeV/c) in-plane jets", 16, 0.5, 16.5, 70, -100, 250, i);
571  fHistTriggerQAOut[i] = BookTH2F("fHistTriggerQAOut", "trigger configuration", "p_{T}^{jet} (GeV/c) out-of-plane jets", 16, 0.5, 16.5, 70, -100, 250, i);
572  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(1, "no trigger");
573  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(2, "kAny");
574  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(3, "kAnyINT");
575  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(4, "kMB");
576  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(5, "kCentral");
577  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(6, "kSemiCentral");
578  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(7, "kEMCEJE");
579  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(8, "kEMCEGA");
580  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(9, "kEMCEJE & kMB");
581  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(10, "kEMCEJE & kCentral");
582  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(11, "kEMCEJE & kSemiCentral");
583  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(12, "kEMCEJE & all min bias");
584  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(13, "kEMCEGA & kMB");
585  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(14, "kEMCEGA & kCentral");
586  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(15, "kEMCEGA & kSemiCentral");
587  fHistTriggerQAIn[i]->GetXaxis()->SetBinLabel(16, "kEMCEGA & all min bias");
588  fHistTriggerQAIn[i]->LabelsOption("v");
589  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(1, "no trigger");
590  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(2, "kAny");
591  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(3, "kAnyINT");
592  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(4, "kMB");
593  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(5, "kCentral");
594  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(6, "kSemiCentral");
595  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(7, "kEMCEJE");
596  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(8, "kEMCEGA");
597  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(9, "kEMCEJE & kMB");
598  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(10, "kEMCEJE & kCentral");
599  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(11, "kEMCEJE & kSemiCentral");
600  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(12, "kEMCEJE & all min bias");
601  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(13, "kEMCEGA & kMB");
602  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(14, "kEMCEGA & kCentral");
603  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(15, "kEMCEGA & kSemiCentral");
604  fHistTriggerQAOut[i]->GetXaxis()->SetBinLabel(16, "kEMCEGA & all min bias");
605  fHistTriggerQAOut[i]->LabelsOption("v");
606  }
607  }
608 
609  if(fFillQAHistograms) {
610  Int_t low(fCentralityClasses->At(0)), up(fCentralityClasses->At(fCentralityClasses->GetSize()-1));
611  Int_t diff(TMath::Abs(up-low));
612  // event plane estimates and quality
613  fHistPsiVZEROAV0M = BookTH2F("fHistPsiVZEROAV0M", "V0M", "#Psi_{2, VZEROA}", diff, low, up, 40, -.5*TMath::Pi(), .5*TMath::Pi());
614  fHistPsiVZEROCV0M = BookTH2F("fHistPsiVZEROCV0M", "V0M", "#Psi_{2, VZEROC}", diff, low, up, 40, -.5*TMath::Pi(), .5*TMath::Pi());
615  fHistPsiVZEROVV0M = BookTH2F("fHistPsiVZEROV0M", "V0M", "#Psi_{2, VZERO}", diff, low, up, 40, -.5*TMath::Pi(), .5*TMath::Pi());
616  fHistPsiTPCV0M = BookTH2F("fHistPsiTPCV0M", "V0M", "#Psi_{2, TRK}", diff, low, up, 40, -.5*TMath::Pi(), .5*TMath::Pi());
617  fHistPsiVZEROATRK = BookTH2F("fHistPsiVZEROATRK", "TRK", "#Psi_{2, VZEROA}", diff, low, up, 40, -.5*TMath::Pi(), .5*TMath::Pi());
618  fHistPsiVZEROCTRK = BookTH2F("fHistPsiVZEROCTRK", "TRK", "#Psi_{2, VZEROC}", diff, low, up, 40, -.5*TMath::Pi(), .5*TMath::Pi());
619  fHistPsiVZEROTRK = BookTH2F("fHistPsiVZEROTRK", "TRK", "#Psi_{2, VZERO}", diff, low, up, 40, -.5*TMath::Pi(), .5*TMath::Pi());
620  fHistPsiTPCTRK = BookTH2F("fHistPsiTPCTRK", "TRK", "#Psi_{2, TRK}", diff, low, up, 40, -.5*TMath::Pi(), .5*TMath::Pi());
621  }
622  // background
623  for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i ++) {
624  fHistRhoPackage[i] = BookTH1F("fHistRhoPackage", "#rho [GeV/c]", 100, 0, 150, i);
625  fHistRho[i] = BookTH1F("fHistRho", "#rho [GeV/c]", 100, 0, 150, i);
626  }
627  fHistRhoVsMult = BookTH2F("fHistRhoVsMult", "multiplicity", "#rho [GeV/c]", 100, 0, 4000, 100, 0, 250);
628  fHistRhoVsCent = BookTH2F("fHistRhoVsCent", "centrality", "#rho [GeV/c]", 100, 0, 100, 100, 0, 250);
629  fHistRhoAVsMult = BookTH2F("fHistRhoAVsMult", "multiplicity", "#rho * A (jet) [GeV/c]", 100, 0, 4000, 100, 0, 50);
630  fHistRhoAVsCent = BookTH2F("fHistRhoAVsCent", "centrality", "#rho * A (jet) [GeV/c]", 100, 0, 100, 100, 0, 50);
631 
632  TString detector("");
633  switch (fDetectorType) {
634  case kTPC : detector+="TPC";
635  break;
636  case kVZEROA : detector+="VZEROA";
637  break;
638  case kVZEROC : detector+="VZEROC";
639  break;
640  case kVZEROComb : detector+="VZEROComb";
641  break;
642  case kFixedEP : detector+="FixedEP";
643  break;
644  default: break;
645  }
646  // delta pt distributions
647  for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i ++) {
648  if(fFillQAHistograms) fHistRCPhiEta[i] = BookTH2F("fHistRCPhiEta", "#phi (RC)", "#eta (RC)", 40, phiMin, phiMax, 40, etaMin, etaMax, i);
649  fHistRhoVsRCPt[i] = BookTH2F("fHistRhoVsRCPt", "p_{t} (RC) [GeV/c]", "#rho * A (RC) [GeV/c]", 100, 0, 300, 100, 0, 350, i);
650  fHistRCPt[i] = BookTH1F("fHistRCPt", "p_{t} (RC) [GeV/c]", 130, -20, 150, i);
651  if(fFillQAHistograms) fHistRCPhiEtaExLJ[i] = BookTH2F("fHistRCPhiEtaExLJ", "#phi (RC)", "#eta (RC)", 40, phiMin, phiMax, 40, etaMin, etaMax, i);
652  fHistDeltaPtDeltaPhi2[i] = BookTH2F("fHistDeltaPtDeltaPhi2", Form("#phi - #Psi_{2, %s}", detector.Data()), "#delta p_{t} [GeV/c]", 40, 0, TMath::Pi(), 400, -70, 130, i);
653  fHistDeltaPtDeltaPhi2Rho0[i] = BookTH2F("fHistDeltaPtDeltaPhi2Rho0", Form("#phi - #Psi_{2, %s}", detector.Data()), "#delta p_{t} [GeV/c]", 40, 0, TMath::Pi(), 400, -70, 130, i);
654  fHistRhoVsRCPtExLJ[i] = BookTH2F("fHistRhoVsRCPtExLJ", "p_{t} (RC) [GeV/c]", "#rho * A (RC) [GeV/c]", 100, 0, 300, 100, 0, 350, i);
655  fHistRCPtExLJ[i] = BookTH1F("fHistRCPtExLJ", "p_{t} (RC) [GeV/c]", 130, -20, 150, i);
656  fHistDeltaPtDeltaPhi2ExLJ[i] = BookTH2F("fHistDeltaPtDeltaPhi2ExLJ", Form("#phi - #Psi_{2, %s}", detector.Data()), "#delta p_{t} [GeV/c]", 40, 0, TMath::Pi(), 400, -70, 130, i);
657  fHistDeltaPtDeltaPhi2ExLJRho0[i] = BookTH2F("fHistDeltaPtDeltaPhi2ExLJRho0", Form("#phi - #Psi_{2, %s}", detector.Data()), "#delta p_{t} [GeV/c]", 40, 0, TMath::Pi(), 400, -70, 130, i);
658  // jet histograms (after kinematic cuts)
659  fHistJetPtRaw[i] = BookTH1F("fHistJetPtRaw", "p_{t, jet} RAW [GeV/c]", 200, -50, 150, i);
660  fHistJetPt[i] = BookTH1F("fHistJetPt", "p_{t, jet} [GeV/c]", 350, -100, 250, i);
661  if(fFillQAHistograms) fHistJetEtaPhi[i] = BookTH2F("fHistJetEtaPhi", "#eta", "#phi", 100, etaMin, etaMax, 100, phiMin, phiMax, i);
662  fHistJetPtArea[i] = BookTH2F("fHistJetPtArea", "p_{t, jet} [GeV/c]", "Area", 175, -100, 250, 30, 0, 0.9, i);
663  fHistJetPtEta[i] = BookTH2F("fHistJetPtEta", "p_{t, jet} [GeV/c]", "Eta", 175, -100, 250, 30, etaMin, etaMax, i);
664  fHistJetPtConstituents[i] = BookTH2F("fHistJetPtConstituents", "p_{t, jet} [GeV/c]", "no. of constituents", 350, -100, 250, 60, 0, 150, i);
665  fHistJetEtaRho[i] = BookTH2F("fHistJetEtaRho", "#eta", "#rho", 100, etaMin, etaMax, 100, 0, 300, i);
666  // in plane and out of plane spectra
667  fHistJetPsi2Pt[i] = BookTH2F("fHistJetPsi2Pt", Form("#phi_{jet} - #Psi_{2, %s}", detector.Data()), "p_{t, jet} [GeV/c]", 40, 0., TMath::Pi(), 350, -100, 250, i);
668  fHistJetLJPsi2Pt[i] = BookTH3F("fHistJetLJPsi2Pt", Form("#phi_{part} - #Psi_{2, %s}", detector.Data()), "p_{t, jet} [GeV/c]", "p_{t, leading track}", 40, 0., TMath::Pi(), 350, -100, 250, 200, 0, 50, i);
669  fHistJetLJPsi2PtRatio[i] = BookTH3F("fHistJetLJPsi2PtRatio", Form("#phi_{part} - #Psi_{2, %s}", detector.Data()), Form("#phi_{jet} - #Psi_{2, %s}", detector.Data()), "p_{t, jet} [GeV/c]", 40, 0., TMath::Pi(), 40, 0., TMath::Pi(), 350, -100, 250, i);
670 
671  fHistJetPsi2PtRho0[i] = BookTH2F("fHistJetPsi2PtRho0", Form("#phi_{jet} - #Psi_{2, %s}", detector.Data()), "p_{t, jet} [GeV/c]", 40, 0., TMath::Pi(), 350, -100, 250, i);
672  // profiles for all correlator permutations which are necessary to calculate each second and third order event plane resolution
673  fProfV2Resolution[i] = new TProfile(Form("fProfV2Resolution_%i", i), Form("fProfV2Resolution_%i", i), 11, -0.5, 10.5);
674  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{VZEROA} - #Psi_{VZEROC}))>");
675  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(4, "<cos(2(#Psi_{VZEROC} - #Psi_{VZEROA}))>");
676  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(5, "<cos(2(#Psi_{VZEROA} - #Psi_{TPC}))>");
677  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(6, "<cos(2(#Psi_{TPC} - #Psi_{VZEROA}))>");
678  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(7, "<cos(2(#Psi_{VZEROC} - #Psi_{TPC}))>");
679  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(8, "<cos(2(#Psi_{TPC} - #Psi_{VZEROC}))>");
680  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(9, "<cos(2(#Psi_{VZERO} - #Psi_{TPC_A}))>");
681  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(10, "<cos(2(#Psi_{VZERO} - #Psi_{TPC_B}))>");
682  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(11, "<cos(2(#Psi_{TPC_A} - #Psi_{TPC_B}))>");
683  fOutputList->Add(fProfV2Resolution[i]);
684  fProfV3Resolution[i] = new TProfile(Form("fProfV3Resolution_%i", i), Form("fProfV3Resolution_%i", i), 11, -0.5, 10.5);
685  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(3, "<cos(3(#Psi_{VZEROA} - #Psi_{VZEROC}))>");
686  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(4, "<cos(3(#Psi_{VZEROC} - #Psi_{VZEROA}))>");
687  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(5, "<cos(3(#Psi_{VZEROA} - #Psi_{TPC}))>");
688  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(6, "<cos(3(#Psi_{TPC} - #Psi_{VZEROA}))>");
689  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(7, "<cos(3(#Psi_{VZEROC} - #Psi_{TPC}))>");
690  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(8, "<cos(3(#Psi_{TPC} - #Psi_{VZEROC}))>");
691  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(9, "<cos(3(#Psi_{VZERO} - #Psi_{TPC_A}))>");
692  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(10, "<cos(3(#Psi_{VZERO} - #Psi_{TPC_B}))>");
693  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(11, "<cos(3(#Psi_{TPC_A} - #Psi_{TPC_B}))>");
694  fOutputList->Add(fProfV3Resolution[i]);
695  }
696  // vn profile
697  Float_t temp[fCentralityClasses->GetSize()];
698  for(Int_t i(0); i < fCentralityClasses->GetSize(); i++) temp[i] = fCentralityClasses->At(i);
699  fProfV2 = new TProfile("fProfV2", "fProfV2", fCentralityClasses->GetSize()-1, temp);
700  fProfV3 = new TProfile("fProfV3", "fProfV3", fCentralityClasses->GetSize()-1, temp);
701  fOutputList->Add(fProfV2);
702  fOutputList->Add(fProfV3);
703  switch (fFitModulationType) {
704  case kQC2 : {
705  fProfV2Cumulant = new TProfile("fProfV2Cumulant", "fProfV2Cumulant", fCentralityClasses->GetSize()-1, temp);
706  fProfV3Cumulant = new TProfile("fProfV3Cumulant", "fProfV3Cumulant", fCentralityClasses->GetSize()-1, temp);
708  fOutputList->Add(fProfV3Cumulant);
709  } break;
710  case kQC4 : {
711  fProfV2Cumulant = new TProfile("fProfV2Cumulant", "fProfV2Cumulant", fCentralityClasses->GetSize()-1, temp);
712  fProfV3Cumulant = new TProfile("fProfV3Cumulant", "fProfV3Cumulant", fCentralityClasses->GetSize()-1, temp);
714  fOutputList->Add(fProfV3Cumulant);
715  } break;
716  default : break;
717  }
718  // for the histograms initialized below, binning is fixed to runnumbers or flags
721  if(fFillQAHistograms) {
722  fHistRunnumbersEta = new TH2F("fHistRunnumbersEta", "fHistRunnumbersEta", fExpectedRuns->GetSize()+1, -.5, fExpectedRuns->GetSize()+.5, 100, -1.1, 1.1);
723  fHistRunnumbersEta->Sumw2();
725  fHistRunnumbersPhi = new TH2F("fHistRunnumbersPhi", "fHistRunnumbersPhi", fExpectedRuns->GetSize()+1, -.5, fExpectedRuns->GetSize()+.5, 100, -0.2, TMath::TwoPi()+0.2);
726  fHistRunnumbersPhi->Sumw2();
728  for(Int_t i(0); i < fExpectedRuns->GetSize(); i++) {
729  fHistRunnumbersPhi->GetXaxis()->SetBinLabel(i+1, Form("%i", fExpectedRuns->At(i)));
730  fHistRunnumbersEta->GetXaxis()->SetBinLabel(i+1, Form("%i", fExpectedRuns->At(i)));
731  }
732  fHistRunnumbersPhi->GetXaxis()->SetBinLabel(fExpectedRuns->GetSize()+1, "undetermined");
733  fHistRunnumbersEta->GetXaxis()->SetBinLabel(fExpectedRuns->GetSize()+1, "undetermined");
734  if(fCollisionType == kPbPb10h) {
735  // control histo to see if the calibration was properly kickstarted
736  fHistRunnumbersCaliInfo = new TH1I("fHistRunnumbersCaliInfo", "fHistRunnumbersCaliInfo", fExpectedRuns->GetSize()+1, -.5, fExpectedRuns->GetSize()+.5);
738  for(Int_t i(0); i < fExpectedRuns->GetSize(); i++) {
739  fHistRunnumbersCaliInfo->GetXaxis()->SetBinLabel(i+1, Form("%i", fExpectedRuns->At(i)));
740  }
741  fHistRunnumbersCaliInfo->GetXaxis()->SetBinLabel(fExpectedRuns->GetSize()+1, "undetermined");
742  }
743  }
744  fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 54, -0.5, 54.5);
745  fHistSwap = new TH1F("fHistSwap", "fHistSwap", 20, 0, TMath::TwoPi());
746  if(fUsePtWeight) fHistSwap->Sumw2();
747 
752  for(Int_t i(0); i < 10; i++) {
753  if(fEventPlaneWeights[i]) {
754  // add the original event plane distribution histogram
755  fOutputList->Add((TH1F*)(fEventPlaneWeights[i]->Clone(Form("EP_distribution_original_cen_%i", i))));
756  // calculate the weights that will actually be used and store them
759  }
760  }
761  // increase readability of output list
762  fOutputList->Sort();
763  // cdf and pdf of chisquare distribution
764  fHistPvalueCDF = BookTH1F("fHistPvalueCDF", "CDF #chi^{2}", 50, 0, 1);
765  fHistPvalueCDFCent = BookTH2F("fHistPvalueCDFCent", "centrality", "p-value", 40, 0, 100, 40, 0, 1);
766  fHistChi2Cent = BookTH2F("fHistChi2Cent", "centrality", "#tilde{#chi^{2}}", 100, 0, 100, 100, 0, 5);
767  fHistPChi2 = BookTH2F("fHistPChi2", "p-value", "#tilde{#chi^{2}}", 1000, 0, 1, 100, 0, 5);
768  fHistKolmogorovTest = BookTH1F("fHistKolmogorovTest", "KolmogorovTest", 50, 0, 1);
769  fHistKolmogorovTestCent = BookTH2F("fHistKolmogorovTestCent", "centrality", "Kolmogorov p", 40, 0, 100, 45, 0, 1);
770  fHistPvalueCDFROOT = BookTH1F("fHistPvalueCDFROOT", "CDF #chi^{2} ROOT", 50, 0, 1);
771  fHistPvalueCDFROOTCent = BookTH2F("fHistPvalueCDFROOTCent", "centrality", "p-value ROOT", 40, 0, 100, 45, 0, 1);
772  fHistChi2ROOTCent = BookTH2F("fHistChi2ROOTCent", "centrality", "#tilde{#chi^{2}}", 40, 0, 100, 45, 0, 5);
773  fHistPChi2Root = BookTH2F("fHistPChi2Root", "p-value", "#tilde{#chi^{2}} ROOT", 1000, 0, 1, 100, 0, 5);
774  fHistPKolmogorov = BookTH2F("fHistPKolmogorov", "p-value", "kolmogorov p",40, 0, 1, 40, 0, 1);
775  fHistRhoStatusCent = BookTH2F("fHistRhoStatusCent", "centrality", "status [-1=lin was better, 0=ok, 1 = failed]", 101, -1, 100, 3, -1.5, 1.5);
776  fHistUndeterminedRunQA = BookTH1F("fHistUndeterminedRunQA", "runnumber", 10, 0, 10);
777 
778  // Mar 24 2016 - add some figures that are missing for the thesis
779  for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i ++) {
780  fHistRhoEtaBC[i] = BookTH2F("fHistRhoEtaBC", "#rho [GeV/c]", "#eta", 100, 0, 150, 50, -1, 1, i);
781  fHistJetPtBC[i] = BookTH1F("fHistJetPtBC", "p_{t, jet} [GeV/c]", 350, -100, 250, i);
782  fHistJetEtaPhiBC[i] = BookTH2F("fHistJetEtaPhiBC", "#eta", "#phi", 100, etaMin, etaMax, 100, phiMin, phiMax, i);
783  fHistJetPtAreaBC[i] = BookTH2F("fHistJetPtAreaBC", "p_{t, jet} [GeV/c]", "Area", 175, -100, 250, 30, 0, 0.3, i);
784  }
785  fHistQxV0aBC = BookTH2F("fHistQxV0aBC", "Q_{x} V0A", "centrality class", 140, -700, 700, 10, -.5, 9.5);
786  fHistQyV0aBC = BookTH2F("fHistQyV0aBC", "Q_{y} V0A", "centrality class", 140, -700, 700, 10, -.5, 9.5);
787  fHistQxV0cBC = BookTH2F("fHistQxV0cBC", "Q_{x} V0C", "centrality class", 140, -700, 700, 10, -.5, 9.5);
788  fHistQyV0cBC = BookTH2F("fHistQyV0cBC", "Q_{y} V0C", "centrality class", 140, -700, 700, 10, -.5, 9.5);
789  fHistQxV0a = BookTH2F("fHistQxV0a", "Q_{x} V0A", "centrality class", 100, -10, 10, 10, -.5, 9.5);
790  fHistQyV0a = BookTH2F("fHistQyV0a", "Q_{y} V0A", "centrality class", 100, -10, 10, 10, -.5, 9.5);
791  fHistQxV0c = BookTH2F("fHistQxV0c", "Q_{x} V0C", "centrality class", 100, -10, 10, 10, -.5, 9.5);
792  fHistQyV0c = BookTH2F("fHistQyV0c", "Q_{y} V0C", "centrality class", 100, -10, 10, 10, -.5, 9.5);
793  fHistMultVsCellBC = BookTH2F("fHistMultVsCellBC", "channel", "multiplicty", 64, -.5, 63.5, 100, 0, 1000);
794  fHistMultVsCell = BookTH2F("fHistMultVsCell", "channel", "multiplicty", 64, -.5, 63.5, 100, 0, 1000);
795  fHistEPBC = BookTH1F("fHistEPBC", "#Psi_{EP, 2}, uncalibrated", 100, -0.5*TMath::Pi(), 0.5*TMath::Pi());
796  fHistEP = BookTH1F("fHistEP", "#Psi_{EP, 2}, calibrated", 100, -0.5*TMath::Pi(), 0.5*TMath::Pi());
797 
798 
799 
800 
801  PostData((fCreateHisto) ? 2 : 1, fOutputList);
802 
803  switch (fRunModeType) {
804  case kLocal : {
805  fOutputListGood = new TList();
806  fOutputListGood->SetOwner(kTRUE);
807  fOutputListBad = new TList();
808  fOutputListBad->SetOwner(kTRUE);
809  PostData((fCreateHisto) ? 3 : 2, fOutputListGood);
810  PostData((fCreateHisto) ? 4 : 3, fOutputListBad);
811  } break;
812  default: break;
813  }
814 
815  // get the containers
816  fTracksCont = GetParticleContainer("Tracks");
817  fClusterCont = GetClusterContainer(0); // get the default cluster container
818  fJetsCont = GetJetContainer("Jets");
819 }
820 //_____________________________________________________________________________
822 {
823  // called for each accepted event (call made from user exec of parent class)
824  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
825  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
826  #endif
827  if(!fTracks||!fJets||!fRho) {
828  if(!fTracks) printf(" > Failed to retrieve fTracks ! < \n");
829  if(!fJets) printf(" > Failed to retrieve fJets ! < \n");
830  if(!fRho) printf(" > Failed to retrieve fRho ! < \n");
831  return kFALSE;
832  }
834  // reject the event if expected data is missing
835  if(!PassesCuts(InputEvent())) return kFALSE;
836  // cache the leading jet within acceptance
838  // set the rho value
839  fLocalRho->SetVal(fRho->GetVal());
840  // place holder arrays for the event planes
841  //
842  // [0][0] psi2a [1,0] psi2c
843  // [0][1] psi3a [1,1] psi3c
844  Double_t vzero[2][2];
845  /* for the combined vzero event plane
846  * [0] psi2 [1] psi3
847  * not fully implmemented yet, use with caution ! */
848  Double_t vzeroComb[2];
849  // [0] psi2 [1] psi3
850  Double_t tpc[2];
851  // evaluate the actual event planes
852  switch (fDetectorType) {
853  case kFixedEP : {
854  // for fixed, fix all ep's to default values
855  tpc[0] = 0.; tpc[1] = 1.;
856  vzero[0][0] = 0.; vzero[0][1] = 1.;
857  vzero[1][0] = 0.; vzero[1][1] = 1.;
858  vzeroComb[0] = 0.; vzeroComb[1] = 1.;
859  } break;
860  default : {
861  // else grab the actual data
865  } break;
866  }
867  Double_t psi2(-1), psi3(-1);
868  // arrays which will hold the fit parameters
869  switch (fDetectorType) { // determine the detector type for the rho fit
870  case kTPC : { psi2 = tpc[0]; psi3 = tpc[1]; } break;
871  case kVZEROA : { psi2 = vzero[0][0]; psi3 = vzero[0][1]; } break;
872  case kVZEROC : { psi2 = vzero[1][0]; psi3 = vzero[1][1]; } break;
873  case kVZEROComb : { psi2 = vzeroComb[0]; psi3 = vzeroComb[1];} break;
874  case kFixedEP : { psi2 = 0.; psi3 = 1.;} break;
875  default : break;
876  }
877  // if requested extract the event plane weight
878  fEventPlaneWeight = 1.; // ALWAYS reset to 1 here to avoid recycling an old weight if the next if-statement fails
880  // get the weight from the corresponding
881  fEventPlaneWeight = fEventPlaneWeights[fInCentralitySelection]->GetBinContent(fEventPlaneWeights[fInCentralitySelection]->FindBin(psi2));
882  }
883  // if requested store the acceptance weights
884  if(fAcceptanceWeights) {
885  Double_t percIn(0.), percOut(0.), percLost(0.);
886  NumericalOverlap(GetJetContainer()->GetJetEtaMin(), GetJetContainer()->GetJetEtaMax(),
887  psi2, percIn, percOut, percLost);
888  fHistCentralityPercIn->Fill(fCent, percIn);
889  fHistCentralityPercOut->Fill(fCent, percOut);
890  fHistCentralityPercLost->Fill(fCent, percLost);
891  }
892  switch (fFitModulationType) { // do the fits
893  case kNoFit : {
894  switch (fCollisionType) {
895  case kPythia : { // background is zero for pp jets
896  fFitModulation->FixParameter(0, 0);
897  fLocalRho->SetVal(0);
898  } break;
899  default : {
900  fFitModulation->FixParameter(0, fLocalRho->GetVal());
901  } break;
902  }
903  } break;
904  case kV2 : { // only v2
905  if(CorrectRho(psi2, psi3)) {
906  fProfV2->Fill(fCent, fFitModulation->GetParameter(3));
907  if(fUserSuppliedR2) {
908  Double_t r(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
909  if(r > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r);
910  }
911  CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
912  }
913  } break;
914  case kV3 : { // only v3
915  if(CorrectRho(psi2, psi3)) {
916  if(fUserSuppliedR3) {
917  Double_t r(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
918  if(r > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r);
919  }
920  fProfV3->Fill(fCent, fFitModulation->GetParameter(3));
921  CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
922  }
923  } break;
924  case kQC2 : { // qc2 analysis
925  if(CorrectRho(psi2, psi3)) {
927  // note for the qc method, resolution is REVERSED to go back to v2obs
928  Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
929  Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
930  if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
931  if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)*r3);
932  }
933  if (fUsePtWeight) { // use weighted weights
934  Double_t dQCnM11 = (fNoEventWeightsForQC) ? 1. : QCnM11();
935  fProfV2->Fill(fCent, fFitModulation->GetParameter(3), dQCnM11);
936  fProfV3->Fill(fCent, fFitModulation->GetParameter(7), dQCnM11);
937  } else {
938  Double_t dQCnM = (fNoEventWeightsForQC) ? 2. : QCnM();
939  fProfV2->Fill(fCent, fFitModulation->GetParameter(3), dQCnM*(dQCnM-1));
940  fProfV3->Fill(fCent, fFitModulation->GetParameter(7), dQCnM*(dQCnM-1));
941  }
942  CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
943  }
944  } break;
945  case kQC4 : {
946  if(CorrectRho(psi2, psi3)) {
948  // note for the qc method, resolution is REVERSED to go back to v2obs
949  Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
950  Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
951  if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
952  if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)*r3);
953  }
954  if (fUsePtWeight) { // use weighted weights
955  fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5)/*, QCnM1111()*/);
956  fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5)/*, QCnM1111()*/);
957  } else {
958  fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
959  fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
960  }
961  }
962  CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
963  } break;
964  default : {
965  if(CorrectRho(psi2, psi3)) {
967  Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
968  Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
969  if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r2);
970  if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)/r3);
971  }
972  fProfV2->Fill(fCent, fFitModulation->GetParameter(3));
973  fProfV3->Fill(fCent, fFitModulation->GetParameter(7));
974  CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
975  }
976  } break;
977  }
978  // if all went well, update the local rho parameter
980  // and only at this point can the leading jet after rho subtraction be evaluated
982  // fill a number of histograms. event qa needs to be filled first as it also determines the runnumber for the track qa
984  if(fFillHistograms) FillHistogramsAfterSubtraction(psi2, vzero, vzeroComb, tpc);
985  // send the output to the connected output container
986  PostData((fCreateHisto) ? 2 : 1, fOutputList);
987  switch (fRunModeType) {
988  case kLocal : {
989  PostData((fCreateHisto) ? 3 : 2, fOutputListGood);
990  PostData((fCreateHisto) ? 4 : 3, fOutputListBad);
991  } break;
992  default: break;
993  }
994  return kTRUE;
995 }
996 //_____________________________________________________________________________
998 {
999  // for stand alone, avoid framework event setup
1000  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1001  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1002  #endif
1003  switch (fCollisionType) {
1004  case kJetFlowMC : {
1005  // need to call ExecOnce as it is not loaded otherwise
1008  } break;
1009 // case kPbPb10h : {
1010 // // bypass framework event selection. additional check for fTracks
1011 // // to avoid the situation where base classes are never initialized
1012 // if(fTracks && fTracks->GetEntriesFast() > 0) AliAnalysisTaskJetV2::Run();
1013 // else AliAnalysisTaskSE::Exec(c);
1014 // } break;
1015  default : {
1016  AliAnalysisTaskSE::Exec(c);
1017  } break;
1018  }
1019 }
1020 //_____________________________________________________________________________
1022 {
1023  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_2
1024  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1025  #endif
1026  // numerically integrate with finite resolution
1027  // idea is the following:
1028  // 1) choose a vector phi
1029  // 2) see if it is in a region of overlap between detector and in/out of plane spectrum
1030  // 3) bookkeep percentages over overlap
1031  Double_t a(psi2 - TMath::Pi()/4.);
1032  // poor man's appproach: fix the frame
1033  if(a < 0) a += TMath::Pi();
1034  // set the rest of the event
1035  Double_t b(a + TMath::Pi()/2.);
1036  Double_t c(b + TMath::Pi()/2.);
1037  Double_t d(c + TMath::Pi()/2.);
1038  Double_t e(d + TMath::Pi()/2.); // may seem mysterious but here for good reasons
1039  // get percetnages
1040  Double_t interval(TMath::TwoPi() / 1000.);
1041  percIn = 0.;
1042  percOut = 0.;
1043  percLost = 0.;
1044  Int_t status(-1);
1045  // automagically do the integration
1046  for(Double_t i = a; i < a+TMath::TwoPi()-interval; i += interval) {
1047  status = OverlapsWithPlane(x1, x2, a, b, c, d, e, i);
1048  if(status == 0 ) percLost += .001;
1049  else if(status == 1 ) percIn += 0.001;
1050  else if(status == 2 ) percOut += 0.001;
1051  }
1052 }
1053 //_____________________________________________________________________________
1055  Double_t x1, Double_t x2, // detector geometry relative to ep
1056  Double_t a, Double_t b, Double_t c, Double_t d, Double_t e, // in-plane, out-of-plane boundaries (see comments)
1057  Double_t phi) // variable
1058 {
1059  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_2
1060  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1061  #endif
1062  // 'numerical integration' of geometric overlap
1063  //
1064  // works as follows: for a given vector phi determines whether
1065  // or not this vector points towards an overlap region of
1066  // detector geometry and plane (in or out)
1067  //
1068  // returns
1069  // 1) if overlap with in plane
1070  // 2) if overlap with out of plane
1071  // 0) if no overlap at all
1072  Int_t overlap(0);
1073  // check for condition in-plane
1074  // conditions are always checked as
1075  // 1) is the angle within in-plane sector?
1076  // 2) is the angle also within detector acceptance?
1077  if(phi > a && phi < b && phi > x1 && phi < x2) overlap = 1;
1078  if(phi > c && phi < d && phi > x1 && phi < x2) overlap = 1;
1079  // likewise for out-of-plane
1080  if(phi > b && phi < c && phi > x1 && phi < x2) overlap = 2;
1081  if(phi > d && phi < e && phi > x1 && phi < x2) overlap = 2;
1082 
1083  // life would be so much easier if the detector was flat instead of cylindrical ....
1084  x1+=TMath::TwoPi();
1085  x2+=TMath::TwoPi();
1086 
1087  if(phi > a && phi < b && phi > x1 && phi < x2) overlap = 1;
1088  if(phi > c && phi < d && phi > x1 && phi < x2) overlap = 1;
1089  // likewise for out-of-plane
1090  if(phi > b && phi < c && phi > x1 && phi < x2) overlap = 2;
1091  if(phi > d && phi < e && phi > x1 && phi < x2) overlap = 2;
1092 
1093  return overlap;
1094 }
1095 //_____________________________________________________________________________
1097 {
1098  // return chi for given resolution to combine event plane estimates from two subevents
1099  // see Phys. Rev. C no. CS6346 (http://arxiv.org/abs/nucl-ex/9805001)
1100  Double_t chi(2.), delta(1.), con((TMath::Sqrt(TMath::Pi()))/(2.*TMath::Sqrt(2)));
1101  for (Int_t i(0); i < 15; i++) {
1102  chi = ((con*chi*TMath::Exp(-chi*chi/4.)*(TMath::BesselI0(chi*chi/4.)+TMath::BesselI1(chi*chi/4.))) < res) ? chi + delta : chi - delta;
1103  delta = delta / 2.;
1104  }
1105  return chi;
1106 }
1107 //_____________________________________________________________________________
1109 {
1110  // get the vzero event plane (a and c separately)
1111  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1112  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1113  #endif
1114  switch (fCollisionType) {
1115  case kPbPb10h : {
1116  // for 10h data, get the calibrated q-vector from the database
1117  Double_t QA2[] = {-999., -999.};
1118  Double_t QA3[] = {-999., -999.};
1119  Double_t QC2[] = {-999., -999.};
1120  Double_t QC3[] = {-999., -999.};
1121  CalculateQvectorVZERO(QA2, QC2, QA3, QC3);
1122  vzero[0][0] = .5*TMath::ATan2(QA2[1], QA2[0]);
1123  vzero[1][0] = .5*TMath::ATan2(QC2[1], QC2[0]);
1124  vzero[0][1] = (1./3.)*TMath::ATan2(QA3[1], QA3[0]);
1125  vzero[1][1] = (1./3.)*TMath::ATan2(QC3[1], QC3[0]);
1126  return; // paranoid return
1127  } break;
1128  default: {
1129  // by default use the ep from the event header (make sure EP selection task is enabeled!)
1130  Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0), h(0);
1131  vzero[0][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, a, b);
1132  vzero[1][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, c, d);
1133  vzero[0][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, e, f);
1134  vzero[1][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, g, h);
1135  return;
1136  }
1137  }
1138 }
1139 //_____________________________________________________________________________
1141 {
1142  // return the combined vzero event plane
1143  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1144  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1145  #endif
1146 
1147  // define some placeholders
1148  Double_t Q2[] = {-999., -999.};
1149  Double_t Q3[] = {-999., -999.};
1150 
1151  switch (fCollisionType) {
1152  // for 10h data call calibration info
1153  case kPbPb10h : {
1154  // get the calibrated q-vectors
1156  comb[0] = .5*TMath::ATan2(Q2[1], Q2[0]);
1157  comb[1] = (1./3.)*TMath::ATan2(Q3[1], Q3[0]);
1158  } break;
1159  default : {
1160  // for all other types use calibrated event plane from the event header
1161  //
1162  // note that the code is a bit messy here. for 10h data retrieving q-vectors of
1163  // the separate vzero detectors and combining the q-vectors have dedicated functions.
1164  // for 11h however this is all done in this function (the lines below)
1165  // reason is that the procedure is much shorter as the calibration is done in another task
1166  //
1167  // define some pleaceholders to the values by reference
1168  Double_t qx2a(0), qy2a(0), qx2c(0), qy2c(0), qx3a(0), qy3a(0), qx3c(0), qy3c(0);
1169  // get the q-vectors by reference
1170  InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, qx2a, qy2a);
1171  InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, qx2c, qy2c);
1172  InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, qx3a, qy3a);
1173  InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, qx3c, qy3c);
1174 
1175  // get cache index and retrieve the chi weights for this centrality
1176  Int_t VZEROcentralityBin(GetVZEROCentralityBin());
1177  Double_t chi2A(1);
1178  Double_t chi2C(1);
1179  Double_t chi3A(1);
1180  Double_t chi3C(1);
1181 
1182  switch (fWeightForVZERO) {
1183  case kChi : {
1184  chi2A = fChi2A->At(VZEROcentralityBin);
1185  chi2C = fChi2C->At(VZEROcentralityBin);
1186  chi3A = fChi3A->At(VZEROcentralityBin);
1187  chi3C = fChi3C->At(VZEROcentralityBin);
1188  } break;
1189  case kSigmaSquared : {
1190  chi2A = fSigma2A->At(VZEROcentralityBin);
1191  chi2C = fSigma2C->At(VZEROcentralityBin);
1192  chi3A = fSigma3A->At(VZEROcentralityBin);
1193  chi3C = fSigma3C->At(VZEROcentralityBin);
1194  chi2A = (chi2A > 0) ? 1./chi2A : 1.;
1195  chi2C = (chi2C > 0) ? 1./chi2C : 1.;
1196  chi3A = (chi3A > 0) ? 1./chi3A : 1.;
1197  chi3C = (chi3C > 0) ? 1./chi3C : 1.;
1198  } break;
1199  default : break;
1200  }
1201 
1202  // combine the vzera and vzeroc signal
1203  Q2[0] = chi2A*chi2A*qx2a+chi2C*chi2C*qx2c;
1204  Q2[1] = chi2A*chi2A*qy2a+chi2C*chi2C*qy2c;
1205  Q3[0] = chi3A*chi3A*qx3a+chi3C*chi3C*qx3c;
1206  Q3[1] = chi3A*chi3A*qy3a+chi3C*chi3C*qy3c;
1207 
1208  comb[0] = .5*TMath::ATan2(Q2[1], Q2[0]);
1209  comb[1] = (1./3.)*TMath::ATan2(Q3[1], Q3[0]);
1210 
1211  Double_t _chi(0), _sigma(0), _none(0);
1212  // if requested do the EP correlation histos
1214  switch (fWeightForVZERO) {
1215  case kNone : {
1216  chi2A = fChi2A->At(VZEROcentralityBin);
1217  chi2C = fChi2C->At(VZEROcentralityBin);
1218  _chi = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1219  chi2A = fSigma2A->At(VZEROcentralityBin);
1220  chi2C = fSigma2C->At(VZEROcentralityBin);
1221  chi2A = (chi2A > 0) ? 1./chi2A : 1.;
1222  chi2C = (chi2C > 0) ? 1./chi2C : 1.;
1223  _sigma = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1224  fHistEPCorrelations[fInCentralitySelection]->Fill(.5*TMath::ATan2(qy2a+qy2c,qx2a+qx2c), _chi, _sigma);
1225  } break;
1226  case kChi : {
1227  _chi = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1228  chi2A = fSigma2A->At(VZEROcentralityBin);
1229  chi2C = fSigma2C->At(VZEROcentralityBin);
1230  chi2A = (chi2A > 0) ? 1./chi2A : 1.;
1231  chi2C = (chi2C > 0) ? 1./chi2C : 1.;
1232  _sigma = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1233  fHistEPCorrelations[fInCentralitySelection]->Fill(.5*TMath::ATan2(qy2a+qy2c,qx2a+qx2c), _chi, _sigma);
1234  } break;
1235  case kSigmaSquared : {
1236  _sigma = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1237  chi2A = fChi2A->At(VZEROcentralityBin);
1238  chi2C = fChi2C->At(VZEROcentralityBin);
1239  _chi = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1240  fHistEPCorrelations[fInCentralitySelection]->Fill(.5*TMath::ATan2(qy2a+qy2c,qx2a+qx2c), _chi, _sigma);
1241  } break;
1242  default : break;
1243  }
1244  _none = .5*TMath::ATan2(qy2a+qy2c,qx2a+qx2c);
1245  fHistEPCorrAvChi[fInCentralitySelection]->Fill(_none, _chi);
1246  fHistEPCorrAvSigma[fInCentralitySelection]->Fill(_none, _sigma);
1247  fHistEPCorrChiSigma[fInCentralitySelection]->Fill(_chi, _sigma);
1248  }
1249  }
1250  }
1251 }
1252 //_____________________________________________________________________________
1254 {
1255  // grab the TPC event plane
1256  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1257  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1258  #endif
1259  fNAcceptedTracks = 0; // reset the track counter
1260  Double_t qx2(0), qy2(0); // for psi2
1261  Double_t qx3(0), qy3(0); // for psi3
1262  if(fTracksCont) {
1263  Float_t excludeInEta = -999;
1264  if(fExcludeLeadingJetsFromFit > 0 ) { // remove the leading jet from ep estimate
1265  if(fLeadingJet) excludeInEta = fLeadingJet->Eta();
1266  }
1267  for(Int_t iTPC(0); iTPC < fTracksCont->GetNEntries(); iTPC++) {
1268  AliVParticle* track = fTracksCont->GetParticle(iTPC);
1269  if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
1270  if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta) < GetJetContainer()->GetJetRadius()*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - GetJetContainer()->GetJetRadius() - GetJetContainer()->GetJetEtaMax() ) > 0 )) continue;
1271  fNAcceptedTracks++;
1272  qx2+= TMath::Cos(2.*track->Phi());
1273  qy2+= TMath::Sin(2.*track->Phi());
1274  qx3+= TMath::Cos(3.*track->Phi());
1275  qy3+= TMath::Sin(3.*track->Phi());
1276  }
1277  }
1278  tpc[0] = .5*TMath::ATan2(qy2, qx2);
1279  tpc[1] = (1./3.)*TMath::ATan2(qy3, qx3);
1280 }
1281 //_____________________________________________________________________________
1283 {
1284  // fill the profiles for the resolution parameters
1285  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1286  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1287  #endif
1288  fProfV2Resolution[fInCentralitySelection]->Fill(2., TMath::Cos(2.*(vzero[0][0] - vzero[1][0])));
1289  fProfV2Resolution[fInCentralitySelection]->Fill(3., TMath::Cos(2.*(vzero[1][0] - vzero[0][0])));
1290  fProfV2Resolution[fInCentralitySelection]->Fill(4., TMath::Cos(2.*(vzero[0][0] - tpc[0])));
1291  fProfV2Resolution[fInCentralitySelection]->Fill(5., TMath::Cos(2.*(tpc[0] - vzero[0][0])));
1292  fProfV2Resolution[fInCentralitySelection]->Fill(6., TMath::Cos(2.*(vzero[1][0] - tpc[0])));
1293  fProfV2Resolution[fInCentralitySelection]->Fill(7., TMath::Cos(2.*(tpc[0] - vzero[1][0])));
1294  fProfV3Resolution[fInCentralitySelection]->Fill(2., TMath::Cos(3.*(vzero[0][0] - vzero[1][0])));
1295  fProfV3Resolution[fInCentralitySelection]->Fill(3., TMath::Cos(3.*(vzero[1][0] - vzero[0][0])));
1296  fProfV3Resolution[fInCentralitySelection]->Fill(4., TMath::Cos(3.*(vzero[0][0] - tpc[0])));
1297  fProfV3Resolution[fInCentralitySelection]->Fill(5., TMath::Cos(3.*(tpc[0] - vzero[0][0])));
1298  fProfV3Resolution[fInCentralitySelection]->Fill(6., TMath::Cos(3.*(vzero[1][0] - tpc[0])));
1299  fProfV3Resolution[fInCentralitySelection]->Fill(7., TMath::Cos(3.*(tpc[0] - vzero[1][0])));
1300  // for the resolution of the combined vzero event plane, use two tpc halves as uncorrelated subdetectors
1301  Double_t qx2a(0), qy2a(0); // for psi2a, negative eta
1302  Double_t qx3a(0), qy3a(0); // for psi3a, negative eta
1303  Double_t qx2b(0), qy2b(0); // for psi2a, positive eta
1304  Double_t qx3b(0), qy3b(0); // for psi3a, positive eta
1305  if(fTracks) {
1306  Int_t iTracks(fTracks->GetEntriesFast());
1307  for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
1308  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
1309  if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
1310  if(track->Eta() < 0 ) {
1311  qx2a+= TMath::Cos(2.*track->Phi());
1312  qy2a+= TMath::Sin(2.*track->Phi());
1313  qx3a+= TMath::Cos(3.*track->Phi());
1314  qy3a+= TMath::Sin(3.*track->Phi());
1315  } else if (track->Eta() > 0) {
1316  qx2b+= TMath::Cos(2.*track->Phi());
1317  qy2b+= TMath::Sin(2.*track->Phi());
1318  qx3b+= TMath::Cos(3.*track->Phi());
1319  qy3b+= TMath::Sin(3.*track->Phi());
1320  }
1321  }
1322  }
1323  Double_t tpca2(.5*TMath::ATan2(qy2a, qx2a));
1324  Double_t tpca3((1./3.)*TMath::ATan2(qy3a, qx3a));
1325  Double_t tpcb2(.5*TMath::ATan2(qy2b, qx2b));
1326  Double_t tpcb3((1./3.)*TMath::ATan2(qy3b, qx3b));
1327  fProfV2Resolution[fInCentralitySelection]->Fill(8., TMath::Cos(2.*(vzeroComb[0] - tpca2)));
1328  fProfV2Resolution[fInCentralitySelection]->Fill(9., TMath::Cos(2.*(vzeroComb[0] - tpcb2)));
1329  fProfV2Resolution[fInCentralitySelection]->Fill(10., TMath::Cos(2.*(tpca2 - tpcb2)));
1330  fProfV3Resolution[fInCentralitySelection]->Fill(8., TMath::Cos(3.*(vzeroComb[1] - tpca3)));
1331  fProfV3Resolution[fInCentralitySelection]->Fill(9., TMath::Cos(3.*(vzeroComb[1] - tpcb3)));
1332  fProfV3Resolution[fInCentralitySelection]->Fill(10., TMath::Cos(3.*(tpca3 - tpcb3)));
1333 }
1334 //_____________________________________________________________________________
1336 {
1337  // return the calibrated 2nd and 3rd order q-vectors for vzeroa and vzeroc
1338  // function takes arrays as arguments, which correspond to vzero info in the following way
1339  //
1340  // Qa2[0] = Qx2 for vzero A Qa2[1] = Qy2 for vzero A (etc)
1341 
1342  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1343  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1344  #endif
1345  // placeholders
1346  Double_t phi(-999.), mult(-999.);
1347  // reset placeholders for Q-vector components
1348  Qa2[0] = 0.; Qc2[0] = 0.; Qa3[0] = 0.; Qc3[0] = 0.;
1349  Qa2[1] = 0.; Qc2[1] = 0.; Qa3[1] = 0.; Qc3[1] = 0.;
1350  // for qa purposes, save also raw signal
1351  Double_t QaX(0), QaY(0), QcX(0), QcY(0);
1352  for(Int_t i(0); i < 64; i++) {
1353  // loop over all scintillators, construct Q-vectors in the same loop
1354  phi = TMath::PiOver4()*(0.5+i%8);
1355  mult = InputEvent()->GetVZEROData()->GetMultiplicity(i);
1356  if(fFillQAHistograms) fHistMultVsCellBC->Fill(i, mult);
1357  // note that disabled rings have already been excluded in ReadVZEROCalibration2010h
1358  if(i < 32) { // v0c side
1359  // fill Q-vectors for v0c side
1360  Qc2[0] += mult*TMath::Cos(2.*phi)*fVZEROCpol/fVZEROgainEqualization->GetBinContent(1+i);
1361  Qc3[0] += mult*TMath::Cos(3.*phi)*fVZEROCpol/fVZEROgainEqualization->GetBinContent(1+i);
1362  Qc2[1] += mult*TMath::Sin(2.*phi)*fVZEROCpol/fVZEROgainEqualization->GetBinContent(1+i);
1363  Qc3[1] += mult*TMath::Sin(3.*phi)*fVZEROCpol/fVZEROgainEqualization->GetBinContent(1+i);
1364  if(fFillQAHistograms) {
1365  fHistMultVsCell->Fill(i, mult*fVZEROCpol/fVZEROgainEqualization->GetBinContent(1+i));
1366  QcX += mult*TMath::Cos(2.*phi);
1367  QcY += mult*TMath::Sin(2.*phi);
1368  }
1369  } else { // v0a side
1370  // fill Q-vectors for v0a side
1371  Qa2[0] += mult*TMath::Cos(2.*phi)*fVZEROApol/fVZEROgainEqualization->GetBinContent(1+i);
1372  Qa3[0] += mult*TMath::Cos(3.*phi)*fVZEROApol/fVZEROgainEqualization->GetBinContent(1+i);
1373  Qa2[1] += mult*TMath::Sin(2.*phi)*fVZEROApol/fVZEROgainEqualization->GetBinContent(1+i);
1374  Qa3[1] += mult*TMath::Sin(3.*phi)*fVZEROApol/fVZEROgainEqualization->GetBinContent(1+i);
1375  if(fFillQAHistograms) {
1376  fHistMultVsCell->Fill(i, mult*fVZEROApol/fVZEROgainEqualization->GetBinContent(1+i));
1377  QaX += mult*TMath::Cos(2.*phi);
1378  QaY += mult*TMath::Sin(2.*phi);
1379  }
1380  }
1381  }
1382  // get the cache index and read the correction terms from the cache
1383  Int_t VZEROcentralityBin(GetVZEROCentralityBin());
1384 
1385  if(fFillQAHistograms) {
1386  // recentering qa
1387  fHistQxV0aBC->Fill(Qa2[0], VZEROcentralityBin);
1388  fHistQyV0aBC->Fill(Qa2[1], VZEROcentralityBin);
1389  fHistQxV0cBC->Fill(Qc2[0], VZEROcentralityBin);
1390  fHistQyV0cBC->Fill(Qc2[0], VZEROcentralityBin);
1391  fHistEPBC->Fill(.5*TMath::ATan2(QaY+QcY, QaX+QcX));
1392  }
1393 
1394  Double_t Qx2amean = fMeanQ[VZEROcentralityBin][1][0];
1395  Double_t Qx2arms = fWidthQ[VZEROcentralityBin][1][0];
1396  Double_t Qy2amean = fMeanQ[VZEROcentralityBin][1][1];
1397  Double_t Qy2arms = fWidthQ[VZEROcentralityBin][1][1];
1398 
1399  Double_t Qx2cmean = fMeanQ[VZEROcentralityBin][0][0];
1400  Double_t Qx2crms = fWidthQ[VZEROcentralityBin][0][0];
1401  Double_t Qy2cmean = fMeanQ[VZEROcentralityBin][0][1];
1402  Double_t Qy2crms = fWidthQ[VZEROcentralityBin][0][1];
1403 
1404  Double_t Qx3amean = fMeanQv3[VZEROcentralityBin][1][0];
1405  Double_t Qx3arms = fWidthQv3[VZEROcentralityBin][1][0];
1406  Double_t Qy3amean = fMeanQv3[VZEROcentralityBin][1][1];
1407  Double_t Qy3arms = fWidthQv3[VZEROcentralityBin][1][1];
1408 
1409  Double_t Qx3cmean = fMeanQv3[VZEROcentralityBin][0][0];
1410  Double_t Qx3crms = fWidthQv3[VZEROcentralityBin][0][0];
1411  Double_t Qy3cmean = fMeanQv3[VZEROcentralityBin][0][1];
1412  Double_t Qy3crms = fWidthQv3[VZEROcentralityBin][0][1];
1413 
1414  // update the weighted q-vectors with the re-centered values
1415  Qa2[0] = (Qa2[0] - Qx2amean)/Qx2arms;
1416  Qa2[1] = (Qa2[1] - Qy2amean)/Qy2arms;
1417  Qc2[0] = (Qc2[0] - Qx2cmean)/Qx2crms;
1418  Qc2[1] = (Qc2[1] - Qy2cmean)/Qy2crms;
1419 
1420  Qa3[0] = (Qa3[0] - Qx3amean)/Qx3arms;
1421  Qa3[1] = (Qa3[1] - Qy3amean)/Qy3arms;
1422  Qc3[0] = (Qc3[0] - Qx3cmean)/Qx3crms;
1423  Qc3[1] = (Qc3[1] - Qy3cmean)/Qy3crms;
1424 
1425  if(fFillQAHistograms) {
1426  // recentering qa
1427  fHistQxV0a->Fill(Qa2[0], VZEROcentralityBin);
1428  fHistQyV0a->Fill(Qa2[1], VZEROcentralityBin);
1429  fHistQxV0c->Fill(Qc2[0], VZEROcentralityBin);
1430  fHistQyV0c->Fill(Qc2[0], VZEROcentralityBin);
1431  fHistEP->Fill(.5*TMath::ATan2(Qa2[1]+Qc2[1], Qa2[0]+Qc2[0]));
1432  }
1433 }
1434 //_____________________________________________________________________________
1436 {
1437  // calculate calibrated q-vector of the combined vzeroa, vzeroc system
1438  // this is somewhat ugly as CalculateQvectorCombinedVZERO is called more than once per event
1439  // but for now it will have to do ...
1440  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1441  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1442  #endif
1443 
1444  // first step: retrieve the q-vectors component-wise per vzero detector
1445  Double_t QA2[] = {-999., -999.};
1446  Double_t QA3[] = {-999., -999.};
1447  Double_t QC2[] = {-999., -999.};
1448  Double_t QC3[] = {-999., -999.};
1449  CalculateQvectorVZERO(QA2, QC2, QA3, QC3);
1450 
1451  // get cache index and retrieve the chi weights for this centrality
1452  Int_t VZEROcentralityBin(GetVZEROCentralityBin());
1453  Double_t chi2A(1);
1454  Double_t chi2C(1);
1455  Double_t chi3A(1);
1456  Double_t chi3C(1);
1457 
1458  switch (fWeightForVZERO) {
1459  case kChi : {
1460  chi2A = fChi2A->At(VZEROcentralityBin);
1461  chi2C = fChi2C->At(VZEROcentralityBin);
1462  chi3A = fChi3A->At(VZEROcentralityBin);
1463  chi3C = fChi3C->At(VZEROcentralityBin);
1464  } break;
1465  case kSigmaSquared : {
1466  chi2A = fSigma2A->At(VZEROcentralityBin);
1467  chi2C = fSigma2C->At(VZEROcentralityBin);
1468  chi3A = fSigma3A->At(VZEROcentralityBin);
1469  chi3C = fSigma3C->At(VZEROcentralityBin);
1470  chi2A = (chi2A > 0) ? 1./chi2A : 1.;
1471  chi2C = (chi2C > 0) ? 1./chi2C : 1.;
1472  chi3A = (chi3A > 0) ? 1./chi3A : 1.;
1473  chi3C = (chi3C > 0) ? 1./chi3C : 1.;
1474  } break;
1475  default : break;
1476  }
1477 
1478  // bookkkeep these guys
1479  Double_t qx2a(QA2[0]), qy2a(QA2[1]), qx2c(QC2[0]), qy2c(QC2[1]);
1480  // combine the vzera and vzeroc signal
1481  Q2[0] = chi2A*chi2A*QA2[0]+chi2C*chi2C*QC2[0];
1482  Q2[1] = chi2A*chi2A*QA2[1]+chi2C*chi2C*QC2[1];
1483  Q3[0] = chi3A*chi3A*QA3[0]+chi3C*chi3C*QC3[0];
1484  Q3[1] = chi3A*chi3A*QA3[1]+chi3C*chi3C*QC3[1];
1485 
1486  Double_t _chi(0), _sigma(0), _none(0);
1487  // if requested do the EP correlation histos
1489  switch (fWeightForVZERO) {
1490  case kNone : {
1491  chi2A = fChi2A->At(VZEROcentralityBin);
1492  chi2C = fChi2C->At(VZEROcentralityBin);
1493  _chi = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1494  chi2A = fSigma2A->At(VZEROcentralityBin);
1495  chi2C = fSigma2C->At(VZEROcentralityBin);
1496  chi2A = (chi2A > 0) ? 1./chi2A : 1.;
1497  chi2C = (chi2C > 0) ? 1./chi2C : 1.;
1498  _sigma = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1499  fHistEPCorrelations[fInCentralitySelection]->Fill(.5*TMath::ATan2(qy2a+qy2c,qx2a+qx2c), _chi, _sigma);
1500  } break;
1501  case kChi : {
1502  _chi = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1503  chi2A = fSigma2A->At(VZEROcentralityBin);
1504  chi2C = fSigma2C->At(VZEROcentralityBin);
1505  chi2A = (chi2A > 0) ? 1./chi2A : 1.;
1506  chi2C = (chi2C > 0) ? 1./chi2C : 1.;
1507  _sigma = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1508  fHistEPCorrelations[fInCentralitySelection]->Fill(.5*TMath::ATan2(qy2a+qy2c,qx2a+qx2c), _chi, _sigma);
1509  } break;
1510  case kSigmaSquared : {
1511  _sigma = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1512  chi2A = fChi2A->At(VZEROcentralityBin);
1513  chi2C = fChi2C->At(VZEROcentralityBin);
1514  _chi = .5*TMath::ATan2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c, chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
1515  fHistEPCorrelations[fInCentralitySelection]->Fill(.5*TMath::ATan2(qy2a+qy2c,qx2a+qx2c), _chi, _sigma);
1516  } break;
1517  default : break;
1518  }
1519  _none = .5*TMath::ATan2(qy2a+qy2c,qx2a+qx2c);
1520  fHistEPCorrAvChi[fInCentralitySelection]->Fill(_none, _chi);
1521  fHistEPCorrAvSigma[fInCentralitySelection]->Fill(_none, _sigma);
1522  fHistEPCorrChiSigma[fInCentralitySelection]->Fill(_chi, _sigma);
1523  }
1524 }
1525 //_____________________________________________________________________________
1527  AliParticleContainer* tracksCont, AliClusterContainer* clusterCont, AliEmcalJet* jet) const
1528 {
1529  // get a random cone
1530  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_2
1531  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1532  #endif
1533  pt = 0; eta = 0; phi = 0;
1534  Float_t etaJet(999), phiJet(999), dJet(999); // no jet: same as jet very far away
1535  if(jet) { // if a leading jet is given, use its kinematic properties to exclude it
1536  etaJet = jet->Eta();
1537  phiJet = jet->Phi();
1538  }
1539  // the random cone acceptance has to equal the jet acceptance
1540  // this also insures safety when runnnig on the semi-good tpc runs for 11h data,
1541  // where jet acceptance is adjusted to reduced acceptance - hence random cone acceptance as well
1542  Float_t minPhi(GetJetContainer()->GetJetPhiMin()), maxPhi(GetJetContainer()->GetJetPhiMax());
1543  if(maxPhi > TMath::TwoPi()) maxPhi = TMath::TwoPi();
1544  if(minPhi < 0 ) minPhi = 0.;
1545  // construct a random cone and see if it's far away enough from the leading jet
1546  Int_t attempts(1000);
1547  while(kTRUE) {
1548  attempts--;
1549  eta = gRandom->Uniform(GetJetContainer()->GetJetEtaMin(), GetJetContainer()->GetJetEtaMax());
1550  phi = gRandom->Uniform(minPhi, maxPhi);
1551 
1552  dJet = TMath::Sqrt((etaJet-eta)*(etaJet-eta)+(phiJet-phi)*(phiJet-phi));
1553  if(dJet > fMinDisanceRCtoLJ) break;
1554  else if (attempts == 0) {
1555  printf(" > No random cone after 1000 tries, giving up ... !\n");
1556  return;
1557  }
1558  }
1559  // get the charged energy (if tracks are provided)
1560  if(tracksCont) {
1561  tracksCont->ResetCurrentID();
1562  AliVParticle* track = tracksCont->GetNextAcceptParticle();
1563  while(track) {
1564  Float_t etaTrack(track->Eta()), phiTrack(track->Phi());
1565  // get distance from cone
1566  if(TMath::Abs(phiTrack-phi) > TMath::Abs(phiTrack - phi + TMath::TwoPi())) phiTrack+=TMath::TwoPi();
1567  if(TMath::Abs(phiTrack-phi) > TMath::Abs(phiTrack - phi - TMath::TwoPi())) phiTrack-=TMath::TwoPi();
1568  if(TMath::Sqrt(TMath::Abs((etaTrack-eta)*(etaTrack-eta)+(phiTrack-phi)*(phiTrack-phi))) <= GetJetRadius()) pt += track->Pt();
1569  track = tracksCont->GetNextAcceptParticle();
1570  }
1571  }
1572  // get the neutral energy (if clusters are provided)
1573  if(clusterCont) {
1574  TLorentzVector momentum;
1575  clusterCont->ResetCurrentID();
1576  AliVCluster* cluster = clusterCont->GetNextAcceptCluster();
1577  while(cluster) {
1578  cluster->GetMomentum(momentum, const_cast<Double_t*>(fVertex));
1579  Float_t etaClus(momentum.Eta()), phiClus(momentum.Phi());
1580  // get distance from cone
1581  if(TMath::Abs(phiClus-phi) > TMath::Abs(phiClus - phi + TMath::TwoPi())) phiClus+=TMath::TwoPi();
1582  if(TMath::Abs(phiClus-phi) > TMath::Abs(phiClus - phi - TMath::TwoPi())) phiClus-=TMath::TwoPi();
1583  if(TMath::Sqrt(TMath::Abs((etaClus-eta)*(etaClus-eta)+(phiClus-phi)*(phiClus-phi))) <= GetJetRadius()) pt += momentum.Pt();
1584  cluster = clusterCont->GetNextAcceptCluster();
1585  }
1586  }
1587 }
1588 //_____________________________________________________________________________
1590  // get the second order q-cumulant, a -999 return will be caught in the qa routine of CorrectRho
1591  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1592  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1593  #endif
1594  Double_t reQ(0), imQ(0), modQ(0), M11(0), M(0);
1595  if(fUsePtWeight) { // for the weighted 2-nd order q-cumulant
1596  QCnQnk(harm, 1, reQ, imQ); // get the weighted 2-nd order q-vectors
1597  modQ = reQ*reQ+imQ*imQ; // get abs Q-squared
1598  M11 = QCnM11(); // equals S2,1 - S1,2
1599  return (M11 > 0) ? ((modQ - QCnS(1,2))/M11) : -999;
1600  } // else return the non-weighted 2-nd order q-cumulant
1601  QCnQnk(harm, 0, reQ, imQ); // get the non-weighted 2-nd order q-vectors
1602  modQ = reQ*reQ+imQ*imQ; // get abs Q-squared
1603  M = QCnM();
1604  return (M > 1) ? (modQ - M)/(M*(M-1)) : -999;
1605 }
1606 //_____________________________________________________________________________
1608  // get the fourth order q-cumulant, a -999 return will be caught in the qa routine of CorrectRho
1609  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1610  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1611  #endif
1612  Double_t reQn1(0), imQn1(0), reQ2n2(0), imQ2n2(0), reQn3(0), imQn3(0), M1111(0), M(0);
1613  Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0); // terms of the calculation
1614  if(fUsePtWeight) { // for the weighted 4-th order q-cumulant
1615  QCnQnk(harm, 1, reQn1, imQn1);
1616  QCnQnk(harm*2, 2, reQ2n2, imQ2n2);
1617  QCnQnk(harm, 3, reQn3, imQn3);
1618  // fill in the terms ...
1619  a = (reQn1*reQn1+imQn1*imQn1)*(reQn1*reQn1+imQn1*imQn1);
1620  b = reQ2n2*reQ2n2 + imQ2n2*imQ2n2;
1621  c = -2.*(reQ2n2*reQn1*reQn1-reQ2n2*imQn1*imQn1+2.*imQ2n2*reQn1*imQn1);
1622  d = 8.*(reQn3*reQn1+imQn3*imQn1);
1623  e = -4.*QCnS(1,2)*(reQn1*reQn1+imQn1*imQn1);
1624  f = -6.*QCnS(1,4);
1625  g = 2.*QCnS(2,2);
1626  M1111 = QCnM1111();
1627  return (M1111 > 0) ? (a+b+c+d+e+f+g)/M1111 : -999;
1628  } // else return the unweighted case
1629  Double_t reQn(0), imQn(0), reQ2n(0), imQ2n(0);
1630  QCnQnk(harm, 0, reQn, imQn);
1631  QCnQnk(harm*2, 0, reQ2n, imQ2n);
1632  // fill in the terms ...
1633  M = QCnM();
1634  if(M < 4) return -999;
1635  a = (reQn*reQn+imQn*imQn)*(reQn*reQn+imQn*imQn);
1636  b = reQ2n*reQ2n + imQ2n*imQ2n;
1637  c = -2.*(reQ2n*reQn*reQn-reQ2n*imQn*imQn+2.*imQ2n*reQn*imQn);
1638  e = -4.*(M-2)*(reQn*reQn+imQn*imQn);
1639  f = 2.*M*(M-3);
1640  return (a+b+c+e+f)/(M*(M-1)*(M-2)*(M-3));
1641 }
1642 //_____________________________________________________________________________
1644  // get the weighted n-th order q-vector, pass real and imaginary part as reference
1645  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1646  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1647  #endif
1648  if(!fTracks) return;
1649  fNAcceptedTracksQCn = 0;
1650  Int_t iTracks(fTracks->GetEntriesFast());
1651  for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
1652  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
1653  if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
1655  // for the unweighted case, k equals zero and the weight doesn't contribute to the equation below
1656  reQ += TMath::Power(track->Pt(), k) * TMath::Cos(((double)n)*track->Phi());
1657  imQ += TMath::Power(track->Pt(), k) * TMath::Sin(((double)n)*track->Phi());
1658  }
1659 }
1660 //_____________________________________________________________________________
1662  TClonesArray* pois, TArrayD* ptBins, Bool_t vpart, Double_t* repn, Double_t* impn,
1663  Double_t *mp, Double_t *reqn, Double_t *imqn, Double_t* mq, Int_t n)
1664 {
1665  // get unweighted differential flow vectors
1666  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1667  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1668  #endif
1669  Int_t iPois(pois->GetEntriesFast());
1670  if(vpart) {
1671  for(Int_t i(0); i < iPois; i++) {
1672  for(Int_t ptBin(0); ptBin < ptBins->GetSize()-1; ptBin++) {
1673  AliVTrack* poi = static_cast<AliVTrack*>(pois->At(i));
1674  if(PassesCuts(poi)) {
1675  if(poi->Pt() >= ptBins->At(ptBin) && poi->Pt() < ptBins->At(ptBin+1)) {
1676  // fill the flow vectors assuming that all poi's are in the rp selection (true by design)
1677  repn[ptBin]+=TMath::Cos(((double)n)*poi->Phi());
1678  impn[ptBin]+=TMath::Sin(((double)n)*poi->Phi());
1679  mp[ptBin]++;
1680  reqn[ptBin]+=TMath::Cos(((double)n)*poi->Phi());
1681  imqn[ptBin]+=TMath::Sin(((double)n)*poi->Phi());
1682  mq[ptBin]++;
1683  }
1684  }
1685  }
1686  }
1687  } else {
1688  for(Int_t i(0); i < iPois; i++) {
1689  for(Int_t ptBin(0); ptBin < ptBins->GetSize()-1; ptBin++) {
1690  AliEmcalJet* poi = static_cast<AliEmcalJet*>(pois->At(i));
1691  if(PassesCuts(poi)) {
1692  Double_t pt(poi->Pt()-poi->Area()*fLocalRho->GetLocalVal(poi->Phi(), GetJetContainer()->GetJetRadius(), fLocalRho->GetVal()));
1693  if(fUse2DIntegration) pt = poi->Pt()-poi->Area()*fLocalRho->GetLocalValInEtaPhi(poi->Phi(), GetJetContainer()->GetJetRadius(), fLocalRho->GetVal());
1694  if(pt >= ptBins->At(ptBin) && pt < ptBins->At(ptBin+1)) {
1695  repn[ptBin]+=TMath::Cos(((double)n)*poi->Phi());
1696  impn[ptBin]+=TMath::Sin(((double)n)*poi->Phi());
1697  mp[ptBin]++; // qn isn't filled, no overlap between poi's and rp's
1698  }
1699  }
1700  }
1701  }
1702  }
1703 }
1704 //_____________________________________________________________________________
1706  // get the weighted ij-th order autocorrelation correction
1707  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1708  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1709  #endif
1710  if(!fTracks || i <= 0 || j <= 0) return -999;
1711  Int_t iTracks(fTracks->GetEntriesFast());
1712  Double_t Sij(0);
1713  for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
1714  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
1715  if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
1716  Sij+=TMath::Power(track->Pt(), j);
1717  }
1718  return TMath::Power(Sij, i);
1719 }
1720 //_____________________________________________________________________________
1722  // get multiplicity for unweighted q-cumulants. function QCnQnk should be called first
1723  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1724  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1725  #endif
1726  return (Double_t) fNAcceptedTracksQCn;
1727 }
1728 //_____________________________________________________________________________
1730  // get multiplicity weights for the weighted two particle cumulant
1731  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1732  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1733  #endif
1734  return (QCnS(2,1) - QCnS(1,2));
1735 }
1736 //_____________________________________________________________________________
1738  // get multiplicity weights for the weighted four particle cumulant
1739  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1740  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1741  #endif
1742  return (QCnS(4,1)-6*QCnS(1,2)*QCnS(2,1)+8*QCnS(1,3)*QCnS(1,1)+3*QCnS(2,2)-6*QCnS(1,4));
1743 }
1744 //_____________________________________________________________________________
1746  // decides how to deal with the situation where c2 or c3 is negative
1747  // returns kTRUE depending on whether or not a modulated rho is used for the jet background
1748  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1749  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1750  #endif
1751  if(TMath::AreEqualAbs(fFitModulation->GetParameter(3), .0, 1e-10) && TMath::AreEqualAbs(fFitModulation->GetParameter(7), .0,1e-10)) {
1752  fFitModulation->SetParameter(7, 0);
1753  fFitModulation->SetParameter(3, 0);
1754  fFitModulation->SetParameter(0, fLocalRho->GetVal());
1755  return kTRUE; // v2 and v3 have physical null values
1756  }
1757  switch (fQCRecovery) {
1758  case kFixedRho : { // roll back to the original rho
1759  fFitModulation->SetParameter(7, 0);
1760  fFitModulation->SetParameter(3, 0);
1761  fFitModulation->SetParameter(0, fLocalRho->GetVal());
1762  return kFALSE; // rho is forced to be fixed
1763  }
1764  case kNegativeVn : {
1765  Double_t c2(fFitModulation->GetParameter(3));
1766  Double_t c3(fFitModulation->GetParameter(7));
1767  if( c2 < 0 ) c2 = -1.*TMath::Sqrt(-1.*c2);
1768  if( c3 < 0 ) c3 = -1.*TMath::Sqrt(-1.*c3);
1769  fFitModulation->SetParameter(3, c2);
1770  fFitModulation->SetParameter(7, c3);
1771  return kTRUE; // is this a physical quantity ?
1772  }
1773  case kTryFit : {
1774  fitModulationType tempType(fFitModulationType); // store temporarily
1776  fFitModulation->SetParameter(7, 0);
1777  fFitModulation->SetParameter(3, 0);
1778  Bool_t pass(CorrectRho(psi2, psi3)); // do the fit and all quality checks
1779  fFitModulationType = tempType; // roll back for next event
1780  return pass;
1781  }
1782  default : return kFALSE;
1783  }
1784  return kFALSE;
1785 }
1786 //_____________________________________________________________________________
1788 {
1789  // get rho' -> rho(phi)
1790  // two routines are available, both can be used with or without pt weights
1791  // [1] get vn from q-cumulants or as an integrated value from a user supplied histogram
1792  // in case of cumulants, both cumulants and vn values are stored. in both cases, v2 and v3
1793  // are expected. a check is performed to see if rho has no negative local minimum
1794  // for full description, see Phys. Rev. C 83, 044913
1795  // since the cn distribution has negative values, vn = sqrt(cn) can be imaginary sometimes
1796  // in this case one can either roll back to the 'original' rixed rho, do a fit for vn or take use
1797  // vn = - sqrt(|cn|)
1798  // [2] fitting a fourier expansion to the de/dphi distribution
1799  // the fit can be done with either v2, v3 or a combination.
1800  // in all cases, a cut can be made on the p-value of the chi-squared value of the fit
1801  // and a check can be performed to see if rho has no negative local minimum
1802  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
1803  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1804  #endif
1805  Int_t freeParams(2); // free parameters of the fit (for NDF)
1806  switch (fFitModulationType) { // for approaches where no fitting is required
1807  case kQC2 : {
1808  fFitModulation->FixParameter(4, psi2);
1809  fFitModulation->FixParameter(6, psi3);
1810  fFitModulation->FixParameter(3, CalculateQC2(2)); // set here with cn, vn = sqrt(cn)
1811  fFitModulation->FixParameter(7, CalculateQC2(3));
1812  // first fill the histos of the raw cumulant distribution
1813  if (fUsePtWeight) { // use weighted weights
1814  Double_t dQCnM11 = (fNoEventWeightsForQC) ? 1. : QCnM11();
1815  fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), dQCnM11);
1816  fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), dQCnM11);
1817  } else {
1818  Double_t dQCnM = (fNoEventWeightsForQC) ? 2. : QCnM();
1819  fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), dQCnM*(dQCnM-1));
1820  fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), dQCnM*(dQCnM-1));
1821  }
1822  // then see if one of the cn value is larger than zero and vn is readily available
1823  if(fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) {
1824  fFitModulation->FixParameter(3, TMath::Sqrt(fFitModulation->GetParameter(3)));
1825  fFitModulation->FixParameter(7, TMath::Sqrt(fFitModulation->GetParameter(7)));
1826  } else if (!QCnRecovery(psi2, psi3)) return kFALSE; // try to recover the cumulant, this will set v2 and v3
1827  if(fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) { // general check
1828  fFitModulation->SetParameter(7, 0);
1829  fFitModulation->SetParameter(3, 0);
1830  fFitModulation->SetParameter(0, fLocalRho->GetVal());
1831  return kFALSE;
1832  }
1833  return kTRUE;
1834  } break;
1835  case kQC4 : {
1836  fFitModulation->FixParameter(4, psi2);
1837  fFitModulation->FixParameter(6, psi3);
1838  fFitModulation->FixParameter(3, CalculateQC4(2)); // set here with cn, vn = sqrt(cn)
1839  fFitModulation->FixParameter(7, CalculateQC4(3));
1840  // first fill the histos of the raw cumulant distribution
1841  if (fUsePtWeight) { // use weighted weights
1842  fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM1111()*/);
1843  fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*, QCnM1111()*/);
1844  } else {
1845  fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM1111()*/);
1846  fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*, QCnM1111()*/);
1847  }
1848  // then see if one of the cn value is larger than zero and vn is readily available
1849  if(fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) {
1850  fFitModulation->FixParameter(3, TMath::Sqrt(fFitModulation->GetParameter(3)));
1851  fFitModulation->FixParameter(7, TMath::Sqrt(fFitModulation->GetParameter(7)));
1852  } else if (!QCnRecovery(psi2, psi3)) return kFALSE; // try to recover the cumulant, this will set v2 and v3
1853  if(fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) { // general check
1854  fFitModulation->SetParameter(7, 0);
1855  fFitModulation->SetParameter(3, 0);
1856  fFitModulation->SetParameter(0, fLocalRho->GetVal());
1857  return kFALSE;
1858  }
1859  } break;
1860  case kIntegratedFlow : {
1861  // use v2 and v3 values from an earlier iteration over the data
1862  fFitModulation->FixParameter(3, fUserSuppliedV2->GetBinContent(fUserSuppliedV2->GetXaxis()->FindBin(fCent)));
1863  fFitModulation->FixParameter(4, psi2);
1864  fFitModulation->FixParameter(6, psi3);
1865  fFitModulation->FixParameter(7, fUserSuppliedV3->GetBinContent(fUserSuppliedV3->GetXaxis()->FindBin(fCent)));
1866  if(fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) {
1867  fFitModulation->SetParameter(7, 0);
1868  fFitModulation->SetParameter(3, 0);
1869  fFitModulation->SetParameter(0, fLocalRho->GetVal());
1870  return kFALSE;
1871  }
1872  return kTRUE;
1873  }
1874  default : break;
1875  }
1876  TString detector("");
1877  switch (fDetectorType) {
1878  case kTPC : detector+="TPC";
1879  break;
1880  case kVZEROA : detector+="VZEROA";
1881  break;
1882  case kVZEROC : detector+="VZEROC";
1883  break;
1884  case kVZEROComb : detector+="VZEROComb";
1885  break;
1886  case kFixedEP : detector+="FixedEP";
1887  break;
1888  default: break;
1889  }
1890  Int_t iTracks(fTracks->GetEntriesFast());
1891  Double_t excludeInEta = -999;
1892  Double_t excludeInPhi = -999;
1893  Double_t excludeInPt = -999;
1894  if(iTracks <= 0 || fLocalRho->GetVal() <= 0 ) return kFALSE; // no use fitting an empty event ...
1895  if(fExcludeLeadingJetsFromFit > 0 ) {
1896  if(fLeadingJet) {
1897  excludeInEta = fLeadingJet->Eta();
1898  excludeInPhi = fLeadingJet->Phi();
1899  excludeInPt = fLeadingJet->Pt();
1900  }
1901  }
1902  // check the acceptance of the track selection that will be used
1903  // if one uses e.g. semi-good tpc tracks, accepance in phi is reduced to 0 < phi < 4
1904  // the defaults (-10 < phi < 10) which accept all, are then overwritten
1905  Double_t lowBound(0.), upBound(TMath::TwoPi()); // bounds for fit
1906  if(GetParticleContainer()->GetParticlePhiMin() > lowBound) lowBound = GetParticleContainer()->GetParticlePhiMin();
1907  if(GetParticleContainer()->GetParticlePhiMax() < upBound) upBound = GetParticleContainer()->GetParticlePhiMax();
1908  fHistSwap->Reset(); // clear the histogram
1909  TH1F _tempSwap; // on stack for quick access
1910  TH1F _tempSwapN; // on stack for quick access, bookkeeping histogram
1912  if(fNAcceptedTracks < 49) fNAcceptedTracks = 49; // avoid aliasing effects
1913  _tempSwap = TH1F("_tempSwap", "_tempSwap", TMath::CeilNint(TMath::Sqrt(fNAcceptedTracks)), lowBound, upBound);
1914  if(fUsePtWeightErrorPropagation) _tempSwapN = TH1F("_tempSwapN", "_tempSwapN", TMath::CeilNint(TMath::Sqrt(fNAcceptedTracks)), lowBound, upBound);
1915  if(fUsePtWeight) _tempSwap.Sumw2();
1916  }
1917  else _tempSwap = *fHistSwap; // now _tempSwap holds the desired histo
1918  // non poissonian error when using pt weights
1919  Double_t totalpts(0.), totalptsquares(0.), totalns(0.);
1920  for(Int_t i(0); i < iTracks; i++) {
1921  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
1922  if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta) < GetJetContainer()->GetJetRadius()*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - GetJetContainer()->GetJetRadius() - GetJetContainer()->GetJetEtaMax() ) > 0 )) continue;
1923  if(!PassesCuts(track) || track->Pt() > fSoftTrackMaxPt || track->Pt() < fSoftTrackMinPt) continue;
1924  if(fUsePtWeight) {
1925  _tempSwap.Fill(track->Phi(), track->Pt());
1927  totalpts += track->Pt();
1928  totalptsquares += track->Pt()*track->Pt();
1929  totalns += 1;
1930  _tempSwapN.Fill(track->Phi());
1931  }
1932  }
1933  else _tempSwap.Fill(track->Phi());
1934  }
1936  // in the case of pt weights overwrite the poissonian error estimate which is assigned by root by a more sophisticated appraoch
1937  // the assumption here is that the bin error will be dominated by the uncertainty in the mean pt in a bin and in the uncertainty
1938  // of the number of tracks in a bin, the first of which will be estimated from the sample standard deviation of all tracks in the
1939  // event, for the latter use a poissonian estimate. the two contrubitions are assumed to be uncorrelated
1940  if(totalns < 2) return kFALSE; // not one track passes the cuts > 2 avoids possible division by 0 later on
1941  for(Int_t l = 0; l < _tempSwap.GetNbinsX(); l++) {
1942  if(_tempSwapN.GetBinContent(l+1) == 0) {
1943  _tempSwap.SetBinContent(l+1,0);
1944  _tempSwap.SetBinError(l+1,0);
1945  }
1946  else {
1947  Double_t vartimesnsq = totalptsquares*totalns - totalpts*totalpts;
1948  Double_t variance = vartimesnsq/(totalns*(totalns-1.));
1949  Double_t SDOMSq = variance / _tempSwapN.GetBinContent(l+1);
1950  Double_t SDOMSqOverMeanSq = SDOMSq * _tempSwapN.GetBinContent(l+1) * _tempSwapN.GetBinContent(l+1) / (_tempSwapN.GetBinContent(l+1) * _tempSwapN.GetBinContent(l+1));
1951  Double_t poissonfrac = 1./_tempSwapN.GetBinContent(l+1);
1952  Double_t vartotalfrac = SDOMSqOverMeanSq + poissonfrac;
1953  Double_t vartotal = vartotalfrac * _tempSwap.GetBinContent(l+1) * _tempSwap.GetBinContent(l+1);
1954  if(vartotal > 0.0001) _tempSwap.SetBinError(l+1,TMath::Sqrt(vartotal));
1955  else {
1956  _tempSwap.SetBinContent(l+1,0);
1957  _tempSwap.SetBinError(l+1,0);
1958  }
1959  }
1960  }
1961  }
1962  fFitModulation->SetParameter(0, fLocalRho->GetVal());
1963  switch (fFitModulationType) {
1964  case kNoFit : {
1965  fFitModulation->FixParameter(0, fLocalRho->GetVal() );
1966  freeParams = 0;
1967  } break;
1968  case kV2 : {
1969  fFitModulation->FixParameter(4, psi2);
1970  freeParams = 1;
1971  } break;
1972  case kV3 : {
1973  fFitModulation->FixParameter(4, psi3);
1974  freeParams = 1;
1975  } break;
1976  case kCombined : {
1977  fFitModulation->FixParameter(4, psi2);
1978  fFitModulation->FixParameter(6, psi3);
1979  freeParams = 2;
1980  } break;
1981  case kFourierSeries : {
1982  // in this approach, an explicit calculation will be made of vn = sqrt(xn^2+yn^2)
1983  // where x[y] = Integrate[r(phi)cos[sin](n phi)dphi, 0, 2pi]
1984  Double_t cos2(0), sin2(0), cos3(0), sin3(0), sumPt(0);
1985  for(Int_t i(0); i < iTracks; i++) {
1986  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
1987  if(!PassesCuts(track) || track->Pt() > fSoftTrackMaxPt || track->Pt() < fSoftTrackMinPt) continue;
1988  sumPt += track->Pt();
1989  cos2 += track->Pt()*TMath::Cos(2*PhaseShift(track->Phi()-psi2));
1990  sin2 += track->Pt()*TMath::Sin(2*PhaseShift(track->Phi()-psi2));
1991  cos3 += track->Pt()*TMath::Cos(3*PhaseShift(track->Phi()-psi3));
1992  sin3 += track->Pt()*TMath::Sin(3*PhaseShift(track->Phi()-psi3));
1993  }
1994  fFitModulation->SetParameter(3, TMath::Sqrt(cos2*cos2+sin2*sin2)/fLocalRho->GetVal());
1995  fFitModulation->SetParameter(4, psi2);
1996  fFitModulation->SetParameter(6, psi3);
1997  fFitModulation->SetParameter(7, TMath::Sqrt(cos3*cos3+sin3*sin3)/fLocalRho->GetVal());
1998  } break;
1999  default : break;
2000  }
2001  if(fRunToyMC) {
2002  // toy mc, just here to check procedure, azimuthal profile is filled from hypothesis so p-value distribution should be flat
2003  Int_t _bins = _tempSwap.GetXaxis()->GetNbins();
2004  TF1* _tempFit = new TF1("temp_fit_kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi());
2005  _tempFit->SetParameter(0, fFitModulation->GetParameter(0)); // normalization
2006  _tempFit->SetParameter(3, 0.1); // v2
2007  _tempFit->FixParameter(1, 1.); // constant
2008  _tempFit->FixParameter(2, 2.); // constant
2009  _tempFit->FixParameter(5, 3.); // constant
2010  _tempFit->FixParameter(4, fFitModulation->GetParameter(4));
2011  _tempFit->FixParameter(6, fFitModulation->GetParameter(6));
2012  _tempFit->SetParameter(7, 0.1); // v3
2013  _tempSwap.Reset(); // rese bin content
2014  for(int _binsI = 0; _binsI < _bins*_bins; _binsI++) _tempSwap.Fill(_tempFit->GetRandom());
2015  }
2016  _tempSwap.Fit(fFitModulation, fFitModulationOptions.Data(), "", lowBound, upBound);
2017  // the quality of the fit is evaluated from 1 - the cdf of the chi square distribution
2018  // three methods are available, all with their drawbacks. all are stored, one is selected to do the cut
2019  Int_t NDF(_tempSwap.GetXaxis()->GetNbins()-freeParams);
2020  if(NDF == 0 || (float)NDF <= 0.) return kFALSE;
2021  Double_t CDF(1.-ChiSquareCDF(NDF, ChiSquare(_tempSwap, fFitModulation)));
2022  Double_t CDFROOT(1.-ChiSquareCDF(NDF, fFitModulation->GetChisquare()));
2023  Double_t CDFKolmogorov(KolmogorovTest(/*_tempSwap, fFitModulation*/));
2024  // fill the values and centrality correlation (redundant but easy on the eyes)
2025  fHistPvalueCDF->Fill(CDF);
2026  fHistPvalueCDFCent->Fill(fCent, CDF);
2027  fHistPvalueCDFROOT->Fill(CDFROOT);
2028  fHistPvalueCDFROOTCent->Fill(fCent, CDFROOT);
2029  fHistKolmogorovTest->Fill(CDFKolmogorov);
2030  fHistChi2ROOTCent->Fill(fCent, fFitModulation->GetChisquare()/((float)NDF));
2031  fHistChi2Cent->Fill(fCent, ChiSquare(_tempSwap, fFitModulation)/((float)NDF));
2032  fHistKolmogorovTestCent->Fill(fCent, CDFKolmogorov);
2033  fHistPChi2Root->Fill(CDFROOT, fFitModulation->GetChisquare()/((float)NDF));
2034  fHistPChi2->Fill(CDF, ChiSquare(_tempSwap, fFitModulation)/((float)NDF));
2035  fHistPKolmogorov->Fill(CDF, CDFKolmogorov);
2036 
2037  // variable CDF is used for making cuts, so we fill it with the selected p-value
2038  switch (fFitGoodnessTest) {
2039  case kChi2ROOT : {
2040  CDF = CDFROOT;
2041  } break;
2042  case kChi2Poisson : break; // CDF is already CDF
2043  case kKolmogorov : {
2044  CDF = CDFKolmogorov;
2045  } break;
2046  default: break;
2047  }
2048 
2049  if(fFitControl) {
2050  // as an additional quality check, see if fitting a control fit has a higher significance
2051  _tempSwap.Fit(fFitControl, fFitModulationOptions.Data(), "", lowBound, upBound);
2052  Double_t CDFControl(-1.);
2053  switch (fFitGoodnessTest) {
2054  case kChi2ROOT : {
2055  CDFControl = 1.-ChiSquareCDF(fFitControl->GetNDF(), fFitModulation->GetChisquare());
2056  } break;
2057  case kChi2Poisson : {
2058  CDFControl = 1.-ChiSquareCDF(fFitControl->GetNDF(), ChiSquare(_tempSwap, fFitModulation));
2059  } break;
2060  case kKolmogorov : {
2061  CDFControl = KolmogorovTest(/*_tempSwap, fFitControl*/);
2062  } break;
2063  default: break;
2064  }
2065  if(CDFControl > CDF) {
2066  CDF = -1.; // control fit is more significant, so throw out the 'old' fit
2067  fHistRhoStatusCent->Fill(fCent, -1);
2068  }
2069  }
2070  if(CDF >= fMinPvalue && CDF <= fMaxPvalue && ( fFitModulation->GetMinimum(0, TMath::TwoPi()) > 0)) {
2071  // fit quality. not that although with limited acceptance the fit is performed on just
2072  // part of phase space, the requirement that energy desntiy is larger than zero is applied
2073  // to the FULL spectrum
2074  fHistRhoStatusCent->Fill(fCent, 0.);
2075  // for LOCAL didactic purposes, save the best and the worst fits
2076  // this routine can produce a lot of output histograms (it's not memory 'safe') and will not work on GRID
2077  // since the output will become unmergeable (i.e. different nodes may produce conflicting output)
2078  switch (fRunModeType) {
2079  case kLocal : {
2080  if(fRandom->Uniform(0, 100) > fPercentageOfFits) break;
2081  static Int_t didacticCounterBest(0);
2082  TProfile* didacticProfile = (TProfile*)_tempSwap.Clone(Form("Fit_%i_1-CDF_%.3f_cen_%i_%s", didacticCounterBest, CDF, fInCentralitySelection, detector.Data()));
2083  TF1* didacticFit = (TF1*)fFitModulation->Clone(Form("fit_%i_CDF_%.3f_cen_%i_%s", didacticCounterBest, CDF, fInCentralitySelection, detector.Data()));
2084  switch(fFitModulationType) {
2085  case kCombined : {
2086  // to make a nice picture also plot the separate components (v2 and v3) of the fit
2087  // only done for cobined fit where there are actually components to split ...
2088  TF1* v0(new TF1("dfit_kV2", "[0]", 0, TMath::TwoPi()));
2089  v0->SetParameter(0, didacticFit->GetParameter(0)); // normalization
2090  v0->SetLineColor(kMagenta);
2091  v0->SetLineStyle(7);
2092  didacticProfile->GetListOfFunctions()->Add(v0);
2093  TF1* v2(new TF1("dfit_kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
2094  v2->SetParameter(0, didacticFit->GetParameter(0)); // normalization
2095  v2->SetParameter(3, didacticFit->GetParameter(3)); // v2
2096  v2->FixParameter(1, 1.); // constant
2097  v2->FixParameter(2, 2.); // constant
2098  v2->FixParameter(4, didacticFit->GetParameter(4)); // psi2
2099  v2->SetLineColor(kGreen);
2100  didacticProfile->GetListOfFunctions()->Add(v2);
2101  TF1* v3(new TF1("dfit_kV3", "[0]*([1]+[2]*[3]*TMath::Cos([5]*(x-[4])))", 0, TMath::TwoPi()));
2102  v3->SetParameter(0, didacticFit->GetParameter(0)); // normalization
2103  v3->SetParameter(3, didacticFit->GetParameter(7)); // v3
2104  v3->FixParameter(1, 1.); // constant
2105  v3->FixParameter(2, 2.); // constant
2106  v3->FixParameter(4, didacticFit->GetParameter(6)); // psi3
2107  v3->FixParameter(5, 3.); // constant
2108  v3->SetLineColor(kCyan);
2109  didacticProfile->GetListOfFunctions()->Add(v3);
2110  }
2111  default : break;
2112  }
2113  didacticProfile->GetListOfFunctions()->Add(didacticFit);
2114  didacticProfile->GetYaxis()->SetTitle("#frac{d #sum #it{p}_{T}}{d #varphi} [GeV/#it{c}]");
2115  didacticProfile->GetXaxis()->SetTitle("#varphi");
2116  fOutputListGood->Add(didacticProfile);
2117  didacticCounterBest++;
2118  TH2F* didacticSurface = BookTH2F(Form("surface_%s", didacticProfile->GetName()), "#phi", "#eta", 50, 0, TMath::TwoPi(), 50, -1, 1, -1, kFALSE);
2119  for(Int_t i(0); i < iTracks; i++) {
2120  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
2121  if(PassesCuts(track)) {
2122  if(fUsePtWeight) didacticSurface->Fill(track->Phi(), track->Eta(), track->Pt());
2123  else didacticSurface->Fill(track->Phi(), track->Eta());
2124  }
2125  }
2126  if(fExcludeLeadingJetsFromFit) { // visualize the excluded region
2127  TF2 *f2 = new TF2(Form("%s_LJ", didacticSurface->GetName()),"[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])", 0, TMath::TwoPi(), -1, 1);
2128  f2->SetParameters(excludeInPt/3.,excludeInPhi,.1,excludeInEta,.1);
2129  didacticSurface->GetListOfFunctions()->Add(f2);
2130  }
2131  fOutputListGood->Add(didacticSurface);
2132  } break;
2133  default : break;
2134  }
2135  } else { // if the fit is of poor quality revert to the original rho estimate
2136  switch (fRunModeType) { // again see if we want to save the fit
2137  case kLocal : {
2138  static Int_t didacticCounterWorst(0);
2139  if(fRandom->Uniform(0, 100) > fPercentageOfFits) break;
2140  TProfile* didacticProfile = (TProfile*)_tempSwap.Clone(Form("Fit_%i_1-CDF_%.3f_cen_%i_%s", didacticCounterWorst, CDF, fInCentralitySelection, detector.Data() ));
2141  TF1* didacticFit = (TF1*)fFitModulation->Clone(Form("fit_%i_p_%.3f_cen_%i_%s", didacticCounterWorst, CDF, fInCentralitySelection, detector.Data()));
2142  didacticProfile->GetListOfFunctions()->Add(didacticFit);
2143  fOutputListBad->Add(didacticProfile);
2144  didacticCounterWorst++;
2145  } break;
2146  default : break;
2147  }
2148  switch (fFitModulationType) {
2149  case kNoFit : break; // nothing to do
2150  case kCombined : fFitModulation->SetParameter(7, 0); // no break
2151  case kFourierSeries : fFitModulation->SetParameter(7, 0); // no break
2152  default : { // needs to be done if there was a poor fit
2153  fFitModulation->SetParameter(3, 0);
2154  fFitModulation->SetParameter(0, fLocalRho->GetVal());
2155  } break;
2156  }
2157  if(CDF > -.5) fHistRhoStatusCent->Fill(fCent, 1.);
2158  return kFALSE; // return false if the fit is rejected
2159  }
2160  return kTRUE;
2161 }
2162 //_____________________________________________________________________________
2164 {
2165  // event cuts
2166  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2167  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2168  #endif
2169  switch (fCollisionType) {
2170  case kJetFlowMC : {
2172  return kTRUE;
2173  } break;
2174  case kPbPb10h : {
2175  // ugly hack for 10h data
2176  UInt_t trigger(0);
2177  AliAODEvent* aodEvent = static_cast<AliAODEvent*>(InputEvent());
2178  if(aodEvent) trigger = ((AliVAODHeader*)(aodEvent->GetHeader()))->GetOfflineTrigger();
2179  else return kFALSE;
2180  if((trigger & AliVEvent::kMB) == 0) return kFALSE;
2181  } break;
2182  default : {
2183  if(!event || !AliAnalysisTaskEmcal::IsEventSelected()) return kFALSE;
2184  } break;
2185  }
2186  // aod and esd specific checks
2187  switch (fDataType) {
2188  case kESD: {
2189  AliESDEvent* esdEvent = static_cast<AliESDEvent*>(InputEvent());
2190  if( (!esdEvent) || (TMath::Abs(esdEvent->GetPrimaryVertexSPD()->GetZ() - esdEvent->GetPrimaryVertex()->GetZ()) > .5) ) return kFALSE;
2191  } break;
2192  case kAOD: {
2193  AliAODEvent* aodEvent = static_cast<AliAODEvent*>(InputEvent());
2194  if( (!aodEvent) || (TMath::Abs(aodEvent->GetPrimaryVertexSPD()->GetZ() - aodEvent->GetPrimaryVertex()->GetZ()) > .5) ) return kFALSE;
2195  } break;
2196  default: break;
2197  }
2198  fCent = InputEvent()->GetCentrality()->GetCentralityPercentile("V0M");
2199  if(fCent <= fCentralityClasses->At(0) || fCent >= fCentralityClasses->At(fCentralityClasses->GetSize()-1) || TMath::Abs(fCent-InputEvent()->GetCentrality()->GetCentralityPercentile("TRK")) > 5.) return kFALSE;
2200  // determine centrality class
2202  for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
2203  if(fCent >= fCentralityClasses->At(i) && fCent <= fCentralityClasses->At(1+i)) {
2205  break;
2206  }
2207  }
2208  if(fInCentralitySelection < 0) return kFALSE;
2209  // see if input containers are filled
2210  if(fTracks->GetEntries() < 1) return kFALSE;
2211  if(fRho->GetVal() <= 0 ) return kFALSE;
2212  if(fAnalysisType == AliAnalysisTaskJetV2::kFull && !fClusterCont) return kFALSE;
2213  // last but not least this hideous pile-up cut for 10h data
2214  if(fCollisionType == kPbPb10h) {
2215  Float_t multTPC(0.), multGlob(0.);
2216  AliAODEvent* event = static_cast<AliAODEvent*>(InputEvent());
2217  Int_t nGoodTracks(event->GetNumberOfTracks());
2218  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
2219  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
2220  if(!trackAOD) AliFatal("Not a standard AOD");
2221  if (!trackAOD) continue;
2222  if (!(trackAOD->TestFilterBit(1))) continue;
2223  if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
2224  multTPC++;
2225  }
2226  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
2227  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
2228  if(!trackAOD) AliFatal("Not a standard AOD");
2229  if (!trackAOD) continue;
2230  if (!(trackAOD->TestFilterBit(16))) continue;
2231  if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
2232  Double_t b[2] = {-99., -99.};
2233  Double_t bCov[3] = {-99., -99., -99.};
2234  AliAODTrack copy(*trackAOD);
2235  if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
2236  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
2237  multGlob++;
2238  } //track loop
2239  if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
2240  fHistMultCorAfterCuts->Fill(multGlob, multTPC);
2241  fHistMultvsCentr->Fill(fCent, multTPC);
2242  }
2243  return kTRUE;
2244 }
2245 //_____________________________________________________________________________
2247 {
2248  // fill histograms
2249  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2250  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2251  #endif
2252  // fill histograms. weight is 1 when no procedure is defined
2256  if(fFillQAHistograms) FillWeightedEventPlaneHistograms(vzero, vzeroComb, tpc);
2259 }
2260 //_____________________________________________________________________________
2261 void AliAnalysisTaskJetV2::FillQAHistograms(AliVTrack* vtrack) const
2262 {
2263  // fill qa histograms for pico tracks
2264  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_2
2265  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2266  #endif
2267  if(!vtrack) return;
2268  AliPicoTrack* track = static_cast<AliPicoTrack*>(vtrack);
2269  fHistRunnumbersPhi->Fill(fMappedRunNumber, track->Phi());
2270  fHistRunnumbersEta->Fill(fMappedRunNumber, track->Eta());
2271  Int_t type((int)(track->GetTrackType()));
2272  switch (type) {
2273  case 0:
2274  fHistPicoCat1[fInCentralitySelection]->Fill(track->Eta(), track->Phi());
2275  break;
2276  case 1:
2277  fHistPicoCat2[fInCentralitySelection]->Fill(track->Eta(), track->Phi());
2278  break;
2279  case 2:
2280  fHistPicoCat3[fInCentralitySelection]->Fill(track->Eta(), track->Phi());
2281  break;
2282  default: break;
2283  }
2284 }
2285 //_____________________________________________________________________________
2287 {
2288  // fill qa histograms for events
2289  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2290  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2291  #endif
2292  if(!vevent) return;
2293  fHistVertexz->Fill(fVertex[2]);
2294  fHistCentrality->Fill(fCent);
2295  Int_t runNumber(InputEvent()->GetRunNumber());
2298  if(fExpectedRuns->At(fMappedRunNumber) == runNumber) return;
2299  }
2300  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_
2301  printf("\n > TASK %s CANNOT IDENTIFY RUN - CONFIGURATION COULD BE INCORRECT < \n", GetName());
2302  #endif
2303 }
2304 //_____________________________________________________________________________
2306 {
2307  // fill track histograms
2308  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2309  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2310  #endif
2311  Int_t iTracks(fTracks->GetEntriesFast()), iAcceptedTracks(0);
2312  for(Int_t i(0); i < iTracks; i++) {
2313  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
2314  if(!PassesCuts(track)) continue;
2315  iAcceptedTracks++;
2318  }
2320 }
2321 //_____________________________________________________________________________
2323 {
2324  // fill cluster histograms
2325  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2326  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2327  #endif
2328  if(!fClusterCont) return;
2329  Int_t iClusters(fClusterCont->GetNClusters());
2330  TLorentzVector clusterLorentzVector;
2331  for(Int_t i(0); i < iClusters; i++) {
2332  AliVCluster* cluster = fClusterCont->GetCluster(i);
2333  if (!PassesCuts(cluster)) continue;
2334  cluster->GetMomentum(clusterLorentzVector, const_cast<Double_t*>(fVertex));
2335  fHistClusterPt[fInCentralitySelection]->Fill(clusterLorentzVector.Pt(), fEventPlaneWeight);
2336  fHistClusterEtaPhi[fInCentralitySelection]->Fill(clusterLorentzVector.Eta(), clusterLorentzVector.Phi(), fEventPlaneWeight);
2337  fHistClusterEtaPhiWeighted[fInCentralitySelection]->Fill(clusterLorentzVector.Eta(), clusterLorentzVector.Phi(), clusterLorentzVector.Pt()*fEventPlaneWeight);
2338  }
2339  return;
2340 }
2341 //_____________________________________________________________________________
2343 {
2344  // fill event plane histograms, only called in qa mode
2345  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2346  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2347  #endif
2348  Double_t TRK(InputEvent()->GetCentrality()->GetCentralityPercentile("TRK"));
2349  Double_t V0M(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
2350  fHistPsiVZEROAV0M->Fill(V0M, vzero[0][0], fEventPlaneWeight);
2351  fHistPsiVZEROCV0M->Fill(V0M, vzero[1][0], fEventPlaneWeight);
2352  fHistPsiVZEROVV0M->Fill(V0M, vzeroComb[0], fEventPlaneWeight);
2353  fHistPsiTPCV0M->Fill(V0M, tpc[0], fEventPlaneWeight);
2354  fHistPsiVZEROATRK->Fill(TRK, vzero[0][0], fEventPlaneWeight);
2355  fHistPsiVZEROCTRK->Fill(TRK, vzero[1][0], fEventPlaneWeight);
2356  fHistPsiVZEROTRK->Fill(TRK, vzeroComb[0], fEventPlaneWeight);
2357  fHistPsiTPCTRK->Fill(TRK, tpc[0], fEventPlaneWeight);
2358  // leading jet vs event plane bias
2359  if(fLeadingJet) {
2362  Double_t pt(fLeadingJet->Pt() - fLeadingJet->Area()*rho);
2367  }
2368  // correlation of event planes
2369  fHistPsi2Correlation[fInCentralitySelection]->Fill(tpc[0], vzero[0][0], vzero[1][0], fEventPlaneWeight);
2370 }
2371 //_____________________________________________________________________________
2373 {
2374  // fill rho histograms
2375  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2376  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2377  #endif
2378  fHistRhoPackage[fInCentralitySelection]->Fill(fLocalRho->GetVal(), fEventPlaneWeight); // save the rho estimate from the emcal jet package
2379  // get multiplicity FIXME inefficient
2380  Int_t iJets(fJets->GetEntriesFast());
2381  Double_t rho(fLocalRho->GetLocalVal(TMath::Pi(), TMath::Pi(), fLocalRho->GetVal()));
2382  if(fUse2DIntegration) rho = fLocalRho->GetLocalValInEtaPhi(TMath::Pi(), TMath::Pi(), fLocalRho->GetVal());
2384  fHistRhoVsMult->Fill(fTracks->GetEntries(), rho, fEventPlaneWeight);
2385  fHistRhoVsCent->Fill(fCent, rho, fEventPlaneWeight);
2386  for(Int_t i(0); i < iJets; i++) {
2387  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
2388  if(!PassesCuts(jet)) continue;
2389  fHistRhoAVsMult->Fill(fTracks->GetEntries(), rho * jet->Area(), fEventPlaneWeight);
2390  fHistRhoAVsCent->Fill(fCent, rho * jet->Area(), fEventPlaneWeight);
2391  }
2392 }
2393 //_____________________________________________________________________________
2395 {
2396  // fill delta pt histograms
2397  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2398  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2399  #endif
2400  Int_t i(0);
2401  const Float_t areaRC = GetJetRadius()*GetJetRadius()*TMath::Pi();
2402  // we're retrieved the leading jet, now get a random cone
2403  for(i = 0; i < fMaxCones; i++) {
2404  Float_t pt(0), eta(0), phi(0);
2405  // get a random cone without constraints on leading jet position
2406  CalculateRandomCone(pt, eta, phi, fTracksCont, fClusterCont, 0x0);
2407  if(pt > 0) {
2411 /* if(fFillQAHistograms) {
2412  Double_t temp(fLocalRho->GetLocalValInEtaPhi(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal())*areaRC);
2413  fHistIntegralCorrelations[fInCentralitySelection]->Fill(fLocalRho->GetLocalVal(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal())*areaRC, temp);
2414  if(temp > 0) fProfIntegralCorrelations[fInCentralitySelection]->Fill(temp, fLocalRho->GetLocalVal(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal())*areaRC/temp);
2415  }*/
2419  fHistDeltaPtDeltaPhi2Rho0[fInCentralitySelection]->Fill(PhaseShift(phi-psi2, 2.), pt - areaRC*fLocalRho->GetVal(), fEventPlaneWeight);
2420 
2421  }
2422  // get a random cone excluding leading jet area
2424  if(pt > 0) {
2431  fHistDeltaPtDeltaPhi2ExLJRho0[fInCentralitySelection]->Fill(PhaseShift(phi-psi2, 2.), pt - areaRC*fLocalRho->GetVal(), fEventPlaneWeight);
2432  }
2433  }
2434 }
2435 //_____________________________________________________________________________
2437 {
2438  // fill jet histograms
2439  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2440  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2441  #endif
2442  Int_t iJets(fJets->GetEntriesFast());
2443  UInt_t trigger(0);
2444  if(fFillQAHistograms) {
2445  trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
2446  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2447  PrintTriggerSummary(trigger);
2448  #endif
2449  }
2450  for(Int_t i(0); i < iJets; i++) {
2451  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
2452  if(fFillQAHistograms) {
2453  if(jet) {
2454  // this is a bit redundant, but today i'm lazy
2455  Double_t pt(jet->Pt()), area(jet->Area()), eta(jet->Eta()), phi(jet->Phi());
2457  fHistRhoEtaBC[fInCentralitySelection]->Fill(rho, eta);
2458  fHistJetPtBC[fInCentralitySelection]->Fill(pt-area*rho);
2459  fHistJetEtaPhiBC[fInCentralitySelection]->Fill(eta, phi);
2460  fHistJetPtAreaBC[fInCentralitySelection]->Fill(pt-area*rho,area);
2461  }
2462  }
2463  if(PassesCuts(jet)) {
2464  Double_t pt(jet->Pt()), area(jet->Area()), eta(jet->Eta()), phi(jet->Phi());
2469  if(fFillQAHistograms) {
2471  FillWeightedTriggerQA(PhaseShift(phi-psi2, 2.), pt - area*rho, trigger);
2472  }
2473  fHistJetPtArea[fInCentralitySelection]->Fill(pt-area*rho, area, fEventPlaneWeight);
2474  fHistJetPtEta[fInCentralitySelection]->Fill(pt-area*rho, eta, fEventPlaneWeight);
2475  fHistJetPsi2Pt[fInCentralitySelection]->Fill(PhaseShift(phi-psi2, 2.), pt-area*rho, fEventPlaneWeight);
2476  AliVParticle* lp(GetLeadingTrack(jet));
2477  if(lp) {
2478  fHistJetLJPsi2Pt[fInCentralitySelection]->Fill(PhaseShift(lp->Phi()-psi2, 2.), pt-area*rho, lp->Pt(), fEventPlaneWeight);
2479  fHistJetLJPsi2PtRatio[fInCentralitySelection]->Fill(PhaseShift(lp->Phi()-psi2, 2.), PhaseShift(phi-psi2, 2.), pt-area*rho, fEventPlaneWeight);
2480  }
2481  fHistJetPsi2PtRho0[fInCentralitySelection]->Fill(PhaseShift(phi-psi2, 2.), pt-area*fLocalRho->GetVal(), fEventPlaneWeight);
2484  }
2485  }
2486 }
2487 //_____________________________________________________________________________
2489 {
2490  // fill qa histograms for pico tracks
2491  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_2
2492  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2493  #endif
2494  if(!vtrack) return;
2495  AliPicoTrack* track = static_cast<AliPicoTrack*>(vtrack);
2498  Int_t type((int)(track->GetTrackType()));
2499  switch (type) {
2500  case 0:
2501  fHistPicoCat1[fInCentralitySelection]->Fill(track->Eta(), track->Phi(), fEventPlaneWeight);
2502  break;
2503  case 1:
2504  fHistPicoCat2[fInCentralitySelection]->Fill(track->Eta(), track->Phi(), fEventPlaneWeight);
2505  break;
2506  case 2:
2507  fHistPicoCat3[fInCentralitySelection]->Fill(track->Eta(), track->Phi(), fEventPlaneWeight);
2508  break;
2509  default: break;
2510  }
2511 }
2512 //_____________________________________________________________________________
2514 {
2515  // fill qa histograms for events
2516  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2517  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2518  #endif
2519  if(!vevent) return;
2520  fHistVertexz->Fill(vevent->GetPrimaryVertex()->GetZ());
2521  fHistCentrality->Fill(fCent);
2522  Int_t runNumber(InputEvent()->GetRunNumber());
2525  if(fExpectedRuns->At(fMappedRunNumber) == runNumber) return;
2526  }
2527  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2528  printf("\n > TASK %s CANNOT IDENTIFY RUN - CONFIGURATION COULD BE INCORRECT < \n", GetName());
2529  #endif
2531  // check if cabration was kickstarted properly. this comes down to seeing if there's a difference between the
2532  // current runnumber and the runnumber as used by the calibration. if there's a difference, flag the offending
2533  // runnumber
2535  }
2536 }
2537 //_____________________________________________________________________________
2539 {
2540  // fill the trigger efficiency histograms
2541  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_2
2542  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2543  #endif
2544  // some trigger definitions for readability. the way this routine is set up is as follows
2545  // 1) define combined trigger conditions, e.g. bitwise representation of a combined trigger
2546  // trigger a = 0 0 1
2547  // trigger b = 1 0 0
2548  // combined trigger mask = 1 0 1
2549  // combined trigger is mask is defined using bitwise OR
2550  // 2) check the condition using bitwise AND and equals operator on unsigned integer
2551  // (incoming trigger & mask) == mask
2552  // 2a) which will do, when incoming trigger equals mask
2553  // 1 0 1 & 1 0 1 -> 1 0 1
2554  // when checked against requested mask
2555  // UInt_t(1 0 1) == UInt_t(1 0 1) returns true
2556  // 2b) for an imcompatible trigger, e.g.
2557  // 0 0 1 & 1 0 1 -> 0 0 1
2558  // UInt_t(0 0 1) == UInt_t(1 0 1) returns false
2559 
2560  // preparing the combined trigger masks
2561  UInt_t MB_EMCEJE(AliVEvent::kMB | AliVEvent::kEMCEJE);
2562  UInt_t CEN_EMCEJE(AliVEvent::kCentral | AliVEvent::kEMCEJE);
2563  UInt_t SEM_EMCEJE(AliVEvent::kSemiCentral | AliVEvent::kEMCEJE);
2564  UInt_t ALL_EMCEJE(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral | AliVEvent::kEMCEJE);
2565  UInt_t MB_EMCEGA(AliVEvent::kMB | AliVEvent::kEMCEGA);
2566  UInt_t CEN_EMCEGA(AliVEvent::kCentral | AliVEvent::kEMCEGA);
2567  UInt_t SEM_EMCEGA(AliVEvent::kSemiCentral | AliVEvent::kEMCEGA);
2568  UInt_t ALL_EMCEGA(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral | AliVEvent::kEMCEGA);
2569  // actual routine
2570  if(IsInPlane(dPhi)) {
2571  // in plane bookkeeping of fired triggers. not 'exclusive' so no == necessary
2572  if(trigger == 0) fHistTriggerQAIn[fInCentralitySelection]->Fill(1, pt);
2573  if(trigger & AliVEvent::kAny) fHistTriggerQAIn[fInCentralitySelection]->Fill(2, pt);
2574  if(trigger & AliVEvent::kAnyINT) fHistTriggerQAIn[fInCentralitySelection]->Fill(3, pt);
2575  if(trigger & AliVEvent::kMB) fHistTriggerQAIn[fInCentralitySelection]->Fill(4, pt);
2576  if(trigger & AliVEvent::kCentral) fHistTriggerQAIn[fInCentralitySelection]->Fill(5, pt);
2577  if(trigger & AliVEvent::kSemiCentral) fHistTriggerQAIn[fInCentralitySelection]->Fill(6, pt);
2578  if(trigger & AliVEvent::kEMCEJE) fHistTriggerQAIn[fInCentralitySelection]->Fill(7, pt);
2579  if(trigger & AliVEvent::kEMCEGA) fHistTriggerQAIn[fInCentralitySelection]->Fill(8, pt);
2580  // in plane bookkeeping of trigger combinations (for efficiency)
2581  if((trigger & MB_EMCEJE) == MB_EMCEJE) fHistTriggerQAIn[fInCentralitySelection]->Fill(9, pt);
2582  if((trigger & CEN_EMCEJE) == CEN_EMCEJE) fHistTriggerQAIn[fInCentralitySelection]->Fill(10, pt);
2583  if((trigger & SEM_EMCEJE) == SEM_EMCEJE) fHistTriggerQAIn[fInCentralitySelection]->Fill(11, pt);
2584  if((trigger & ALL_EMCEJE) == ALL_EMCEJE) fHistTriggerQAIn[fInCentralitySelection]->Fill(12, pt);
2585  if((trigger & MB_EMCEGA) == MB_EMCEGA) fHistTriggerQAIn[fInCentralitySelection]->Fill(13, pt);
2586  if((trigger & CEN_EMCEGA) == CEN_EMCEGA) fHistTriggerQAIn[fInCentralitySelection]->Fill(14, pt);
2587  if((trigger & SEM_EMCEGA) == SEM_EMCEGA) fHistTriggerQAIn[fInCentralitySelection]->Fill(15, pt);
2588  if((trigger & ALL_EMCEGA) == ALL_EMCEGA) fHistTriggerQAIn[fInCentralitySelection]->Fill(16, pt);
2589  } else {
2590  // out-of-plane bookkeeping of fired triggers. not 'exclusive' so no == necessary
2591  if(trigger == 0) fHistTriggerQAOut[fInCentralitySelection]->Fill(1, pt);
2592  if(trigger & AliVEvent::kAny) fHistTriggerQAOut[fInCentralitySelection]->Fill(2, pt);
2593  if(trigger & AliVEvent::kAnyINT) fHistTriggerQAOut[fInCentralitySelection]->Fill(3, pt);
2594  if(trigger & AliVEvent::kMB) fHistTriggerQAOut[fInCentralitySelection]->Fill(4, pt);
2595  if(trigger & AliVEvent::kCentral) fHistTriggerQAOut[fInCentralitySelection]->Fill(5, pt);
2596  if(trigger & AliVEvent::kSemiCentral) fHistTriggerQAOut[fInCentralitySelection]->Fill(6, pt);
2597  if(trigger & AliVEvent::kEMCEJE) fHistTriggerQAOut[fInCentralitySelection]->Fill(7, pt);
2598  if(trigger & AliVEvent::kEMCEGA) fHistTriggerQAOut[fInCentralitySelection]->Fill(8, pt);
2599  // out-of-plane bookkeeping of trigger combinations (for efficiency)
2600  if((trigger & MB_EMCEJE) == MB_EMCEJE) fHistTriggerQAOut[fInCentralitySelection]->Fill(9, pt);
2601  if((trigger & CEN_EMCEJE) == CEN_EMCEJE) fHistTriggerQAOut[fInCentralitySelection]->Fill(10, pt);
2602  if((trigger & SEM_EMCEJE) == SEM_EMCEJE) fHistTriggerQAOut[fInCentralitySelection]->Fill(11, pt);
2603  if((trigger & ALL_EMCEJE) == ALL_EMCEJE) fHistTriggerQAOut[fInCentralitySelection]->Fill(12, pt);
2604  if((trigger & MB_EMCEGA) == MB_EMCEGA) fHistTriggerQAOut[fInCentralitySelection]->Fill(13, pt);
2605  if((trigger & CEN_EMCEGA) == CEN_EMCEGA) fHistTriggerQAOut[fInCentralitySelection]->Fill(14, pt);
2606  if((trigger & SEM_EMCEGA) == SEM_EMCEGA) fHistTriggerQAOut[fInCentralitySelection]->Fill(15, pt);
2607  if((trigger & ALL_EMCEGA) == ALL_EMCEGA) fHistTriggerQAOut[fInCentralitySelection]->Fill(16, pt);
2608  }
2609 }
2610 //_____________________________________________________________________________
2612 {
2613  // fill the analysis summary histrogram, saves all relevant analysis settigns
2614  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2615  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2616  #endif
2617  fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fJetRadius");
2618  fHistAnalysisSummary->SetBinContent(2, GetJetContainer()->GetJetRadius());
2619  fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fJetEtaMin");
2620  fHistAnalysisSummary->SetBinContent(3, GetJetContainer()->GetJetEtaMin());
2621  fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fJetEtaMax");
2622  fHistAnalysisSummary->SetBinContent(4, GetJetContainer()->GetJetEtaMax());
2623  fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fJetPhiMin");
2624  fHistAnalysisSummary->SetBinContent(5, GetJetContainer()->GetJetPhiMin());
2625  fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fJetPhiMax");
2626  fHistAnalysisSummary->SetBinContent(6, GetJetContainer()->GetJetPhiMin());
2627  fHistAnalysisSummary->GetXaxis()->SetBinLabel(16, "fForceBeamType");
2628  fHistAnalysisSummary->SetBinContent(16, fForceBeamType);
2629  fHistAnalysisSummary->GetXaxis()->SetBinLabel(17, "fMinCent");
2630  fHistAnalysisSummary->SetBinContent(17, fMinCent);
2631  fHistAnalysisSummary->GetXaxis()->SetBinLabel(18, "fMaxCent");
2632  fHistAnalysisSummary->SetBinContent(18, fMaxCent);
2633  fHistAnalysisSummary->GetXaxis()->SetBinLabel(19, "fMinVz");
2634  fHistAnalysisSummary->SetBinContent(19, fMinVz);
2635  fHistAnalysisSummary->GetXaxis()->SetBinLabel(20, "fMaxVz");
2636  fHistAnalysisSummary->SetBinContent(20, fMaxVz);
2637  fHistAnalysisSummary->GetXaxis()->SetBinLabel(21, "fOffTrigger");
2638  fHistAnalysisSummary->SetBinContent(21, fOffTrigger);
2639  fHistAnalysisSummary->GetXaxis()->SetBinLabel(34, "fitModulationType");
2640  fHistAnalysisSummary->SetBinContent(34, (int)fFitModulationType);
2641  fHistAnalysisSummary->GetXaxis()->SetBinLabel(35, "runModeType");
2642  fHistAnalysisSummary->SetBinContent(35, (int)fRunModeType);
2643  fHistAnalysisSummary->GetXaxis()->SetBinLabel(36, "data type");
2644  fHistAnalysisSummary->SetBinContent(36, (int)fDataType);
2645  fHistAnalysisSummary->GetXaxis()->SetBinLabel(37, "iterator");
2646  fHistAnalysisSummary->SetBinContent(37, 1.);
2647  fHistAnalysisSummary->GetXaxis()->SetBinLabel(38, "fMinPvalue");
2648  fHistAnalysisSummary->SetBinContent(38, fMinPvalue);
2649  fHistAnalysisSummary->GetXaxis()->SetBinLabel(39, "fMaxPvalue");
2650  fHistAnalysisSummary->SetBinContent(39, fMaxPvalue);
2651  fHistAnalysisSummary->GetXaxis()->SetBinLabel(40, "fExcludeLeadingJetsFromFit");
2653  fHistAnalysisSummary->GetXaxis()->SetBinLabel(41, "fRebinSwapHistoOnTheFly");
2654  fHistAnalysisSummary->SetBinContent(41, (int)fRebinSwapHistoOnTheFly);
2655  fHistAnalysisSummary->GetXaxis()->SetBinLabel(42, "fUsePtWeight");
2656  fHistAnalysisSummary->SetBinContent(42, (int)fUsePtWeight);
2657  fHistAnalysisSummary->GetXaxis()->SetBinLabel(44, "fSoftTrackMinPt");
2658  fHistAnalysisSummary->SetBinContent(44, fSoftTrackMinPt);
2659  fHistAnalysisSummary->GetXaxis()->SetBinLabel(45, "fSoftTrackMaxPt");
2660  fHistAnalysisSummary->SetBinContent(45, fSoftTrackMaxPt);
2661  fHistAnalysisSummary->GetXaxis()->SetBinLabel(46, "fMaxCones");
2662  fHistAnalysisSummary->SetBinContent(46, fMaxCones);
2663  fHistAnalysisSummary->GetXaxis()->SetBinLabel(47, "used rho");
2664  fHistAnalysisSummary->GetXaxis()->SetBinLabel(48, "used small rho");
2665 }
2666 //_____________________________________________________________________________
2668 {
2669  // terminate
2670  switch (fRunModeType) {
2671  case kLocal : {
2672  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2673  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2674  #endif
2675  AliAnalysisTaskJetV2::Dump();
2676  for(Int_t i(0); i < fHistAnalysisSummary->GetXaxis()->GetNbins(); i++) printf( " > flag: %s \t content %.2f \n", fHistAnalysisSummary->GetXaxis()->GetBinLabel(1+i), fHistAnalysisSummary->GetBinContent(1+i));
2677  } break;
2678  default : break;
2679  }
2680 }
2681 //_____________________________________________________________________________
2683 {
2684  // set modulation fit
2685  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2686  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2687  #endif
2688  if (fFitModulation) delete fFitModulation;
2689  fFitModulation = fit;
2690 }
2691 //_____________________________________________________________________________
2693 {
2694  // set control fit
2695  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2696  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2697  #endif
2698  if (fFitControl) delete fFitControl;
2699  if (c) {
2700  fFitControl = new TF1("controlFit", "pol0", 0, TMath::TwoPi());
2701  } else fFitControl = 0x0;
2702 }
2703 //_____________________________________________________________________________
2705 {
2706  // INTERFACE METHOD FOR OUTPUTFILE
2707  // get the detector resolution, user has ownership of the returned histogram
2708  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2709  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2710  #endif
2711  if(!fOutputList) {
2712  printf(" > Please add fOutputList first < \n");
2713  return 0x0;
2714  }
2715  TH1F* r(0x0);
2716  (cen) ? r = new TH1F("R", "R", cen->GetSize()-1, cen->GetArray()) : r = new TH1F("R", "R", 10, 0, 10);
2717  if(!cen) r->GetXaxis()->SetTitle("number of centrality bin");
2718  r->GetYaxis()->SetTitle(Form("Resolution #Psi_{%i}", h));
2719  for(Int_t i(0); i < 10; i++) {
2720  TProfile* temp((TProfile*)fOutputList->FindObject(Form("fProfV%iResolution_%i", h, i)));
2721  if(!temp) break;
2722  Double_t a(temp->GetBinContent(3)), b(temp->GetBinContent(5)), c(temp->GetBinContent(7));
2723  Double_t d(temp->GetBinContent(9)), e(temp->GetBinContent(10)), f(temp->GetBinContent(11));
2724  Double_t _a(temp->GetBinError(3)), _b(temp->GetBinError(5)), _c(temp->GetBinError(7));
2725  Double_t _d(temp->GetBinError(9)), _e(temp->GetBinError(10)), _f(temp->GetBinError(11));
2726  Double_t error(0);
2727  if(a <= 0 || b <= 0 || c <= 0 || d <= 0 || e <= 0 || f <= 0) continue;
2728  switch (det) {
2729  case kVZEROA : {
2730  r->SetBinContent(1+i, TMath::Sqrt((a*b)/c));
2731  if(i==0) r->SetNameTitle("VZEROA resolution", "VZEROA resolution");
2732  error = TMath::Power((2.*a*TMath::Sqrt((a*b)/c))/3.,2.)*_a*_a+TMath::Power((2.*b*TMath::Sqrt((a*b)/c))/3.,2.)*_b*_b+TMath::Power(2.*c*TMath::Sqrt((a*b)/c),2.)*_c*_c;
2733  if(error > 0.) error = TMath::Sqrt(error);
2734  r->SetBinError(1+i, error);
2735  } break;
2736  case kVZEROC : {
2737  r->SetBinContent(1+i, TMath::Sqrt((a*c)/b));
2738  error = TMath::Power((2.*a*TMath::Sqrt((a*c)/b))/3.,2.)*_a*_a+TMath::Power((2.*b*TMath::Sqrt((a*c)/b)),2.)*_b*_b+TMath::Power(2.*c*TMath::Sqrt((a*c)/b)/3.,2.)*_c*_c;
2739  if(error > 0.) error = TMath::Sqrt(error);
2740  if(i==0) r->SetNameTitle("VZEROC resolution", "VZEROC resolution");
2741  r->SetBinError(1+i, error);
2742  } break;
2743  case kTPC : {
2744  r->SetBinContent(1+i, TMath::Sqrt((b*c)/a));
2745  if(i==0) r->SetNameTitle("TPC resolution", "TPC resolution");
2746  r->SetBinError(1+i, TMath::Sqrt(_a*_a+_b*_b+_c*_c));
2747  } break;
2748  case kVZEROComb : {
2749  r->SetBinContent(1+i, TMath::Sqrt((d*e)/f));
2750  if(i==0) r->SetNameTitle("VZEROComb resolution", "VZEROComb resolution");
2751  r->SetBinError(1+i, TMath::Sqrt(_d*_d+_e*_e+_f*_f));
2752  } break;
2753  default : break;
2754  }
2755  }
2756  return r;
2757 }
2758 //_____________________________________________________________________________
2760 {
2761  // INTERFACE METHOD FOR OUTPUT FILE
2762  // correct the supplied differential vn histogram v for detector resolution
2763  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2764  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2765  #endif
2766  TH1F* r(GetResolutionFromOutputFile(det, h, cen));
2767  if(!r) {
2768  printf(" > Couldn't find resolution < \n");
2769  return 0x0;
2770  }
2771  Double_t res(1./r->GetBinContent(1+r->FindBin(c)));
2772  TF1* line = new TF1("line", "pol0", 0, 200);
2773  line->SetParameter(0, res);
2774  v->Multiply(line);
2775  return v;
2776 }
2777 //_____________________________________________________________________________
2779 {
2780  // INTERFACE METHOD FOR OUTPUT FILE
2781  // correct the supplied intetrated vn histogram v for detector resolution
2782  // integrated vn must have the same centrality binning as the resolotion correction
2783  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2784  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2785  #endif
2786  TH1F* r(GetResolutionFromOutputFile(det, h, cen));
2787  v->Divide(v, r);
2788  return v;
2789 }
2790 //_____________________________________________________________________________
2791 TH1F* AliAnalysisTaskJetV2::GetDifferentialQC(TProfile* refCumulants, TProfile* diffCumlants, TArrayD* ptBins, Int_t h)
2792 {
2793  // get differential QC
2794  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2795  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2796  #endif
2797  Double_t r(refCumulants->GetBinContent(h-1)); // v2 reference flow
2798  if(r > 0) r = TMath::Sqrt(r);
2799  TH1F* qc = new TH1F(Form("QC2v%i", h), Form("QC2v%i", h), ptBins->GetSize()-1, ptBins->GetArray());
2800  Double_t a(0), b(0), c(0); // dummy variables
2801  for(Int_t i(0); i < ptBins->GetSize(); i++) {
2802  if(r > 0) {
2803  a = diffCumlants->GetBinContent(1+i);
2804  b = diffCumlants->GetBinError(1+i);
2805  c = a/r;
2806  qc->SetBinContent(1+i, c);
2807  (a <= 0 || b <= 0) ? qc->SetBinError(1+i, b) : qc->SetBinError(1+i, TMath::Sqrt(c*c*b*b/(a*a)));
2808  }
2809  }
2810  return qc;
2811 }
2812 //_____________________________________________________________________________
2814 {
2815  // necessary for calibration of 10h vzero event plane. code copied from flow package
2816  // (duplicate, but i didn't want to introduce an ulgy dependency )
2817  // this function is only called when the runnumber changes
2818  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2819  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2820  #endif
2821 
2822  // 1) check if the proper chi weights for merging vzero a and vzero c ep are present
2823  // if not, use sane defaults. centrality binning is equal to that given in the fVZEROcentralityBin snippet
2824  //
2825  // when the user wants to, set the weights to 1 (effectively disabling them)
2826  // chi values can be calculated using the static helper function
2827  // AliAnalysisTaskJetV2::CalculateEventPlaneChi(Double_t res) where res is the event plane
2828  // resolution in a given centrality bin
2829  // the resolutions that were used for these defaults are
2830  Double_t chiC2[] = {0.771423, 1.10236, 1.38116, 1.48077, 1.31964, 1.10236, 0.674622, 0.600403, 0.273865};
2831  Double_t chiA2[] = {0.582214, 0.674622, 0.832214, 0.873962, 0.832214, 0.771423, 0.637146, 0.424255, 0.257385};
2832  Double_t chiC3[] = {0.493347, 0.493347, 0.458557, 0.407166, 0.356628, 0.273865, 0.176208, 6.10352e-05, 6.10352e-05};
2833  Double_t chiA3[] = {0.356628, 0.373474, 0.356628, 0.306702, 0.24115, 0.192322, 0.127869, 6.10352e-05, 6.10352e-05};
2834 
2835  if(!fChi2A) fChi2A = new TArrayD(9, chiA2);
2836  if(!fChi2C) fChi2C = new TArrayD(9, chiC2);
2837  if(!fChi3A) fChi3A = new TArrayD(9, chiA3);
2838  if(!fChi3C) fChi3C = new TArrayD(9, chiC3);
2839 
2840  Double_t sigmaC2[] = {0.000210563,0.000554248,0.00126934,0.00138031,0.00124522,0.000948494,0.00115442,0.000626186,0.000161246};
2841  Double_t sigmaA2[] = {0.000195393,0.000509235,0.00112734,0.00121416,0.00110601,0.00086572,0.0010805,0.000579927,0.00013517};
2842  Double_t sigmaC3[] = {0.000131573,0.000317261,0.000783971,0.000885244,0.000763271,0.000542612,0.000647701,0.000524767,0};
2843  Double_t sigmaA3[] = {0.000123304,0.000293338,0.000714463,0.000798547,0.00069079,0.000503398,0.000615878,0.000489984,0};
2844 
2845  if(!fSigma2A) fSigma2A = new TArrayD(9, sigmaA2);
2846  if(!fSigma2C) fSigma2C = new TArrayD(9, sigmaC2);
2847  if(!fSigma3A) fSigma3A = new TArrayD(9, sigmaA3);
2848  if(!fSigma3C) fSigma3C = new TArrayD(9, sigmaC3);
2849 
2850  // 2) check if the database file is open, if not, open it
2851  if(!fOADB || fOADB->IsZombie()) fOADB = TFile::Open("$ALICE_PHYSICS/OADB/PWGCF/VZERO/VZEROcalibEP.root");
2852  if(fOADB->IsZombie()) {
2853  printf("OADB file $ALICE_PHYSICS/OADB/PWGCF/VZERO/VZEROcalibEP.root cannot be opened, CALIBRATION FAILED !");
2854  return;
2855  }
2856 
2857  AliOADBContainer *cont = (AliOADBContainer*) fOADB->Get("hMultV0BefCorr");
2858  if(!cont){
2859  // see if database is readable
2860  printf("OADB object hMultV0BefCorr is not available in the file\n");
2861  return;
2862  }
2863  Int_t run(fRunNumber);
2864  if(!(cont->GetObject(run))){
2865  // if the run isn't recognized fall back to a default run
2866  printf("OADB object hMultVZEROBefCorr is not available for run %i (used default run 137366)\n",run);
2867  run = 137366;
2868  }
2869  // step 3) get the proper multiplicity weights from the vzero signal
2870  fVZEROgainEqualization = ((TH2F*)cont->GetObject(run))->ProfileX();
2871  if(!fVZEROgainEqualization) {
2872  AliFatal(Form("%s: Fatal error, couldn't read fVZEROgainEqualization from OADB object < \n", GetName()));
2873  return;
2874  }
2875 
2876  TF1* fpol0 = new TF1("fpol0","pol0");
2877  fVZEROgainEqualization->Fit(fpol0, "N0", "", 0, 31);
2878  fVZEROCpol = fpol0->GetParameter(0);
2879  fVZEROgainEqualization->Fit(fpol0, "N0", "", 32, 64);
2880  fVZEROApol = fpol0->GetParameter(0);
2881 
2882  // step 4) extract the information to re-weight the q-vectors
2883  for(Int_t iside=0;iside<2;iside++){
2884  for(Int_t icoord=0;icoord<2;icoord++){
2885  for(Int_t i=0;i < 9;i++){
2886  char namecont[100];
2887  if(iside==0 && icoord==0)
2888  snprintf(namecont,100,"hQxc2_%i",i);
2889  else if(iside==1 && icoord==0)
2890  snprintf(namecont,100,"hQxa2_%i",i);
2891  else if(iside==0 && icoord==1)
2892  snprintf(namecont,100,"hQyc2_%i",i);
2893  else if(iside==1 && icoord==1)
2894  snprintf(namecont,100,"hQya2_%i",i);
2895 
2896  cont = (AliOADBContainer*) fOADB->Get(namecont);
2897  if(!cont){
2898  printf("OADB object %s is not available in the file\n",namecont);
2899  return;
2900  }
2901 
2902  if(!(cont->GetObject(run))){
2903  printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
2904  run = 137366;
2905  }
2906 
2907  // store info for all centralities to cache
2908  fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
2909  fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
2910 
2911  //for v3
2912  if(iside==0 && icoord==0)
2913  snprintf(namecont,100,"hQxc3_%i",i);
2914  else if(iside==1 && icoord==0)
2915  snprintf(namecont,100,"hQxa3_%i",i);
2916  else if(iside==0 && icoord==1)
2917  snprintf(namecont,100,"hQyc3_%i",i);
2918  else if(iside==1 && icoord==1)
2919  snprintf(namecont,100,"hQya3_%i",i);
2920 
2921  cont = (AliOADBContainer*) fOADB->Get(namecont);
2922  if(!cont){
2923  printf("OADB object %s is not available in the file\n",namecont);
2924  return;
2925  }
2926 
2927  if(!(cont->GetObject(run))){
2928  printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
2929  run = 137366;
2930  }
2931  // store info for all centralities to cache
2932  fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
2933  fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
2934  }
2935  }
2936  }
2937  // cleanup. the opened file is closed in the destructor, otherwise fVZEROgainEqualization is no longer available
2938  delete fpol0;
2939  // for qa store the runnumber that is currently used for calibration purposes
2940  fRunNumberCaliInfo = run;
2941 }
2942 //_____________________________________________________________________________i
2944 {
2945  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2946  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2947  #endif
2948 
2949  // 1) check if the proper chi weights for merging vzero a and vzero c ep are present
2950  // if not, use sane defaults. centrality binning is equal to that given in the fVZEROcentralityBin snippet
2951  //
2952  // when the user wants to, set the weights to 1 (effectively disabling them)
2953  // chi values can be calculated using the static helper function
2954  // AliAnalysisTaskJetV2::CalculateEventPlaneChi(Double_t res) where res is the event plane
2955  // resolution in a given centrality bin
2956  // the resolutions that were used for these defaults are
2957  Double_t chiC2[] = {0.771423, 1.10236, 1.38116, 1.48077, 1.31964, 1.10236, 0.674622, 0.600403, 0.273865};
2958  Double_t chiA2[] = {0.582214, 0.674622, 0.832214, 0.873962, 0.832214, 0.771423, 0.637146, 0.424255, 0.257385};
2959  Double_t chiC3[] = {0.493347, 0.493347, 0.458557, 0.407166, 0.356628, 0.273865, 0.176208, 6.10352e-05, 6.10352e-05};
2960  Double_t chiA3[] = {0.356628, 0.373474, 0.356628, 0.306702, 0.24115, 0.192322, 0.127869, 6.10352e-05, 6.10352e-05};
2961 
2962  if(!fChi2A) fChi2A = new TArrayD(9, chiA2);
2963  if(!fChi2C) fChi2C = new TArrayD(9, chiC2);
2964  if(!fChi3A) fChi3A = new TArrayD(9, chiA3);
2965  if(!fChi3C) fChi3C = new TArrayD(9, chiC3);
2966 
2967  Double_t sigmaC2[] = {7.50161e-05,0.000186685,0.000283528,0.000251427,0.000258122,2.26943e-05,0,0,0};
2968  Double_t sigmaA2[] = {0.000633027,0.000598435,0.000520023,0.000602312,0.00141679,0.00351296,0,0,0};
2969  Double_t sigmaC3[] = {4.69125e-05,0.000106922,0.000177552,0.000149093,0.000149436,0,0,0,0};
2970  Double_t sigmaA3[] = {0.000651813,0.000686852,0.000713499,0.000759663,0.00153532,0,0,0,0};
2971 
2972  if(!fSigma2A) fSigma2A = new TArrayD(9, sigmaA2);
2973  if(!fSigma2C) fSigma2C = new TArrayD(9, sigmaC2);
2974  if(!fSigma3A) fSigma3A = new TArrayD(9, sigmaA3);
2975  if(!fSigma3C) fSigma3C = new TArrayD(9, sigmaC3);
2976 }
2977 //_____________________________________________________________________________
2979 {
2980  // return cache index number corresponding to the event centrality
2981  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
2982  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2983  #endif
2984  Float_t v0Centr(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
2985  if(v0Centr < 5) return 0;
2986  else if(v0Centr < 10) return 1;
2987  else if(v0Centr < 20) return 2;
2988  else if(v0Centr < 30) return 3;
2989  else if(v0Centr < 40) return 4;
2990  else if(v0Centr < 50) return 5;
2991  else if(v0Centr < 60) return 6;
2992  else if(v0Centr < 70) return 7;
2993  else return 8;
2994 }
2995 //_____________________________________________________________________________
2997  // return pointer to the highest pt jet (before background subtraction) within acceptance
2998  // only rudimentary cuts are applied on this level, hence the implementation outside of
2999  // the framework
3000  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
3001  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
3002  #endif
3003  Int_t iJets(fJets->GetEntriesFast());
3004  Double_t pt(0);
3005  AliEmcalJet* leadingJet(0x0);
3006  if(!localRho) {
3007  for(Int_t i(0); i < iJets; i++) {
3008  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
3009  if(!PassesSimpleCuts(jet)) continue;
3010  if(jet->Pt() > pt) {
3011  leadingJet = jet;
3012  pt = leadingJet->Pt();
3013  }
3014  }
3015  return leadingJet;
3016  } else {
3017  // return leading jet after background subtraction
3018  Double_t rho(0);
3019  for(Int_t i(0); i < iJets; i++) {
3020  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
3021  if(!PassesSimpleCuts(jet)) continue;
3022  rho = localRho->GetLocalVal(jet->Phi(), GetJetContainer()->GetJetRadius(), localRho->GetVal());
3023  if(fUse2DIntegration) rho = localRho->GetLocalValInEtaPhi(jet->Phi(), GetJetContainer()->GetJetRadius(), localRho->GetVal());
3024  if((jet->Pt()-jet->Area()*rho) > pt) {
3025  leadingJet = jet;
3026  pt = (leadingJet->Pt()-jet->Area()*rho);
3027  }
3028  }
3029  return leadingJet;
3030  }
3031  return 0x0;
3032 }
3033 //_____________________________________________________________________________
3035  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
3036  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
3037  #endif
3038 
3039  // find and return the leading constituent of the jet
3040  Double_t maxPt(-1.);
3041  Int_t iTracks(jet->GetNumberOfTracks());
3042  AliVParticle* leadingTrack(0x0);
3043  for(Int_t i(0); i < iTracks; i++) {
3044  AliVParticle* vp(static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray())));
3045  if(vp && (vp->Pt() > maxPt)) {
3046  maxPt = vp->Pt();
3047  leadingTrack = vp;
3048  }
3049  }
3050  return leadingTrack;
3051 }
3052 //_____________________________________________________________________________
3054 {
3055  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
3056  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
3057  #endif
3058 
3059  // get event weights distribution from event plane distribution
3060  TH1F* temp((TH1F*)hist->Clone(Form("EP_weights_cen_%i", c)));
3061  Double_t integral(hist->Integral()/hist->GetNbinsX());
3062  // loop over bins and extract the weights
3063  for(Int_t i(0); i < hist->GetNbinsX(); i++) {
3064  temp->SetBinError(1+i, 0.); // uncertainty is irrelevant
3065  temp->SetBinContent(1+i, integral/hist->GetBinContent(1+i));
3066  }
3067  return temp;
3068 }
3069 //_____________________________________________________________________________
3071 {
3072  #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
3073  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
3074  #endif
3075 
3076  // test function to print binary representation of given trigger mask
3077  // trigger mask is represented by 32 bits (hardcoded as it is an UInt_t )
3078  TString triggerName[] = { // trigger names and their corresponding bits. some bits have multiple names
3079  "kMB", // 0
3080  "kINT7", // 1
3081  "kMUON", // 2
3082  "kHighMult", // 3
3083  "kEMC1", // 4
3084  "kCINT5", // 5
3085  "kCMUS5 kMUSPB", // 6
3086  "kMUSH7 kMUSHPB", // 7
3087  "kMUL7 kMuonLikePB", // 8
3088  "kMUU7 kMuonUnlikePB", // 9
3089  "kEMC7 kEMC8", // 10
3090  "kMUS7", // 11
3091  "kPHI1", // 12
3092  "kPHI7 kPHI8 kPHOSPb", // 13
3093  "kEMCEJE", // 14
3094  "kEMCEGA", // 15
3095  "kCentral", // 16
3096  "kSemiCentral", // 17
3097  "kDG5", // 18
3098  "kZED", // 19
3099  "kSPI7 kSPI", // 20
3100  "kINT8", // 21
3101  "kMuonSingleLowPt", // 22
3102  "kMuonSingleHighPt8", // 23
3103  "kMuonLikeLowPt8", // 24
3104  "kMuonUnlikeLowPt8", // 25
3105  "kMuonUnlikeLowPt0", // 26
3106  "kUserDefined", // 27
3107  "kTRD"}; // 28
3108  TString notTriggered = "not fired";
3109  printf(" > trigger is %u \n ", trigger);
3110 
3111  // extract which triggers have been fired exactly and print summary of bits
3112  for (Int_t i(0); i < 29; i++) printf("[bit %i]\t [%u] [%s]\n", i, (trigger & ((UInt_t)1 << i)) ? 1U : 0U, (trigger & ((UInt_t)1 << i)) ? triggerName[i].Data() : notTriggered.Data());
3113 
3114  // print accepted trigger combinations
3115  printf(" ====== accepted trigger combinations ======= \n");
3116  UInt_t MB_EMCEJE(AliVEvent::kMB | AliVEvent::kEMCEJE);
3117  UInt_t CEN_EMCEJE(AliVEvent::kCentral | AliVEvent::kEMCEJE);
3118  UInt_t SEM_EMCEJE(AliVEvent::kSemiCentral | AliVEvent::kEMCEJE);
3119  UInt_t ALL_EMCEJE(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral | AliVEvent::kEMCEJE);
3120  UInt_t MB_EMCEGA(AliVEvent::kMB | AliVEvent::kEMCEGA);
3121  UInt_t CEN_EMCEGA(AliVEvent::kCentral | AliVEvent::kEMCEGA);
3122  UInt_t SEM_EMCEGA(AliVEvent::kSemiCentral | AliVEvent::kEMCEGA);
3123  UInt_t ALL_EMCEGA(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral | AliVEvent::kEMCEGA);
3124  if(trigger == 0) printf("(trigger == 0)\n");
3125  if(trigger & AliVEvent::kAny) printf("(trigger & AliVEvent::kAny)\n");
3126  if(trigger & AliVEvent::kAnyINT) printf("(trigger & AliVEvent::kAnyINT\n");
3127  if(trigger & AliVEvent::kMB) printf("(trigger & AliVEvent::kMB)\n");
3128  if(trigger & AliVEvent::kCentral) printf("(trigger & AliVEvent::kCentral)\n");
3129  if(trigger & AliVEvent::kSemiCentral) printf("(trigger & AliVEvent::kSemiCentral)\n");
3130  if(trigger & AliVEvent::kEMCEJE) printf("(trigger & AliVEvent::kEMCEJE)\n");
3131  if(trigger & AliVEvent::kEMCEGA) printf("(trigger & AliVEvent::kEMCEGA)\n");
3132  if((trigger & MB_EMCEJE) == MB_EMCEJE) printf("(trigger & MB_EMCEJE) == MB_EMCEJE)\n");
3133  if((trigger & CEN_EMCEJE) == CEN_EMCEJE) printf("(trigger & CEN_EMCEJE) == CEN_EMCEJE)\n");
3134  if((trigger & SEM_EMCEJE) == SEM_EMCEJE) printf("(trigger & SEM_EMCEJE) == SEM_EMCEJE)\n");
3135  if((trigger & ALL_EMCEJE) == ALL_EMCEJE) printf("(trigger & ALL_EMCEJE) == ALL_EMCEJE)\n");
3136  if((trigger & MB_EMCEGA) == MB_EMCEGA) printf("(trigger & MB_EMCEGA) == MB_EMCEGA)\n");
3137  if((trigger & CEN_EMCEGA) == CEN_EMCEGA) printf("(trigger & CEN_EMCEGA) == CEN_EMCEGA)\n");
3138  if((trigger & SEM_EMCEGA) == SEM_EMCEGA) printf("(trigger & SEM_EMCEGA) == SEM_EMCEGA)\n");
3139  if((trigger & ALL_EMCEGA) == ALL_EMCEGA) printf("(trigger & ALL_EMCEGA) == ALL_EMCEGA)\n");
3140 }
3141 //_____________________________________________________________________________
3143 {
3144  // function for simple illustration of in-plane, out-of-plane method
3145 
3146  // azimuthal distribution
3147  TF1* dNdphi = new TF1("dNdphi", "1.+2.*([0]*TMath::Cos(2.*(x-[1]))+[2]*TMath::Cos(3.*(x-[3]))+[4]*TMath::Cos(4.*(x-[5])))", 0, 2*TMath::Pi());
3148 
3149  // set harmonics
3150  dNdphi->SetParameter(0, v2); // v2
3151  dNdphi->SetParameter(2, v3); // v3
3152  dNdphi->SetParameter(4, v4); // v4
3153  Double_t in(0), out(0), r(0);
3154 
3155  for(Int_t i(0); i < nEvents; i ++) {
3156  // orthogonal event planes
3157  dNdphi->SetParameter(1, gRandom->Uniform(-TMath::Pi()/2.,TMath::Pi()/2.));
3158  dNdphi->SetParameter(3, gRandom->Uniform(-TMath::Pi()/3.,TMath::Pi()/3.));
3159  dNdphi->SetParameter(5, gRandom->Uniform(-TMath::Pi()/4.,TMath::Pi()/4.));
3160 
3161  // ep loop
3162  Double_t qx(0), qy(0);
3163  for(Int_t j(0); j < 100; j++) {
3164  Double_t x = dNdphi->GetRandom(0, TMath::TwoPi());
3165  qx+=TMath::Cos(2.*x);
3166  qy+=TMath::Sin(2.*x);
3167  }
3168  Double_t ep(TMath::ATan2(qy,qx)/2.);
3169 
3170  // track loop
3171  for(Int_t j(0); j < 500; j++) {
3172  Double_t x(dNdphi->GetRandom(0, TMath::TwoPi())-ep);
3173  x = PhaseShift(x, 2);
3174  // determine which plane it is in
3175  (x > TMath::Pi()/4. && x < 3*TMath::Pi()/4.) ? out++ : in++;
3176  }
3177  r += TMath::Cos(2.*(ep-dNdphi->GetParameter(1)));
3178  }
3179 
3180  r/=100000;
3181  cout << " event plane resolution is: " << r << endl;
3182 
3183  Double_t pre = TMath::Pi()/(r*4.);
3184  Double_t ratio = pre*((in-out)/(in+out));
3185  Double_t eout = TMath::Sqrt(out);
3186  Double_t ein = TMath::Sqrt(in);
3187  Double_t error2 = (4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout;
3188  error2 = error2*pre*pre;
3189  if(error2 > 0) error2 = TMath::Sqrt(error2);
3190 
3191  dNdphi->SetTitle("total");
3192  dNdphi->DrawCopy();
3193  cout << "in: " << in << "\t out: " << out << endl;
3194  cout << "v2: " << ratio << "\t error: " << error2 << endl;
3195 
3196  TF1* dNdphi2 = new TF1("dNdphi", "1.+2.*([0]*TMath::Cos(2.*(x-[1])))", 0, 2*TMath::Pi());
3197  TF1* dNdphi3 = new TF1("dNdphi", "1.+2.*([0]*TMath::Cos(3.*(x-[1])))", 0, 2*TMath::Pi());
3198  TF1* dNdphi4 = new TF1("dNdphi", "1.+2.*([0]*TMath::Cos(4.*(x-[1])))", 0, 2*TMath::Pi());
3199 
3200  dNdphi2->SetParameter(0, dNdphi->GetParameter(0));
3201  dNdphi2->SetParameter(1, dNdphi->GetParameter(1));
3202  dNdphi2->SetLineColor(kBlue);
3203  dNdphi2->SetLineStyle(7);
3204  dNdphi2->SetTitle("v_{2}");
3205  dNdphi2->DrawCopy("same");
3206 
3207  dNdphi3->SetParameter(0, dNdphi->GetParameter(2));
3208  dNdphi3->SetParameter(1, dNdphi->GetParameter(3));
3209  dNdphi3->SetLineColor(kGreen);
3210  dNdphi3->SetLineStyle(7);
3211  dNdphi3->SetTitle("v_{3}");
3212  dNdphi3->DrawCopy("same");
3213 
3214  dNdphi4->SetParameter(0, dNdphi->GetParameter(4));
3215  dNdphi4->SetParameter(1, dNdphi->GetParameter(5));
3216  dNdphi4->SetLineColor(kMagenta);
3217  dNdphi4->SetLineStyle(7);
3218  dNdphi4->SetTitle("v_{4}");
3219  dNdphi4->DrawCopy("same");
3220 }
3221 //_____________________________________________________________________________
3222 //Float_t AliAnalysisTaskJetV2::GetCentrality() const
3223 //{
3224 // // return centrality percentile using new framework
3225 // // return -1 when something goes wrong
3226 // #ifdef ALIANALYSISTASKJETV2_DEBUG_FLAG_1
3227 // printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
3228 // #endif
3229 // AliMultSelection *multSelection = 0x0;
3230 // if(!InputEvent()) return -1.;
3231 // multSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
3232 // if(multSelection) multSelection->GetMultiplicityPercentile("V0M");
3233 // return -1.;
3234 //}
TProfile * fProfV3
resolution parameters for v2
TH2F * fHistPsiVZEROVV0M
psi 2 from vzero c
TH2F * fHistJetPtConstituents[10]
jet pt versus eta (temp control)
Double_t Area() const
Definition: AliEmcalJet.h:117
TH2F * fHistQyV0c
qx v0a before cuts
TH2F * fHistQxV0a
qx v0a before cuts
virtual AliVParticle * GetNextAcceptParticle()
TH2F * fHistQyV0a
qx v0a before cuts
double Double_t
Definition: External.C:58
fitModulationType fFitModulationType
accepted tracks for QCn
Definition: External.C:260
TH2F * BookTH2F(const char *name, const char *x, const char *y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t c=-1, Bool_t append=kTRUE)
TProfile * fProfV2Cumulant
extracted v2
Definition: External.C:236
void FillWeightedTriggerQA(Double_t dPhi, Double_t pt, UInt_t trigger)
TH1F * BookTH1F(const char *name, const char *x, Int_t bins, Double_t min, Double_t max, Int_t c=-1, Bool_t append=kTRUE)
const char * title
Definition: MakeQAPdf.C:26
TH3F * fHistPsiVZEROALeadingJet[10]
correlation tpc EP, LJ pt
TH2F * fHistRCPhiEta[10]
rho vs eta before cuts
TH1F * fHistJetPtBC[10]
pt of found jets (background subtracted)
AliJetContainer * GetJetContainer(Int_t i=0) const
TH2F * fHistTriggerQAIn[10]
eta phi emcal clusters, pt weighted
AliEmcalJet * GetLeadingJet(AliLocalRhoParameter *localRho=0x0)
UInt_t fOffTrigger
offline trigger for event selection
TH2F * fHistRhoVsRCPtExLJ[10]
random cone eta and phi, excl leading jet
TH1F * fHistUndeterminedRunQA
status of rho as function of centrality
TH2F * fHistRhoAVsCent
rho * A vs multiplicity for all jets
TH2F * fHistJetPtEta[10]
jet pt versus area before cuts
Double_t Eta() const
Definition: AliEmcalJet.h:108
Double_t fMinCent
min centrality for event selection
TProfile * fProfIntegralCorrelations[10]
correlate polar or local integral
TH2F * fHistJetPsi2PtRho0[10]
ratio of leading track v2 to jet v2
TH3F * fHistEPCorrelations[10]
psi 2 from tpc
void FillWeightedTrackHistograms() const
Double_t Phi() const
Definition: AliEmcalJet.h:104
Double_t GetLocalVal(Double_t phi, Double_t r, Double_t n) const
void CalculateEventPlaneCombinedVZERO(Double_t *comb) const
TH2F * fHistQyV0aBC
qx v0a before cuts
TH2F * fHistMultCorAfterCuts
accepted verte
Float_t fMeanQ[9][2][2]
event plane dependence of jet pt vs rho_0
void CalculateEventPlaneResolution(Double_t vzero[2][2], Double_t *vzeroComb, Double_t *tpc)
Double_t GetJetEtaMax() const
TH1F * fHistKolmogorovTest
correlation p value and reduced chi2
TH2F * fHistQxV0cBC
qx v0a before cuts
void FillHistogramsAfterSubtraction(Double_t psi2, Double_t vzero[2][2], Double_t *vzeroComb, Double_t *tpc)
TH2F * fHistPvalueCDFROOTCent
pdf value of chisquare p
TProfile * fHistCentralityPercIn
accepted centrality
TCanvas * c
Definition: TestFitELoss.C:172
TH1F * fHistRCPt[10]
rho * A vs rcpt
TH2F * fHistDeltaPtDeltaPhi2Rho0[10]
dpt vs dphi (psi2 - phi)
TH2F * fHistPsiVZEROAV0M
trigger qa out of plane
Double_t GetJetRadius() const
TF1 * fFitModulation
centrality bin
static void DoSimpleSimulation(Int_t nEvents=100000, Float_t v2=0.02, Float_t v3=0.04, Float_t v4=0.03)
TH3F * fHistJetLJPsi2PtRatio[10]
event plane dependence of jet pt and leading track pt
void QCnQnk(Int_t n, Int_t k, Double_t &reQ, Double_t &imQ)
Float_t fSoftTrackMinPt
temp cache for rho pointer
TH2F * fHistRunnumbersEta
run numbers averaged phi
static Bool_t IsInPlane(Double_t dPhi)
TRandom * gRandom
TH1F * fHistClusterPt[10]
pico tracks wo spd hit wo refit, constrained
TH2F * fHistChi2Cent
p value vs centrality
Float_t fWidthQv3[9][2][2]
recentering
TH2F * fHistPChi2Root
reduced chi2 from ROOT, centrality correlation
virtual void Exec(Option_t *)
TH1F * fHistRhoPackage[10]
geometric correlation of leading jet w/wo bkg subtraction
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:127
Container for particles within the EMCAL framework.
TH3F * fHistPsiVZEROCLeadingJet[10]
correlation vzeroa EP, LJ pt
TH2F * fHistEPCorrAvChi[10]
ep correlations
TH2F * fHistChi2ROOTCent
p value versus centrlaity from root
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:126
void FillWeightedDeltaPtHistograms(Double_t psi2) const
static Double_t PhaseShift(Double_t x)
void FillWeightedClusterHistograms() const
TH2F * fHistPicoCat3[10]
pico tracks wo spd hit w refit, constrained
AliJetContainer * fJetsCont
cluster container
Float_t fVZEROCpol
calibration info per disc
fitGoodnessTest fFitGoodnessTest
void CalculateQvectorCombinedVZERO(Double_t Q2[2], Double_t Q3[2]) const
void CalculateEventPlaneTPC(Double_t *tpc)
void FillWeightedJetHistograms(Double_t psi2)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
void FillWeightedQAHistograms(AliVTrack *vtrack) const
TString fLocalRhoName
name for local rho
TProfile * fHistCentralityPercOut
centrality versus perc in
TList * fOutputListBad
output list for local analysis
Double_t Eta() const
Definition: AliPicoTrack.h:37
static Double_t ChiSquare(TH1 &histo, TF1 *func)
TH2F * fHistPsiVZEROTRK
psi 2 from vzero c
TH2F * fHistIntegralCorrelations[10]
ep corr
TH2F * fHistPKolmogorov
KolmogorovTest value, centrality correlation.
TH2F * fHistRhoAVsMult
rho veruss centrality
int Int_t
Definition: External.C:63
TH2F * fHistRhoVsRCPt[10]
random cone eta and phi
TH1 * fVZEROgainEqualization
recentering
AliVParticle * GetLeadingTrack(AliEmcalJet *jet)
TH3F * fHistPsiTPCLeadingJet[10]
same qa lot
TH2F * fHistTriggerQAOut[10]
trigger qa in plane
Definition: External.C:204
Bool_t fAttachToEvent
is the analysis initialized?
Double_t CalculateQC2(Int_t harm)
static Double_t CalculateEventPlaneChi(Double_t res)
unsigned int UInt_t
Definition: External.C:33
TH2F * fHistQyV0cBC
qx v0a before cuts
TH2F * fHistPsiTPCTRK
psi 2 from combined vzero
float Float_t
Definition: External.C:68
Double_t fMaxVz
max vertex for event selection
void FillQAHistograms(AliVTrack *vtrack) const
Float_t fMinDisanceRCtoLJ
undetermined run QA
TH1F * fHistAnalysisSummary
output list for local analysis
void QCnDiffentialFlowVectors(TClonesArray *pois, TArrayD *ptBins, Bool_t vpart, Double_t *repn, Double_t *impn, Double_t *mp, Double_t *reqn, Double_t *imqn, Double_t *mq, Int_t n)
TH2F * fHistLeadingJetBackground[10]
correlation of event planes
TProfile * fProfV2Resolution[10]
v2 cumulant
AliRhoParameter * fRho
! event rho
Bool_t CorrectRho(Double_t psi2, Double_t psi3)
TH2F * fHistMultVsCellBC
qx v0a before cuts
AliEmcalJet * fLeadingJet
jets
static void NumericalOverlap(Double_t x1, Double_t x2, Double_t psi2, Double_t &percIn, Double_t &percOut, Double_t &percLost)
TH2F * fHistQxV0c
qx v0a before cuts
void SetTrackPhiLimits(Double_t min, Double_t max, Int_t c=0)
Bool_t PassesCuts(AliVParticle *track) const
virtual AliVParticle * GetParticle(Int_t i=-1) const
TH1F * fHistJetPtRaw[10]
dpt vs dphi, excl leading jet, rho_0
TH2F * fHistPvalueCDFCent
cdf value of chisquare p
AliRhoParameter * fCachedRho
Double_t GetLocalValInEtaPhi(Double_t phi, Double_t r, Double_t n, Int_t gran=20) const
void CalculateQvectorVZERO(Double_t Qa2[2], Double_t Qc2[2], Double_t Qa3[2], Double_t Qc3[2]) const
Int_t fNAcceptedTracksQCn
number of accepted tracks
BeamType fForceBeamType
forced beam type
TH1F * fHistPvalueCDFROOT
calibration info per runnumber
AliClusterContainer * GetClusterContainer(Int_t i=0) const
TH2F * fHistPsiVZEROCV0M
psi 2 from vzero a
TH2F * fHistPicoCat2[10]
pico tracks spd hit and refit
TH2F * fHistJetEtaRho[10]
jet pt versus number of constituents
Int_t fInCentralitySelection
mapped runnumer (for QA)
TH2F * fHistRhoVsCent
rho versus multiplicity
Int_t GetNClusters() const
Double_t fMinVz
min vertex for event selection
Double_t Phi() const
Definition: AliPicoTrack.h:33
TH2F * fHistRCPhiEtaExLJ[10]
dpt vs dphi, rho_0
static Double_t ChiSquareCDF(Int_t ndf, Double_t x)
TH1F * fHistPicoTrackPt[10]
resolution parameters for v3
Double_t fCent
!event centrality
AliLocalRhoParameter * fLocalRho
! local event rho
TH1F * fHistJetPt[10]
jet pt - no background subtraction
TH2F * fHistKolmogorovTestCent
KolmogorovTest value.
Int_t fRunNumberCaliInfo
current runnumber (for QA and jet, track selection)
Int_t fNAcceptedTracks
leading jet after background subtraction
Byte_t GetTrackType() const
Definition: AliPicoTrack.h:42
TH1F * fHistSwap
analysis summary
AliClusterContainer * fClusterCont
tracks
TH2F * fHistJetPtAreaBC[10]
jet pt versus area
TH2F * fHistDeltaPtDeltaPhi2ExLJ[10]
rcpt, excl leading jet
static TH1F * GetEventPlaneWeights(TH1F *hist, Int_t c)
TH1F * fHistPicoTrackMult[10]
pt of all charged tracks
TH1F * fHistVertexz
centrality versus perc lost
AliEmcalJet * fLeadingJetAfterSub
leading jet
AliVCluster * GetCluster(Int_t i) const
TClonesArray * fJets
! jets
TH2F * fHistEPCorrAvSigma[10]
ep corr
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
TH2F * fHistJetPtArea[10]
eta and phi correlation before cuts
TH1F * fHistRho[10]
rho as estimated by emcal jet package
TH2F * fHistDeltaPtDeltaPhi2[10]
rcpt
TH2F * fHistRhoEtaBC[10]
rho * A vs centrality for all jets
TH2F * fHistPsiTPCV0M
psi 2 from combined vzero
void SetJetPhiLimits(Float_t min, Float_t max, Int_t c=0)
Double_t Pt() const
Definition: AliEmcalJet.h:96
TH2F * fHistJetPsi2Pt[10]
jet eta versus rho
TH1F * CorrectForResolutionDiff(TH1F *v, detectorType detector, TArrayD *cen, Int_t c, Int_t h=2)
Float_t fMeanQv3[9][2][2]
recentering
TH1F * fHistEPBC
fHistMultVsCellBC
void CalculateEventPlaneVZERO(Double_t vzero[2][2]) const
TH2F * fHistClusterEtaPhi[10]
pt emcal clusters
Double_t QCnS(Int_t i, Int_t j)
TH2F * fHistJetEtaPhiBC[10]
eta and phi correlation
virtual Bool_t IsEventSelected()
TH2F * fHistPChi2
reduced chi2, centrlaity correlation
Double_t CalculateQC4(Int_t harm)
Float_t fWidthQ[9][2][2]
recentering
Float_t GetJetRadius() const
TH2F * fHistRunnumbersPhi
QA profile of centralty vs multiplicity.
TH2F * fHistPicoCat1[10]
multiplicity of accepted pico tracks
Float_t nEvents[nProd]
TClonesArray * fTracks
!tracks
TH1F * fHistPvalueCDF
correlation p value and reduced chi2
ClassImp(AliAnalysisTaskJetV2) AliAnalysisTaskJetV2
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:147
TH2F * fHistMultvsCentr
QA profile global and tpc multiplicity after outlier cut.
Double_t GetParticlePhiMax() const
Double_t KolmogorovTest() const
Double_t fVertex[3]
!event vertex
static Int_t OverlapsWithPlane(Double_t x1, Double_t x2, Double_t a, Double_t b, Double_t c, Double_t d, Double_t e, Double_t phi)
Bool_t QCnRecovery(Double_t psi2, Double_t psi3)
Float_t fVZEROApol
equalization histo
Bool_t fCreateHisto
whether or not create histograms
virtual void UserCreateOutputObjects()
void SetMakeGeneralHistograms(Bool_t g)
AliParticleContainer * fTracksCont
the actual weight of an event
Double_t GetParticlePhiMin() const
TH3F * fHistPsiVZEROCombLeadingJet[10]
correlation vzeroc EP, LJ pt
Base task in the EMCAL jet framework.
TH2F * fHistRhoStatusCent
p value vs kolmogorov value
TH2F * fHistPsiVZEROATRK
psi 2 from tpc
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
TH2F * fHistPsiVZEROCTRK
psi 2 from vzero a
Bool_t PassesSimpleCuts(AliEmcalJet *jet)
Int_t GetVZEROCentralityBin() const
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
TH1F * GetResolutionFromOutputFile(detectorType detector, Int_t h=2, TArrayD *c=0x0)
TH1I * fHistRunnumbersCaliInfo
run numbers averaged eta
const char Option_t
Definition: External.C:48
TH2F * fHistEPCorrChiSigma[10]
ep corr
TProfile * fProfV2
swap histogram
TH2F * fHistClusterEtaPhiWeighted[10]
eta phi emcal clusters
TProfile * fProfV3Cumulant
extracted v3
virtual void Terminate(Option_t *option)
void FillWeightedEventPlaneHistograms(Double_t vzero[2][2], Double_t *vzeroComb, Double_t *tpc) const
TList * fOutputListGood
output list
TProfile * fHistCentralityPercLost
centrality versus perc out
bool Bool_t
Definition: External.C:53
static void PrintTriggerSummary(UInt_t trigger)
void FillAnalysisSummaryHistogram() const
TH2F * fHistDeltaPtDeltaPhi2ExLJRho0[10]
dpt vs dphi, excl leading jet
TH3F * fHistJetLJPsi2Pt[10]
event plane dependence of jet pt
Double_t fMaxCent
max centrality for event selection
TProfile * fProfV3Resolution[10]
v3 cumulant
void CalculateRandomCone(Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer *tracksCont, AliClusterContainer *clusterCont=0x0, AliEmcalJet *jet=0x0) const
Container structure for EMCAL clusters.
Int_t fMappedRunNumber
runnumber of the cached calibration info
AliVCluster * GetNextAcceptCluster()
TH1F * fHistRCPtExLJ[10]
rho * A vs rcpt, excl leading jet
TH2F * fHistJetEtaPhi[10]
jet pt before area cut
TH3F * fHistPsi2Correlation[10]
correlation vzerocomb EP, LJ pt
TH1F * CorrectForResolutionInt(TH1F *v, detectorType detector, TArrayD *cen, Int_t h=2)
TH2F * fHistRhoVsMult
background
TH3F * BookTH3F(const char *name, const char *x, const char *y, const char *z, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t binsz, Double_t minz, Double_t maxz, Int_t c=-1, Bool_t append=kTRUE)
TH1F * GetDifferentialQC(TProfile *refCumulants, TProfile *diffCumlants, TArrayD *ptBins, Int_t h)
TArrayD * fChi2A
calibration info per disc
TH2F * fHistMultVsCell
fHistEP