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