]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGCF/Correlations/DPhi/TriggerPID/AliTwoParticlePIDCorr.cxx
Fix for TOF new calib task for CPass
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / TriggerPID / AliTwoParticlePIDCorr.cxx
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 #include "AliTwoParticlePIDCorr.h"
17 #include "AliVParticle.h"
18 #include "TFormula.h"
19 #include "TAxis.h"
20 #include "TChain.h"
21 #include "TTree.h"
22 #include "TH1F.h"
23 #include "TH2F.h"
24 #include "TH3F.h"
25 #include "TProfile.h"
26 #include "TList.h"
27 #include "TFile.h"
28 #include "TGrid.h"
29 #include "TExMap.h"
30 #include "AliCentrality.h"
31 #include "Riostream.h"
32
33 #include "AliTHn.h"    
34 #include "AliCFContainer.h"
35 #include "THn.h"
36 #include "THnSparse.h"
37 #include "TBits.h"
38 #include <TSpline.h>
39 #include <AliPID.h>
40 #include "AliESDpid.h"
41 #include "AliAODpidUtil.h"
42 #include <AliPIDResponse.h>
43 #include "AliPIDCombined.h"   
44
45 #include <AliAnalysisManager.h>
46 #include <AliInputEventHandler.h>
47 #include "AliAODInputHandler.h"
48
49 #include "AliAnalysisTaskSE.h"
50 #include "AliAnalysisManager.h"
51 #include "AliCentrality.h"
52
53 #include "AliVEvent.h"
54 #include "AliAODEvent.h"
55 #include "AliAODTrack.h"
56 #include "AliVTrack.h"
57
58 #include "THnSparse.h"
59
60 #include "AliAODMCHeader.h"
61 #include "AliAODMCParticle.h"
62 #include "AliMCEventHandler.h"
63 #include "AliMCEvent.h"
64 #include "AliMCParticle.h"
65 #include "TParticle.h"
66 #include <TDatabasePDG.h>
67 #include <TParticlePDG.h>
68
69 #include "AliGenCocktailEventHeader.h"
70 #include "AliGenEventHeader.h"
71 #include "AliCollisionGeometry.h"
72 #include "AliOADBContainer.h"
73
74 #include "AliEventPoolManager.h"
75 #include "AliAnalysisUtils.h"
76 using namespace AliPIDNameSpace;
77 using namespace std;
78
79 ClassImp(AliTwoParticlePIDCorr)
80 ClassImp(LRCParticlePID)
81
82 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
83 const char * kDetectorName[]={"ITS","TPC","TOF"} ;
84 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
85 //Source code::dphicorrelations,VnV0, TaskBFpsi, AliHelperPID, 
86
87 //________________________________________________________________________
88 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
89 :AliAnalysisTaskSE(),
90   fOutput(0),
91    fOutputList(0),
92   fList(0),
93   fCentralityMethod("V0A"),
94   fPPVsMultUtils(kFALSE),
95   fSampleType("pPb"),
96  fRequestEventPlane(kFALSE),
97   fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default)
98   trkVtx(0),
99   zvtx(0),
100   fFilterBit(768),
101   fTrackStatus(0),
102   fSharedClusterCut(-1),
103  fSharedTPCmapCut(-1),
104  fSharedfraction_Pair_cut(-1),
105   fVertextype(1),
106  skipParticlesAbove(0),
107   fzvtxcut(10.0),
108   ffilltrigassoUNID(kFALSE),
109   ffilltrigUNIDassoID(kFALSE),
110   ffilltrigIDassoUNID(kTRUE),
111   ffilltrigIDassoID(kFALSE),
112   ffilltrigIDassoIDMCTRUTH(kFALSE),
113   fMaxNofMixingTracks(50000),
114   fPtOrderMCTruth(kTRUE),
115  fPtOrderDataReco(kTRUE),
116   fWeightPerEvent(kFALSE),
117   fTriggerSpeciesSelection(kFALSE),
118   fAssociatedSpeciesSelection(kFALSE),
119  fRandomizeReactionPlane(kFALSE),
120   fTriggerSpecies(SpPion),
121   fAssociatedSpecies(SpPion),
122   fCustomBinning(""),
123   fBinningString(""),
124   fSelectHighestPtTrig(kFALSE),
125   fcontainPIDtrig(kTRUE),
126   fcontainPIDasso(kFALSE),
127   SetChargeAxis(0),
128   frejectPileUp(kFALSE),
129   fminPt(0.2),
130   fmaxPt(20.0),
131   fmineta(-0.8),
132   fmaxeta(0.8),
133   fselectprimaryTruth(kTRUE),
134   fonlyprimarydatareco(kFALSE),
135   fdcacut(kFALSE),
136   fdcacutvalue(3.0),
137   ffillhistQAReco(kFALSE),
138   ffillhistQATruth(kFALSE),
139   kTrackVariablesPair(0),
140   fminPtTrig(0),
141   fmaxPtTrig(0),
142   fminPtComboeff(2.0),
143   fmaxPtComboeff(4.0), 
144   fminPtAsso(0),
145   fmaxPtAsso(0),
146  fmincentmult(0),
147  fmaxcentmult(0), 
148  fPriHistShare(0),
149   fhistcentrality(0),
150   fEventCounter(0),
151   fEtaSpectrasso(0),
152   fphiSpectraasso(0),
153   MCtruthpt(0),
154   MCtrutheta(0),
155   MCtruthphi(0),
156   MCtruthpionpt(0),
157   MCtruthpioneta(0),
158   MCtruthpionphi(0),
159   MCtruthkaonpt(0),
160   MCtruthkaoneta(0),
161   MCtruthkaonphi(0),
162   MCtruthprotonpt(0),
163   MCtruthprotoneta(0),
164   MCtruthprotonphi(0),
165   fPioncont(0),
166   fKaoncont(0),
167   fProtoncont(0),
168   fUNIDcont(0),
169   fEventno(0),
170   fEventnobaryon(0),
171   fEventnomeson(0),
172  fhistJetTrigestimate(0),
173 fTwoTrackDistancePtdip(0x0),
174 fTwoTrackDistancePtdipmix(0x0),
175   fCentralityCorrelation(0x0),
176  fHistVZEROAGainEqualizationMap(0),
177   fHistVZEROCGainEqualizationMap(0),
178  fHistVZEROChannelGainEqualizationMap(0),
179 fCentralityWeights(0),
180  fHistCentStats(0x0),
181  fHistRefmult(0x0),
182  fHistEQVZEROvsTPCmultiplicity(0x0),
183     fHistEQVZEROAvsTPCmultiplicity(0x0),
184     fHistEQVZEROCvsTPCmultiplicity(0x0),
185     fHistVZEROCvsEQVZEROCmultiplicity(0x0),
186     fHistVZEROAvsEQVZEROAmultiplicity(0x0),
187     fHistVZEROCvsVZEROAmultiplicity(0x0),
188     fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
189     fHistVZEROSignal(0x0),
190 fHistEventPlaneTruth(0x0),
191 fHistPsiMinusPhi(0x0),
192 fEventPlanePID(0x0),
193 evplaneMC(999.),
194  fgPsi2v0a(999.),
195     fgPsi2v0c(999.),
196     fgPsi2tpc(999.),
197     fgPsi3v0a(999.),
198     fgPsi3v0c(999.),
199     fgPsi3tpc(999.),
200     fgPsi2v0aMC(999.),
201     fgPsi2v0cMC(999.),
202     fgPsi2tpcMC(999.),
203     fgPsi3v0aMC(999.),
204     fgPsi3v0cMC(999.),
205     fgPsi3tpcMC(999.),
206  gReactionPlane(999.),
207   fV2(kTRUE),
208  fV3(kFALSE),
209  fIsAfter2011(kTRUE),
210   fRun(-1),
211   fNcluster(70),
212  fEPdet("V0A"),  
213  fMultV0(NULL),
214   fV0Cpol(100),
215   fV0Apol(100),
216  fHResTPCv0A2(NULL),
217 fHResTPCv0C2(NULL),
218 fHResv0Cv0A2(NULL),
219 fHResTPCv0A3(NULL),
220 fHResTPCv0C3(NULL),
221 fHResv0Cv0A3(NULL),
222  fHResMA2(NULL),
223 fHResMC2(NULL),
224 fHResAC2(NULL),
225 fHResMA3(NULL),
226 fHResMC3(NULL),
227 fHResAC3(NULL),
228 fPhiRPTPC(NULL),
229 fPhiRPTPCv3(NULL),
230 fPhiRPv0A(NULL),
231 fPhiRPv0C(NULL),
232 fPhiRPv0Av3(NULL),
233 fPhiRPv0Cv3(NULL),
234  fControlConvResoncances(0),
235   fHistoTPCdEdx(0x0),
236   fHistoTOFbeta(0x0),
237   fTPCTOFPion3d(0),
238   fTPCTOFKaon3d(0),
239   fTPCTOFProton3d(0),
240   fPionPt(0),
241   fPionEta(0),
242   fPionPhi(0),
243   fKaonPt(0),
244   fKaonEta(0),
245   fKaonPhi(0),
246   fProtonPt(0),
247   fProtonEta(0),
248   fProtonPhi(0),
249   fCorrelatonTruthPrimary(0),
250   fCorrelatonTruthPrimarymix(0),
251   fTHnCorrUNID(0),
252   fTHnCorrUNIDmix(0),
253   fTHnCorrID(0),
254   fTHnCorrIDmix(0),
255   fTHnCorrIDUNID(0),
256   fTHnCorrIDUNIDmix(0),
257   fTHnTrigcount(0),
258   fTHnTrigcountMCTruthPrim(0),
259   fPoolMgr(0x0),
260   fArrayMC(0),
261   fAnalysisType("AOD"), 
262   fefffilename(""),
263  ftwoTrackEfficiencyCutDataReco(kTRUE),
264 fTwoTrackCutMinRadius(0.8),
265 fTwoTrackCutMaxRadius(2.5),
266   twoTrackEfficiencyCutValue(0.02),
267   fPID(NULL),
268  fPIDCombined(NULL),
269  eventno(0),
270   fPtTOFPIDmin(0.5),
271   fPtTOFPIDmax(4.0),
272   fRequestTOFPID(kTRUE),
273   fPIDType(NSigmaTPCTOF),
274  fFIllPIDQAHistos(kTRUE),
275   fNSigmaPID(3),
276   fBayesCut(0.8),
277  fdiffPIDcutvalues(kFALSE),
278  fPIDCutval1(0.0),
279  fPIDCutval2(0.0),
280  fPIDCutval3(0.0),
281  fPIDCutval4(0.0),
282  fHighPtKaonNSigmaPID(-1),
283  fHighPtKaonSigma(3.5),
284   fUseExclusiveNSigma(kFALSE),
285   fRemoveTracksT0Fill(kFALSE),
286 fSelectCharge(0),
287 fTriggerSelectCharge(0),
288 fAssociatedSelectCharge(0),
289 fTriggerRestrictEta(-1),
290 fEtaOrdering(kFALSE),
291 fCutConversions(kFALSE),
292 fCutResonances(kFALSE),
293 fRejectResonanceDaughters(-1),
294   fOnlyOneEtaSide(0),
295 fInjectedSignals(kFALSE),
296   fRemoveWeakDecays(kFALSE),
297 fRemoveDuplicates(kFALSE),
298   fapplyTrigefficiency(kFALSE),
299   fapplyAssoefficiency(kFALSE),
300   ffillefficiency(kFALSE),
301   fmesoneffrequired(kFALSE),
302   fkaonprotoneffrequired(kFALSE),
303 fAnalysisUtils(0x0),
304   fDCAXYCut(0)     
305
306
307  for ( Int_t i = 0; i < 16; i++) { 
308     fHistQA[i] = NULL;
309   }
310
311  for ( Int_t i = 0; i < 6; i++ ){
312     fTrackHistEfficiency[i] = NULL;
313     effcorection[i]=NULL;
314     //effmap[i]=NULL;
315   }
316  for ( Int_t i = 0; i < 2; i++ ){
317    fTwoTrackDistancePt[i]=NULL;
318    fTwoTrackDistancePtmix[i]=NULL;
319 }
320
321  for(Int_t ipart=0;ipart<NSpecies;ipart++)
322     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
323       fnsigmas[ipart][ipid]=999.;
324
325  for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
326
327   for(Int_t i = 0; i != 2; ++i)
328     for(Int_t j = 0; j != 2; ++j)
329       for(Int_t iC = 0; iC < 9; iC++){
330         fMeanQ[iC][i][j] = 0.;
331         fWidthQ[iC][i][j] = 1.;
332         fMeanQv3[iC][i][j] = 0.;
333         fWidthQv3[iC][i][j] = 1.;
334     }
335
336   }
337 //________________________________________________________________________
338 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
339   :AliAnalysisTaskSE(name),
340  fOutput(0),
341    fOutputList(0),
342    fList(0),
343  fCentralityMethod("V0A"),
344   fPPVsMultUtils(kFALSE),
345   fSampleType("pPb"),
346  fRequestEventPlane(kFALSE),
347   fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default)
348   trkVtx(0),
349   zvtx(0),
350   fFilterBit(768),
351   fTrackStatus(0),
352   fSharedClusterCut(-1),
353   fSharedTPCmapCut(-1),
354  fSharedfraction_Pair_cut(-1),
355   fVertextype(1),
356    skipParticlesAbove(0),
357   fzvtxcut(10.0),
358   ffilltrigassoUNID(kFALSE),
359   ffilltrigUNIDassoID(kFALSE),
360   ffilltrigIDassoUNID(kTRUE),
361   ffilltrigIDassoID(kFALSE),
362   ffilltrigIDassoIDMCTRUTH(kFALSE),
363   fMaxNofMixingTracks(50000),
364   fPtOrderMCTruth(kTRUE),
365   fPtOrderDataReco(kTRUE),
366   fWeightPerEvent(kFALSE),
367   fTriggerSpeciesSelection(kFALSE),
368   fAssociatedSpeciesSelection(kFALSE),
369    fRandomizeReactionPlane(kFALSE),
370   fTriggerSpecies(SpPion),
371   fAssociatedSpecies(SpPion),
372   fCustomBinning(""),
373   fBinningString(""),
374   fSelectHighestPtTrig(kFALSE),
375   fcontainPIDtrig(kTRUE),
376   fcontainPIDasso(kFALSE),
377   SetChargeAxis(0),
378   frejectPileUp(kFALSE),
379   fminPt(0.2),
380   fmaxPt(20.0),
381   fmineta(-0.8),
382   fmaxeta(0.8),
383   fselectprimaryTruth(kTRUE),
384   fonlyprimarydatareco(kFALSE),
385    fdcacut(kFALSE),
386   fdcacutvalue(3.0),
387   ffillhistQAReco(kFALSE),
388   ffillhistQATruth(kFALSE),
389  kTrackVariablesPair(0),
390   fminPtTrig(0),
391   fmaxPtTrig(0),
392   fminPtComboeff(2.0),
393   fmaxPtComboeff(4.0), 
394   fminPtAsso(0),
395   fmaxPtAsso(0),
396    fmincentmult(0),
397    fmaxcentmult(0),
398    fPriHistShare(0),
399   fhistcentrality(0),
400   fEventCounter(0),
401   fEtaSpectrasso(0),
402   fphiSpectraasso(0),
403   MCtruthpt(0),
404   MCtrutheta(0),
405   MCtruthphi(0),
406   MCtruthpionpt(0),
407   MCtruthpioneta(0),
408   MCtruthpionphi(0),
409   MCtruthkaonpt(0),
410   MCtruthkaoneta(0),
411   MCtruthkaonphi(0),
412   MCtruthprotonpt(0),
413   MCtruthprotoneta(0),
414   MCtruthprotonphi(0),
415   fPioncont(0),
416   fKaoncont(0),
417   fProtoncont(0),
418    fUNIDcont(0),
419   fEventno(0),
420   fEventnobaryon(0),
421   fEventnomeson(0),
422   fhistJetTrigestimate(0),
423 fTwoTrackDistancePtdip(0x0),
424 fTwoTrackDistancePtdipmix(0x0),
425   fCentralityCorrelation(0x0),
426  fHistVZEROAGainEqualizationMap(0),
427   fHistVZEROCGainEqualizationMap(0),
428    fHistVZEROChannelGainEqualizationMap(0),
429 fCentralityWeights(0),
430   fHistCentStats(0x0),
431   fHistRefmult(0x0),
432     fHistEQVZEROvsTPCmultiplicity(0x0),
433     fHistEQVZEROAvsTPCmultiplicity(0x0),
434     fHistEQVZEROCvsTPCmultiplicity(0x0),
435     fHistVZEROCvsEQVZEROCmultiplicity(0x0),
436     fHistVZEROAvsEQVZEROAmultiplicity(0x0),
437     fHistVZEROCvsVZEROAmultiplicity(0x0),
438     fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
439     fHistVZEROSignal(0x0),
440 fHistEventPlaneTruth(0x0),
441    fHistPsiMinusPhi(0x0),
442 fEventPlanePID(0x0),
443 evplaneMC(999.),
444  fgPsi2v0a(999.),
445     fgPsi2v0c(999.),
446     fgPsi2tpc(999.),
447     fgPsi3v0a(999.),
448     fgPsi3v0c(999.),
449     fgPsi3tpc(999.),
450     fgPsi2v0aMC(999.),
451     fgPsi2v0cMC(999.),
452     fgPsi2tpcMC(999.),
453     fgPsi3v0aMC(999.),
454     fgPsi3v0cMC(999.),
455     fgPsi3tpcMC(999.),
456    gReactionPlane(999.),
457  fV2(kTRUE),
458  fV3(kFALSE),
459  fIsAfter2011(kTRUE),
460   fRun(-1),
461   fNcluster(70),
462    fEPdet("V0A"),  
463  fMultV0(NULL),
464   fV0Cpol(100),
465   fV0Apol(100),
466  fHResTPCv0A2(NULL),
467 fHResTPCv0C2(NULL),
468 fHResv0Cv0A2(NULL),
469 fHResTPCv0A3(NULL),
470 fHResTPCv0C3(NULL),
471 fHResv0Cv0A3(NULL),
472  fHResMA2(NULL),
473 fHResMC2(NULL),
474 fHResAC2(NULL),
475 fHResMA3(NULL),
476 fHResMC3(NULL),
477 fHResAC3(NULL),
478 fPhiRPTPC(NULL),
479 fPhiRPTPCv3(NULL),
480 fPhiRPv0A(NULL),
481 fPhiRPv0C(NULL),
482 fPhiRPv0Av3(NULL),
483 fPhiRPv0Cv3(NULL),
484   fControlConvResoncances(0), 
485   fHistoTPCdEdx(0x0),
486   fHistoTOFbeta(0x0),
487   fTPCTOFPion3d(0),
488   fTPCTOFKaon3d(0),
489   fTPCTOFProton3d(0),
490   fPionPt(0),
491   fPionEta(0),
492   fPionPhi(0),
493   fKaonPt(0),
494   fKaonEta(0),
495   fKaonPhi(0),
496   fProtonPt(0),
497   fProtonEta(0),
498   fProtonPhi(0),
499   fCorrelatonTruthPrimary(0),
500  fCorrelatonTruthPrimarymix(0),
501   fTHnCorrUNID(0),
502   fTHnCorrUNIDmix(0),
503   fTHnCorrID(0),
504   fTHnCorrIDmix(0),
505   fTHnCorrIDUNID(0),
506   fTHnCorrIDUNIDmix(0),
507   fTHnTrigcount(0),
508   fTHnTrigcountMCTruthPrim(0),
509   fPoolMgr(0x0),
510   fArrayMC(0),
511   fAnalysisType("AOD"),
512   fefffilename(""),
513   ftwoTrackEfficiencyCutDataReco(kTRUE),
514 fTwoTrackCutMinRadius(0.8),
515 fTwoTrackCutMaxRadius(2.5),
516   twoTrackEfficiencyCutValue(0.02),
517   fPID(NULL),
518   fPIDCombined(NULL),
519   eventno(0),
520  fPtTOFPIDmin(0.5),
521   fPtTOFPIDmax(4.0),
522   fRequestTOFPID(kTRUE),
523   fPIDType(NSigmaTPCTOF),
524   fFIllPIDQAHistos(kTRUE),
525   fNSigmaPID(3),
526   fBayesCut(0.8),
527  fdiffPIDcutvalues(kFALSE),
528  fPIDCutval1(0.0),
529  fPIDCutval2(0.0),
530  fPIDCutval3(0.0),
531  fPIDCutval4(0.0),
532 fHighPtKaonNSigmaPID(-1),
533  fHighPtKaonSigma(3.5),
534   fUseExclusiveNSigma(kFALSE),
535   fRemoveTracksT0Fill(kFALSE),
536 fSelectCharge(0),
537 fTriggerSelectCharge(0),
538 fAssociatedSelectCharge(0),
539 fTriggerRestrictEta(-1),
540 fEtaOrdering(kFALSE),
541 fCutConversions(kFALSE),
542 fCutResonances(kFALSE),
543 fRejectResonanceDaughters(-1),
544   fOnlyOneEtaSide(0),
545 fInjectedSignals(kFALSE),
546   fRemoveWeakDecays(kFALSE),
547 fRemoveDuplicates(kFALSE),
548   fapplyTrigefficiency(kFALSE),
549   fapplyAssoefficiency(kFALSE),
550   ffillefficiency(kFALSE),
551  fmesoneffrequired(kFALSE),
552  fkaonprotoneffrequired(kFALSE),
553    fAnalysisUtils(0x0),
554   fDCAXYCut(0)    
555   // The last in the above list should not have a comma after it
556      
557 {
558   
559    for ( Int_t i = 0; i < 16; i++) { 
560     fHistQA[i] = NULL;
561   }
562  
563 for ( Int_t i = 0; i < 6; i++ ){
564     fTrackHistEfficiency[i] = NULL;
565     effcorection[i]=NULL;
566     //effmap[i]=NULL;
567   }
568
569 for ( Int_t i = 0; i < 2; i++ ){
570    fTwoTrackDistancePt[i]=NULL;
571    fTwoTrackDistancePtmix[i]=NULL;
572 }
573
574  for(Int_t ipart=0;ipart<NSpecies;ipart++)
575     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
576       fnsigmas[ipart][ipid]=999.;
577
578    for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
579
580   for(Int_t i = 0; i != 2; ++i)
581     for(Int_t j = 0; j != 2; ++j)
582       for(Int_t iC = 0; iC < 9; iC++){
583         fMeanQ[iC][i][j] = 0.;
584         fWidthQ[iC][i][j] = 1.;
585         fMeanQv3[iC][i][j] = 0.;
586         fWidthQv3[iC][i][j] = 1.;
587     }
588   
589   // Constructor
590   // Define input and output slots here (never in the dummy constructor)
591   // Input slot #0 works with a TChain - it is connected to the default input container
592   // Output slot #1 writes into a TH1 container
593  
594   DefineOutput(1, TList::Class());                                        // for output list
595   DefineOutput(2, TList::Class());
596
597 }
598
599 //________________________________________________________________________
600 AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
601 {
602   // Destructor. Clean-up the output list, but not the histograms that are put inside
603   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
604   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
605     delete fOutput;
606
607   }
608
609 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
610     delete fOutputList;
611
612   }
613
614   if (fPID) delete fPID;
615   if (fPIDCombined) delete fPIDCombined;
616
617   }
618 //________________________________________________________________________
619
620 //////////////////////////////////////////////////////////////////////////////////////////////////
621
622 TH2F* AliTwoParticlePIDCorr::GetHistogram2D(const char * name){
623   // returns histo named name
624   return (TH2F*) fOutputList->FindObject(name);
625 }
626
627 //////////////////////////////////////////////////////////////////////////////////////////////////
628
629 Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
630
631 {
632         //
633         // Puts the argument in the range [-pi/2,3 pi/2].
634         //
635         
636         if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
637         if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();      
638
639         return DPhi;
640         
641 }
642 //________________________________________________________________________
643 void AliTwoParticlePIDCorr::UserCreateOutputObjects()
644 {
645   // Create histograms
646   // Called once (on the worker node)
647   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
648   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
649   fPID = inputHandler->GetPIDResponse();
650
651   //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
652
653 //get the efficiency correction map
654
655 // global switch disabling the reference 
656   // (to avoid "Replacing existing TH1" if several wagons are created in train)
657   Bool_t oldStatus = TH1::AddDirectoryStatus();
658   TH1::AddDirectory(kFALSE);
659
660   const Int_t nPsiTOF = 10;  
661   const Int_t nCentrBin = 9;  
662
663
664   fOutput = new TList();
665   fOutput->SetOwner();  // IMPORTANT!  
666
667   fOutputList = new TList;
668   fOutputList->SetOwner();
669   fOutputList->SetName("PIDQAList");
670
671   fList = new TList;
672   fList->SetOwner();
673   fList->SetName("EPQAList");
674   
675   fEventCounter = new TH1F("fEventCounter","EventCounter", 19, 0.5,19.5);
676   fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
677   fEventCounter->GetXaxis()->SetBinLabel(3,"After PileUP Cut");//only for Data
678   fEventCounter->GetXaxis()->SetBinLabel(5,"Have A Vertex");
679   fEventCounter->GetXaxis()->SetBinLabel(7,"After vertex Cut");
680   fEventCounter->GetXaxis()->SetBinLabel(9,"Getting centrality");
681   fEventCounter->GetXaxis()->SetBinLabel(11,"After centrality flattening");
682   fEventCounter->GetXaxis()->SetBinLabel(13,"Within 0-100% centrality");
683   fEventCounter->GetXaxis()->SetBinLabel(15,"Event Analyzed");
684   //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
685   fOutput->Add(fEventCounter);
686   
687 fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
688 fOutput->Add(fEtaSpectrasso);
689
690 fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
691 fOutput->Add(fphiSpectraasso);
692
693  if(fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
694       fOutput->Add(fCentralityCorrelation);
695  }
696
697 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")
698   {
699  TString gCentName[8] = {"V0A","V0C","V0M","V0AEq","V0CEq","V0MEq","CL1","ZNA"};
700   fHistCentStats = new TH2F("fHistCentStats",
701                              "Centrality statistics;;Cent percentile",
702                             8,-0.5,7.5,220,-5,105);
703   for(Int_t i = 1; i <= 8; i++){
704     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
705     //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
706   }
707   fOutput->Add(fHistCentStats);
708   }
709
710 if(fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
711   {
712 fhistcentrality=new TH1F("fhistcentrality","referencemultiplicity",30001,-0.5,30000.5);
713 fOutput->Add(fhistcentrality);
714   }
715  else{
716 fhistcentrality=new TH1F("fhistcentrality","centrality",220,-5,105);
717 fOutput->Add(fhistcentrality);
718  }
719
720 if(fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
721   {
722 TString gmultName[4] = {"V0A_MANUAL","V0C_MANUAL","V0M_MANUAL","TRACKS_MANUAL"};
723   fHistRefmult = new TH2F("fHistRefmult",
724                              "Reference multiplicity",
725                             4,-0.5,3.5,10000,0,20000);
726   for(Int_t i = 1; i <= 4; i++){
727     fHistRefmult->GetXaxis()->SetBinLabel(i,gmultName[i-1].Data());
728     //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
729   }
730   fOutput->Add(fHistRefmult);
731
732  if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
733  //TPC vs EQVZERO multiplicity
734     fHistEQVZEROvsTPCmultiplicity = new TH2F("fHistEQVZEROvsTPCmultiplicity","EqVZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
735     fHistEQVZEROvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO multiplicity (a.u.)");
736     fHistEQVZEROvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
737     fOutput->Add(fHistEQVZEROvsTPCmultiplicity);
738
739
740     fHistEQVZEROAvsTPCmultiplicity = new TH2F("fHistEQVZEROAvsTPCmultiplicity","EqVZERO_A vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
741     fHistEQVZEROAvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
742     fHistEQVZEROAvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
743     fOutput->Add(fHistEQVZEROAvsTPCmultiplicity);
744
745
746     fHistEQVZEROCvsTPCmultiplicity = new TH2F("fHistEQVZEROCvsTPCmultiplicity","EqVZERO_C vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
747     fHistEQVZEROCvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
748     fHistEQVZEROCvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
749     fOutput->Add(fHistEQVZEROCvsTPCmultiplicity);
750
751  //EQVZERO vs VZERO multiplicity
752   fHistVZEROCvsEQVZEROCmultiplicity = new TH2F("fHistVZEROCvsEQVZEROCmultiplicity","EqVZERO_C vs VZERO_C multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
753     fHistVZEROCvsEQVZEROCmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
754     fHistVZEROCvsEQVZEROCmultiplicity->GetYaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
755     fOutput->Add(fHistVZEROCvsEQVZEROCmultiplicity);
756
757
758 fHistVZEROAvsEQVZEROAmultiplicity = new TH2F("fHistVZEROAvsEQVZEROAmultiplicity","EqVZERO_A vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
759     fHistVZEROAvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
760     fHistVZEROAvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
761     fOutput->Add(fHistVZEROAvsEQVZEROAmultiplicity);
762
763
764   //VZEROC vs VZEROA multiplicity
765 fHistVZEROCvsVZEROAmultiplicity = new TH2F("fHistVZEROCvsVZEROAmultiplicity","VZERO_C vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
766     fHistVZEROCvsVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
767     fHistVZEROCvsVZEROAmultiplicity->GetYaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
768     fOutput->Add(fHistVZEROCvsVZEROAmultiplicity);
769
770
771
772   //EQVZEROC vs EQVZEROA multiplicity
773 fHistEQVZEROCvsEQVZEROAmultiplicity = new TH2F("fHistEQVZEROCvsEQVZEROAmultiplicity","EqVZERO_C vs EqVZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
774     fHistEQVZEROCvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
775     fHistEQVZEROCvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
776     fOutput->Add(fHistEQVZEROCvsEQVZEROAmultiplicity);
777
778  fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
779   fOutput->Add(fHistVZEROSignal);
780  }
781 }
782
783  if(fRequestEventPlane){
784 //Event plane
785  
786   fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
787   fList->Add(fHistPsiMinusPhi);
788
789   fEventPlanePID = new TH3F("fEventPlanePID",";centrality;eventplane;PID",20,0.0,100.0,4,-0.5,3.5,4,-0.5,3.5);
790   fList->Add(fEventPlanePID);
791
792  }
793  
794 if(fCutConversions || fCutResonances)
795     {
796 fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
797  fOutput->Add(fControlConvResoncances);
798     }
799
800 fHistoTPCdEdx = new TH2F("fHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
801 fOutputList->Add(fHistoTPCdEdx);
802 fHistoTOFbeta = new TH2F(Form("fHistoTOFbeta"), ";p_{T} (GeV/c);v/c",100, 0., fmaxPt, 500, 0.1, 1.1);
803   fOutputList->Add(fHistoTOFbeta);
804   
805    fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
806    fOutputList->Add(fTPCTOFPion3d);
807   
808    fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
809    fOutputList->Add(fTPCTOFKaon3d);
810
811    fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
812    fOutputList->Add(fTPCTOFProton3d);
813
814 if(ffillhistQAReco)
815     {
816     fPionPt = new TH1F("fPionPt","p_{T} distribution",200,0.,10.);
817  fOutputList->Add(fPionPt);
818     fPionEta= new TH1F("fPionEta","#eta distribution",360,-1.8,1.8);
819  fOutputList->Add(fPionEta);
820     fPionPhi = new TH1F("fPionPhi","#phi distribution",340,0,6.8);
821  fOutputList->Add(fPionPhi);
822   
823     fKaonPt = new TH1F("fKaonPt","p_{T} distribution",200,0.,10.);
824  fOutputList->Add(fKaonPt);
825     fKaonEta= new TH1F("fKaonEta","#eta distribution",360,-1.8,1.8);
826  fOutputList->Add(fKaonEta);
827     fKaonPhi = new TH1F("fKaonPhi","#phi distribution",340,0,6.8);
828  fOutputList->Add(fKaonPhi);
829   
830     fProtonPt = new TH1F("fProtonPt","p_{T} distribution",200,0.,10.);
831  fOutputList->Add(fProtonPt);
832     fProtonEta= new TH1F("fProtonEta","#eta distribution",360,-1.8,1.8);
833  fOutputList->Add(fProtonEta);
834     fProtonPhi= new TH1F("fProtonPhi","#phi distribution",340,0,6.8);
835  fOutputList->Add(fProtonPhi);
836     }
837
838   fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
839   fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
840   fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);  
841   fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx  After Cut ", 50, -5., 5.);
842   fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
843   fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
844   fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
845   fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);   
846   fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
847   fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
848   fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
849   fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
850   fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
851  fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
852  fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
853  fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
854     
855 for(Int_t i = 0; i < 16; i++)
856     {
857       fOutput->Add(fHistQA[i]);
858     }
859
860     fPriHistShare = new TH1F ("fPriHistShare","Shared clusters, primaries;#shared clusters;counts",160,0,160);
861     fOutput->Add(fPriHistShare);
862
863    Int_t eventplaneaxis=0;
864
865    if (fRequestEventPlane) eventplaneaxis=1;
866
867    kTrackVariablesPair=6+SetChargeAxis+eventplaneaxis;
868
869    if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
870  
871  if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
872  
873  if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8+SetChargeAxis+eventplaneaxis;
874  
875  
876 // two particle histograms
877   Int_t anaSteps   = 1;       // analysis steps
878   const char* title = "d^{2}N_{ch}/d#varphid#eta";
879
880   Int_t iBinPair[kTrackVariablesPair];         // binning for track variables
881   Double_t* dBinsPair[kTrackVariablesPair];    // bins for track variables  
882   TString* axisTitlePair;  // axis titles for track variables
883   axisTitlePair=new TString[kTrackVariablesPair];
884
885  TString defaultBinningStr;
886   defaultBinningStr =   "eta: -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0\n"
887     "p_t_assoc: 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 8.0,10.0\n"
888     "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
889     "p_t_eff:0.0,0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0,5.5, 6.0, 7.0, 8.0,9.0,10.0\n"
890     "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
891   "delta_phi: -1.570796, -1.483530, -1.396263, -1.308997, -1.221730, -1.134464, -1.047198, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.087266, 0.0, 0.087266, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.047198, 1.134464, 1.221730, 1.308997, 1.396263, 1.483530, 1.570796, 1.658063, 1.745329, 1.832596, 1.919862, 2.007129, 2.094395, 2.181662, 2.268928, 2.356194, 2.443461, 2.530727, 2.617994, 2.705260, 2.792527, 2.879793, 2.967060, 3.054326, 3.141593, 3.228859, 3.316126, 3.403392, 3.490659, 3.577925, 3.665191, 3.752458, 3.839724, 3.926991, 4.014257, 4.101524, 4.188790, 4.276057, 4.363323, 4.450590, 4.537856, 4.625123, 4.712389\n" // this binning starts at -pi/2 and is modulo 3 
892         "delta_eta: -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,2.1, 2.2, 2.3, 2.4\n"
893       "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
894
895  if(fRequestEventPlane){
896    defaultBinningStr += "eventPlane: -0.5,0.5,1.5,2.5,3.5\n"; // Event Plane Bins (Psi: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest))
897   }
898
899   if(fcontainPIDtrig){
900       defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
901   }
902   if(fcontainPIDasso){
903       defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
904   }
905  
906   if(SetChargeAxis==2){
907       defaultBinningStr += "TrigCharge: -2.0,0.0,2.0\n"; // course
908       defaultBinningStr += "AssoCharge: -2.0,0.0,2.0\n"; // course
909   }
910  // =========================================================
911   // Customization (adopted from AliUEHistograms)
912   // =========================================================
913
914   TObjArray* lines = defaultBinningStr.Tokenize("\n");
915   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
916   {
917     TString line(lines->At(i)->GetName());
918     TString tag = line(0, line.Index(":")+1);
919     if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
920       fBinningString += line + "\n";
921     else
922       AliInfo(Form("Using custom binning for %s", tag.Data()));
923   }
924   delete lines;
925   fBinningString += fCustomBinning;
926   
927   AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
928
929  //  =========================================================
930   // Now set the bins
931   // =========================================================
932
933     dBinsPair[0]       = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
934     axisTitlePair[0]   = " multiplicity";
935
936     dBinsPair[1]     = GetBinning(fBinningString, "vertex", iBinPair[1]);
937     axisTitlePair[1]  = "v_{Z} (cm)"; 
938
939     dBinsPair[2]     = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
940     axisTitlePair[2]    = "p_{T,trig.} (GeV/c)"; 
941
942     dBinsPair[3]     = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
943     axisTitlePair[3]    = "p_{T,assoc.} (GeV/c)";
944
945     dBinsPair[4]       = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
946     axisTitlePair[4]   = "#Delta#eta"; 
947
948     dBinsPair[5]       = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
949     axisTitlePair[5]   = "#Delta#varphi (rad)";  
950
951     Int_t dim_val=6;
952
953     if(fRequestEventPlane){
954     dBinsPair[dim_val]       = GetBinning(fBinningString, "eventPlane", iBinPair[dim_val]);
955     axisTitlePair[dim_val]   = "#varphi - #Psi_{2} (a.u.)";
956     dim_val=7;
957     }
958
959     if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
960     dBinsPair[dim_val]       = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val]);
961     axisTitlePair[dim_val]   = "TrigCharge";
962
963     dBinsPair[dim_val+1]       = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+1]);
964     axisTitlePair[dim_val+1]   = "AssoCharge";
965     }
966
967  if(fcontainPIDtrig && !fcontainPIDasso){
968     dBinsPair[dim_val]       = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
969     axisTitlePair[dim_val]   = "PIDTrig"; 
970     if(SetChargeAxis==2){
971     dBinsPair[dim_val+1]       = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
972     axisTitlePair[dim_val+1]   = "TrigCharge";
973
974     dBinsPair[dim_val+2]       = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
975     axisTitlePair[dim_val+2]   = "AssoCharge";
976     }
977  }
978
979  if(!fcontainPIDtrig && fcontainPIDasso){
980     dBinsPair[dim_val]       = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val]);
981     axisTitlePair[dim_val]   = "PIDAsso"; 
982
983  if(SetChargeAxis==2){
984     dBinsPair[dim_val+1]       = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
985     axisTitlePair[dim_val+1]   = "TrigCharge";
986
987     dBinsPair[dim_val+2]       = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
988     axisTitlePair[dim_val+2]   = "AssoCharge";
989     }
990  }
991
992 if(fcontainPIDtrig && fcontainPIDasso){
993
994     dBinsPair[dim_val]       = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
995     axisTitlePair[dim_val]   = "PIDTrig";
996
997     dBinsPair[dim_val+1]       = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val+1]);
998     axisTitlePair[dim_val+1]   = "PIDAsso";
999
1000     if(SetChargeAxis==2){
1001     dBinsPair[dim_val+2]       = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+2]);
1002     axisTitlePair[dim_val+2]   = "TrigCharge";
1003
1004     dBinsPair[dim_val+3]       = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+3]);
1005     axisTitlePair[dim_val+3]   = "AssoCharge";
1006     }
1007  }
1008         
1009         Int_t nEtaBin = -1;
1010         Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
1011         
1012         Int_t nPteffbin = -1;
1013         Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
1014
1015
1016         fminPtTrig=dBinsPair[2][0];
1017         fmaxPtTrig=dBinsPair[2][iBinPair[2]];
1018         fminPtAsso=dBinsPair[3][0];
1019         fmaxPtAsso=dBinsPair[3][iBinPair[3]];
1020         fmincentmult=dBinsPair[0][0];
1021         fmaxcentmult=dBinsPair[0][iBinPair[0]];
1022
1023         //event pool manager
1024 Int_t MaxNofEvents=1000;
1025 const Int_t NofVrtxBins=10+(1+10)*2;
1026 Double_t ZvrtxBins[NofVrtxBins+1]={ -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10, 
1027                                        90,  92,  94,  96,  98, 100, 102, 104, 106, 108, 110, 
1028                                     190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210}; 
1029
1030 if(fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
1031    {
1032  const Int_t NofCentBins=10;
1033  Double_t CentralityBins[NofCentBins+1]={0.,9.,14.,19.,26.,34.,44.,58.,80.,500.,1000.};//Is This binning is fine for pp, or we don't require them....
1034 if(fRequestEventPlane){
1035     // Event plane angle (Psi) bins
1036   /*
1037     Double_t* psibins = NULL;
1038     Int_t nPsiBins; 
1039     psibins = GetBinning(fBinningString, "eventPlane", nPsiBins);
1040   */
1041  const Int_t  nPsiBins=6;
1042  Double_t psibins[nPsiBins+1]={0.0*TMath::DegToRad(), 30.0*TMath::DegToRad(), 60.0*TMath::DegToRad(), 90.0*TMath::DegToRad(), 120.0*TMath::DegToRad(),150.0*TMath::DegToRad(),180.1*TMath::DegToRad()};
1043 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1044 // if(psibins)  delete [] psibins; 
1045                                     }
1046
1047  else{
1048  const Int_t  nPsiBinsd=1;
1049  Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1050 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1051
1052 // fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
1053  }  
1054 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1055
1056    }
1057  else
1058    {
1059  const Int_t  NofCentBins=15;
1060 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
1061  if(fRequestEventPlane){
1062     // Event plane angle (Psi) bins
1063    /*
1064     Double_t* psibins = NULL;
1065     Int_t nPsiBins; 
1066     psibins = GetBinning(fBinningString, "eventPlane", nPsiBins);
1067    */
1068  const Int_t  nPsiBins=6;
1069  Double_t psibins[nPsiBins+1]={0.0*TMath::DegToRad(), 30.0*TMath::DegToRad(), 60.0*TMath::DegToRad(), 90.0*TMath::DegToRad(), 120.0*TMath::DegToRad(),150.0*TMath::DegToRad(),180.1*TMath::DegToRad()};
1070 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1071 // if(psibins)  delete [] psibins; 
1072                                     }
1073
1074  else{
1075 const Int_t  nPsiBinsd=1;
1076  Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1077 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1078
1079 //fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
1080  }  
1081 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1082    }
1083
1084  
1085    if(!fPoolMgr){
1086       AliError("Event Mixing required, but Pool Manager not initialized...");
1087       return;
1088     }
1089
1090         //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
1091         //fmaxPtComboeff=fmaxPtTrig;
1092 //THnSparses for calculation of efficiency
1093
1094  if((fAnalysisType =="MCAOD") && ffillefficiency) {
1095 TString Histrename;
1096   Int_t effbin[4];
1097   effbin[0]=iBinPair[0];
1098   effbin[1]=iBinPair[1];
1099   effbin[2]=nPteffbin;
1100   effbin[3]=nEtaBin;
1101   Int_t effsteps=5;//for each species type::primMCParticles(0),primRecoTracksMatched(1),allRecoTracksMatched(2),primRecoTracksMatchedPID(3),allRecoTracksMatchedPID(4)
1102 for(Int_t jj=0;jj<6;jj++)//PID type binning
1103     {
1104      if(jj==5) effsteps=3;//for unidentified particles
1105   Histrename="fTrackHistEfficiency";Histrename+=jj;
1106   fTrackHistEfficiency[jj] = new AliTHn(Histrename.Data(), "Tracking efficiency", effsteps, 4, effbin);
1107   fTrackHistEfficiency[jj]->SetBinLimits(0, dBinsPair[0]);
1108   fTrackHistEfficiency[jj]->SetVarTitle(0, "Centrality");
1109   fTrackHistEfficiency[jj]->SetBinLimits(1, dBinsPair[1]);
1110   fTrackHistEfficiency[jj]->SetVarTitle(1, "zvtx");
1111   fTrackHistEfficiency[jj]->SetBinLimits(2, Pteff);
1112   fTrackHistEfficiency[jj]->SetVarTitle(2, "p_{T} (GeV/c)");
1113   fTrackHistEfficiency[jj]->SetBinLimits(3, EtaBin);
1114   fTrackHistEfficiency[jj]->SetVarTitle(3, "#eta");
1115   fOutput->Add(fTrackHistEfficiency[jj]);
1116     }
1117  }
1118
1119 //AliThns for Correlation plots(data &  MC)
1120  
1121      if(ffilltrigassoUNID)
1122        {
1123     fTHnCorrUNID = new AliTHn("fTHnCorrUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1124 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1125     fTHnCorrUNID->SetBinLimits(j, dBinsPair[j]);
1126     fTHnCorrUNID->SetVarTitle(j, axisTitlePair[j]);
1127   }
1128   fOutput->Add(fTHnCorrUNID);
1129
1130  fTHnCorrUNIDmix = new AliTHn("fTHnCorrUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1131 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1132     fTHnCorrUNIDmix->SetBinLimits(j, dBinsPair[j]);
1133     fTHnCorrUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1134   }
1135   fOutput->Add(fTHnCorrUNIDmix);
1136        }
1137
1138      if(ffilltrigIDassoID)
1139        {
1140 fTHnCorrID = new AliTHn("fTHnCorrID", title, anaSteps, kTrackVariablesPair, iBinPair);
1141 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1142     fTHnCorrID->SetBinLimits(j, dBinsPair[j]);
1143     fTHnCorrID->SetVarTitle(j, axisTitlePair[j]);
1144   }
1145   fOutput->Add(fTHnCorrID);
1146
1147 fTHnCorrIDmix = new AliTHn("fTHnCorrIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1148 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1149     fTHnCorrIDmix->SetBinLimits(j, dBinsPair[j]);
1150     fTHnCorrIDmix->SetVarTitle(j, axisTitlePair[j]);
1151   }
1152   fOutput->Add(fTHnCorrIDmix);
1153        }
1154
1155      if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
1156        {
1157 fTHnCorrIDUNID = new AliTHn("fTHnCorrIDUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1158 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1159     fTHnCorrIDUNID->SetBinLimits(j, dBinsPair[j]);
1160     fTHnCorrIDUNID->SetVarTitle(j, axisTitlePair[j]);
1161   }
1162   fOutput->Add(fTHnCorrIDUNID);
1163
1164
1165 fTHnCorrIDUNIDmix = new AliTHn("fTHnCorrIDUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1166 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1167     fTHnCorrIDUNIDmix->SetBinLimits(j, dBinsPair[j]);
1168     fTHnCorrIDUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1169   }
1170   fOutput->Add(fTHnCorrIDUNIDmix);
1171        }
1172
1173
1174
1175   //ThnSparse for Correlation plots(truth MC)
1176      if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
1177
1178 fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
1179 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1180     fCorrelatonTruthPrimary->SetBinLimits(j, dBinsPair[j]);
1181     fCorrelatonTruthPrimary->SetVarTitle(j, axisTitlePair[j]);
1182   }
1183   fOutput->Add(fCorrelatonTruthPrimary);
1184
1185
1186 fCorrelatonTruthPrimarymix = new AliTHn("fCorrelatonTruthPrimarymix", title, anaSteps, kTrackVariablesPair, iBinPair);
1187 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1188     fCorrelatonTruthPrimarymix->SetBinLimits(j, dBinsPair[j]);
1189     fCorrelatonTruthPrimarymix->SetVarTitle(j, axisTitlePair[j]);
1190   }
1191   fOutput->Add(fCorrelatonTruthPrimarymix);     
1192  }
1193
1194     //binning for trigger no. counting
1195
1196      Int_t ChargeAxis=0;
1197      if(SetChargeAxis==2) ChargeAxis=1;
1198
1199         Int_t* fBinst;
1200         Int_t dims=3+ChargeAxis+eventplaneaxis;
1201         if(fcontainPIDtrig) dims=4+ChargeAxis+eventplaneaxis;
1202         fBinst= new Int_t[dims];
1203    Double_t* dBinsTrig[dims];    // bins for track variables  
1204    TString* axisTitleTrig;  // axis titles for track variables
1205    axisTitleTrig=new TString[dims];
1206
1207         for(Int_t i=0; i<3;i++)
1208           {
1209             fBinst[i]=iBinPair[i];
1210             dBinsTrig[i]=dBinsPair[i];
1211             axisTitleTrig[i]=axisTitlePair[i];
1212           }
1213         Int_t dim_val_trig=3;
1214     if(fRequestEventPlane){
1215       fBinst[dim_val_trig]=iBinPair[6];//if fRequestEventPlane=TRUE, dim_val already becomes 7.
1216       dBinsTrig[dim_val_trig]=dBinsPair[6];
1217       axisTitleTrig[dim_val_trig]=axisTitlePair[6];
1218       dim_val_trig=4;
1219     }
1220
1221 if(!fcontainPIDtrig && !fcontainPIDasso && ChargeAxis==1){
1222 fBinst[dim_val_trig]=iBinPair[dim_val];
1223 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1224 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1225     }
1226
1227 if(fcontainPIDtrig && !fcontainPIDasso){
1228 fBinst[dim_val_trig]=iBinPair[dim_val];
1229 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1230 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val]; 
1231     if(ChargeAxis==1){
1232 fBinst[dim_val_trig+1]=iBinPair[dim_val+1];
1233 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+1];
1234 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+1];
1235     }
1236  }
1237
1238  if(!fcontainPIDtrig && fcontainPIDasso){
1239  if(ChargeAxis==1){
1240     fBinst[dim_val_trig]=iBinPair[dim_val+1];
1241 dBinsTrig[dim_val_trig]=dBinsPair[dim_val+1];
1242 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val+1];
1243     }
1244  }
1245
1246 if(fcontainPIDtrig && fcontainPIDasso){
1247   fBinst[dim_val_trig]=iBinPair[dim_val];
1248 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1249 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val]; 
1250     if(ChargeAxis==1){
1251 fBinst[dim_val_trig+1]=iBinPair[dim_val+2];
1252 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+2];
1253 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+2];
1254     }
1255     }
1256  
1257   //ThSparse for trigger counting(data & reco MC)
1258   if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
1259           {
1260             fTHnTrigcount = new  AliTHn("fTHnTrigcount", "fTHnTrigcount", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1261    for(Int_t i=0; i<dims;i++){
1262     fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
1263     fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
1264   } 
1265   fOutput->Add(fTHnTrigcount);
1266           }
1267   
1268   if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
1269   //AliTHns for trigger counting(truth MC)
1270   fTHnTrigcountMCTruthPrim = new  AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1271  for(Int_t i=0; i<dims;i++){
1272     fTHnTrigcountMCTruthPrim->SetBinLimits(i, dBinsTrig[i]);
1273     fTHnTrigcountMCTruthPrim->SetVarTitle(i, axisTitleTrig[i]);
1274   } 
1275   fOutput->Add(fTHnTrigcountMCTruthPrim);
1276  }
1277
1278 if(fAnalysisType=="MCAOD"){
1279   if(ffillhistQATruth)
1280     {
1281   MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
1282   fOutputList->Add(MCtruthpt);
1283
1284   MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
1285   fOutputList->Add(MCtrutheta);
1286
1287   MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
1288   fOutputList->Add(MCtruthphi);
1289
1290   MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
1291   fOutputList->Add(MCtruthpionpt);
1292
1293   MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
1294   fOutputList->Add(MCtruthpioneta);
1295
1296   MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
1297   fOutputList->Add(MCtruthpionphi);
1298
1299   MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
1300   fOutputList->Add(MCtruthkaonpt);
1301
1302   MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
1303   fOutputList->Add(MCtruthkaoneta);
1304
1305   MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
1306   fOutputList->Add(MCtruthkaonphi);
1307
1308   MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
1309   fOutputList->Add(MCtruthprotonpt);
1310
1311   MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
1312   fOutputList->Add(MCtruthprotoneta);
1313
1314   MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
1315   fOutputList->Add(MCtruthprotonphi);
1316     }
1317  fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
1318   fOutputList->Add(fPioncont);
1319
1320  fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
1321   fOutputList->Add(fKaoncont);
1322
1323  fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
1324   fOutputList->Add(fProtoncont);
1325
1326 fUNIDcont=new TH2F("fUNIDcont","fUNIDcont",10,-0.5,9.5,100,0.,10.);
1327   fOutputList->Add(fUNIDcont);
1328   }
1329
1330 fEventno=new TH2F("fEventno","fEventno",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1331  fEventno->GetXaxis()->SetTitle("Centrality");
1332  fEventno->GetYaxis()->SetTitle("Z_Vtx");
1333 fOutput->Add(fEventno);
1334 fEventnobaryon=new TH2F("fEventnobaryon","fEventnobaryon",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1335  fEventnobaryon->GetXaxis()->SetTitle("Centrality");
1336  fEventnobaryon->GetYaxis()->SetTitle("Z_Vtx");
1337 fOutput->Add(fEventnobaryon);
1338 fEventnomeson=new TH2F("fEventnomeson","fEventnomeson",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1339  fEventnomeson->GetXaxis()->SetTitle("Centrality");
1340  fEventnomeson->GetYaxis()->SetTitle("Z_Vtx");
1341 fOutput->Add(fEventnomeson);
1342
1343 fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
1344 fOutput->Add(fhistJetTrigestimate);
1345
1346    fTwoTrackDistancePtdip = new TH3F("fTwoTrackDistancePtdip", ";#Delta#eta;#Delta#varphi;#Delta p_{T}", 36, -1.8, 1.8, 72,-TMath::Pi()/2, 3*TMath::Pi()/2, 40, 0, 10);
1347   fOutput->Add(fTwoTrackDistancePtdip);
1348
1349 fTwoTrackDistancePtdipmix = new TH3F("fTwoTrackDistancePtdipmix", ";#Delta#eta;#Delta#varphi;#Delta p_{T}", 36, -1.8, 1.8, 72,-TMath::Pi()/2, 3*TMath::Pi()/2, 40, 0, 10);
1350   fOutput->Add(fTwoTrackDistancePtdipmix);
1351
1352   TString Histttrname;
1353 for(Int_t jj=0;jj<2;jj++)// PID type binning
1354     {
1355   Histttrname="fTwoTrackDistancePt";Histttrname+=jj;
1356   fTwoTrackDistancePt[jj] = new TH3F(Histttrname.Data(), ";#Delta#eta;#Delta#varphi^{*}_{min};#Delta p_{T}", 100, -0.15, 0.15, 100, -0.05, 0.05, 20, 0, 10);
1357   fOutput->Add(fTwoTrackDistancePt[jj]);
1358
1359  Histttrname="fTwoTrackDistancePtmix";Histttrname+=jj;
1360   fTwoTrackDistancePtmix[jj] = new TH3F(Histttrname.Data(), ";#Delta#eta;#Delta#varphi^{*}_{min};#Delta p_{T}", 100, -0.15, 0.15, 100, -0.05, 0.05, 20, 0, 10);
1361   fOutput->Add(fTwoTrackDistancePtmix[jj]);
1362     }
1363 //Mixing
1364 //DefineEventPool();
1365
1366   if(fapplyTrigefficiency || fapplyAssoefficiency)
1367    {
1368      const Int_t nDimt = 4;//       cent zvtx  pt   eta
1369      Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
1370      Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
1371      Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
1372
1373   TString Histrexname;
1374 for(Int_t jj=0;jj<6;jj++)// PID type binning
1375     {
1376   Histrexname="effcorection";Histrexname+=jj;
1377   effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
1378   effcorection[jj]->Sumw2(); 
1379   effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
1380   effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
1381   effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
1382   effcorection[jj]->GetAxis(1)->SetTitle("zvtx"); 
1383   effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);  
1384   effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1385   effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
1386   effcorection[jj]->GetAxis(3)->SetTitle("#eta");
1387   fOutput->Add(effcorection[jj]);
1388     }
1389 // TFile *fsifile = new TFile(fefffilename,"READ");
1390
1391  if (TString(fefffilename).BeginsWith("alien:"))
1392     TGrid::Connect("alien:");
1393  TFile *fileT=TFile::Open(fefffilename);
1394  TString Nameg;
1395 for(Int_t jj=0;jj<6;jj++)//type binning
1396     {
1397 Nameg="effmap";Nameg+=jj;
1398 //effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
1399 effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
1400
1401 //effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
1402     }
1403 //fsifile->Close();
1404 fileT->Close();
1405
1406    }
1407
1408   //*************************************************************EP plots***********************************************//
1409   if(fRequestEventPlane){
1410   // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
1411   // v2
1412   fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
1413   fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
1414   fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
1415
1416   fList->Add(fHResTPCv0A2);
1417   fList->Add(fHResTPCv0C2);
1418   fList->Add(fHResv0Cv0A2);
1419
1420   // v3
1421   fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
1422   fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
1423   fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
1424
1425   fList->Add(fHResTPCv0A3);
1426   fList->Add(fHResTPCv0C3);
1427   fList->Add(fHResv0Cv0A3);
1428
1429   // MC as in the dataEP resolution (but using MC tracks)
1430   if(fAnalysisType == "MCAOD"  && fV2){
1431     fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
1432     fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
1433     fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
1434     fList->Add(fHResMA2); 
1435     fList->Add(fHResMC2); 
1436     fList->Add(fHResAC2); 
1437   }
1438   if(fAnalysisType == "MCAOD" && fV3){
1439     fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
1440     fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
1441     fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
1442     fList->Add(fHResMA3); 
1443     fList->Add(fHResMC3); 
1444     fList->Add(fHResAC3); 
1445   }
1446
1447
1448   // V0A and V0C event plane distributions
1449   //v2 
1450   fPhiRPTPC = new TH2F("fPhiRPTPCv2","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1451   fPhiRPTPCv3 = new TH2F("fPhiRPTPCv3","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1452   fList->Add(fPhiRPTPC);
1453   fList->Add(fPhiRPTPCv3);
1454
1455   fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1456   fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1457   fList->Add(fPhiRPv0A);
1458   fList->Add(fPhiRPv0C);
1459
1460   //v3
1461   fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1462   fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1463   fList->Add(fPhiRPv0Av3);
1464   fList->Add(fPhiRPv0Cv3);
1465
1466   fHistEventPlaneTruth = new TH2F("fHistEventPlaneTruth","#phi distribution of EP MCTRUTHheader;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1467   fList->Add(fHistEventPlaneTruth);
1468
1469   }
1470     
1471   //*****************************************************PIDQA histos*****************************************************//
1472
1473  
1474   //nsigma plot
1475   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1476     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1477       Double_t miny=-30;
1478       Double_t maxy=30;
1479       if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1480       TH2F *fHistoNSigma=new TH2F(Form("NSigma_%d_%d",ipart,ipid),Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1481       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1482       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1483       fOutputList->Add(fHistoNSigma);
1484     }
1485   }
1486   
1487   //nsigmaRec plot
1488   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1489     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1490       Double_t miny=-10;
1491       Double_t maxy=10;
1492       if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1493       TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
1494                                   Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1495       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1496       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1497       fOutputList->Add(fHistoNSigma);
1498     }
1499   }
1500
1501   //BayesRec plot
1502   if(fPIDType==Bayes){//use bayesianPID
1503     fPIDCombined = new AliPIDCombined();
1504     fPIDCombined->SetDefaultTPCPriors();//****************************************Need to know about it
1505
1506   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1507     Double_t miny=0.;
1508     Double_t maxy=1;
1509     TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
1510                                Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1511     fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1512     fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
1513     fOutputList->Add(fHistoBayes);
1514
1515
1516    TH2F *fHistoBayesTPC=new TH2F(Form("probBayes_TPC_%d",ipart),
1517                                Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1518     fHistoBayesTPC->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1519     fHistoBayesTPC->GetYaxis()->SetTitle(Form("Bayes prob TPC %s",kParticleSpeciesName[ipart]));
1520     fOutputList->Add(fHistoBayesTPC);
1521
1522   TH2F *fHistoBayesTOF=new TH2F(Form("probBayes_TOF_%d",ipart),
1523                                Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1524     fHistoBayesTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1525     fHistoBayesTOF->GetYaxis()->SetTitle(Form("Bayes prob TOF %s",kParticleSpeciesName[ipart]));
1526     fOutputList->Add(fHistoBayesTOF);
1527
1528  TH2F *fHistoBayesTPCTOF=new TH2F(Form("probBayes_TPCTOF_%d",ipart),
1529                                Form("probability for Tracks as  %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1530     fHistoBayesTPCTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1531     fHistoBayesTPCTOF->GetYaxis()->SetTitle(Form("Bayes prob TPCTOF %s",kParticleSpeciesName[ipart]));
1532     fOutputList->Add(fHistoBayesTPCTOF);
1533   }
1534   }
1535
1536   //nsigma separation power plot 
1537     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1538  Double_t miny=0;
1539  Double_t maxy=10;
1540    TH2F *Pi_Ka_sep=new TH2F(Form("Pi_Ka_sep_%d",ipid),
1541                                Form("Pi_Ka separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1542     Pi_Ka_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1543     Pi_Ka_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1544     fOutputList->Add(Pi_Ka_sep);
1545
1546    TH2F *Pi_Pr_sep=new TH2F(Form("Pi_Pr_sep_%d",ipid),
1547                                Form("Pi_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1548     Pi_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1549     Pi_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1550     fOutputList->Add(Pi_Pr_sep);
1551
1552     TH2F *Ka_Pr_sep=new TH2F(Form("Ka_Pr_sep_%d",ipid),
1553                                Form("Ka_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1554     Ka_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1555     Ka_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1556     fOutputList->Add(Ka_Pr_sep);
1557     }
1558
1559   //nsigmaDC plot
1560   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1561     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1562       Double_t miny=-10;
1563       Double_t maxy=10;
1564       if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1565       TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
1566                                   Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1567       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1568       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1569       fOutputList->Add(fHistoNSigma);
1570     }
1571   }
1572   
1573   //nsigmaMC plot
1574  if (fAnalysisType == "MCAOD"){
1575   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1576     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1577       Double_t miny=-30;
1578       Double_t maxy=30;
1579       if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1580       TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
1581                                   Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1582       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1583       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1584       fOutputList->Add(fHistoNSigma);
1585     }
1586   }
1587   }
1588   //PID signal plot
1589   for(Int_t idet=0;idet<fNDetectors;idet++){
1590     for(Int_t ipart=0;ipart<NSpecies;ipart++){
1591       Double_t maxy=500;
1592       if(idet==fTOF)maxy=1.1;
1593       TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
1594       fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1595       fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1596       fOutputList->Add(fHistoPID);
1597     }
1598   }
1599   //PID signal plot, before PID cut
1600   for(Int_t idet=0;idet<fNDetectors;idet++){
1601     Double_t maxy=500;
1602     if(idet==fTOF)maxy=1.1;
1603     TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
1604     fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1605     fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1606     fOutputList->Add(fHistoPID);
1607   }
1608
1609   PostData(1, fOutput);              // Post data for ALL output slots >0 here, to get at least an empty histogram
1610   PostData(2, fOutputList);
1611   if(fRequestEventPlane) PostData(3, fList);
1612   AliInfo("Finished setting up the Output");
1613
1614   TH1::AddDirectory(oldStatus);
1615
1616
1617
1618 }
1619 //-------------------------------------------------------------------------------
1620 void AliTwoParticlePIDCorr::UserExec( Option_t * ){
1621
1622  
1623   if(fAnalysisType == "AOD") {
1624
1625     doAODevent();
1626     
1627   }//AOD--analysis-----
1628
1629   else if(fAnalysisType == "MCAOD") {
1630   
1631     doMCAODevent();
1632     
1633   }
1634   
1635   else return;
1636   
1637 }
1638 //-------------------------------------------------------------------------
1639 void AliTwoParticlePIDCorr::doMCAODevent() 
1640 {
1641   AliVEvent *event = InputEvent();
1642   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1643   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1644   if (!aod) {
1645     AliError("Cannot get the AOD event");
1646     return;
1647   }
1648  
1649 // count all events(physics triggered)   
1650   fEventCounter->Fill(1);
1651
1652     evplaneMC=999.;
1653     fgPsi2v0aMC=999.;
1654     fgPsi2v0cMC=999.;
1655     fgPsi2tpcMC=999.;
1656     fgPsi3v0aMC=999.;
1657     fgPsi3v0cMC=999.;
1658     fgPsi3tpcMC=999.;
1659     gReactionPlane = 999.;
1660
1661  // get centrality object and check quality(valid for p-Pb and Pb-Pb; coming soon for pp 7 TeV)
1662   Double_t cent_v0=-1.0;
1663   Double_t effcent=1.0;
1664   Double_t refmultReco =0.0;
1665
1666 //check the PIDResponse handler
1667   if (!fPID) return;
1668
1669 // get mag. field required for twotrack efficiency cut
1670  Float_t bSign = 0;
1671  bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
1672
1673  //check for TClonesArray(truth track MC information)
1674 fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1675   if (!fArrayMC) {
1676     AliFatal("Error: MC particles branch not found!\n");
1677     return;
1678   }
1679   
1680   //check for AliAODMCHeader(truth event MC information)
1681   AliAODMCHeader *header=NULL;
1682   header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());  
1683   if(!header) {
1684     printf("MC header branch not found!\n");
1685     return;
1686   }
1687  
1688 //Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
1689 Float_t zVtxmc =header->GetVtxZ();
1690  if(TMath::Abs(zVtxmc)>fzvtxcut) return;
1691
1692  // For productions with injected signals, figure out above which label to skip particles/tracks
1693
1694  if (fInjectedSignals)
1695   {
1696     AliGenEventHeader* eventHeader = 0;
1697     Int_t headers = 0;
1698
1699 // AOD
1700       if (!header)
1701       AliFatal("fInjectedSignals set but no MC header found");
1702       
1703       headers = header->GetNCocktailHeaders();
1704       eventHeader = header->GetCocktailHeader(0);
1705
1706  if (!eventHeader)
1707     {
1708       // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing 
1709       // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
1710       AliError("First event header not found. Skipping this event.");
1711       //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1712       return;
1713     }
1714 skipParticlesAbove = eventHeader->NProduced();
1715     AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
1716   }
1717
1718  if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
1719    {
1720  //make the event selection with reco vertex cut and centrality cut and return the value of the centrality
1721      Double_t refmultTruth = GetAcceptedEventMultiplicity(aod,kTRUE);  //incase of ref multiplicity it will return the truth MC ref mullt value; need to determine the ref mult value separately for reco Mc case; in case of centrality this is final and fine
1722      refmultReco = GetAcceptedEventMultiplicity(aod,kFALSE); 
1723      if(refmultTruth<=0 || refmultReco<=0) return;
1724      cent_v0=refmultTruth;
1725    }
1726  else {
1727  cent_v0=GetAcceptedEventMultiplicity(aod,kFALSE); //centrality value; 2nd argument has no meaning
1728  }
1729
1730  if(cent_v0<0.) return;
1731  effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1732
1733   //get the event plane in case of PbPb
1734    if(fRequestEventPlane){
1735    gReactionPlane=GetEventPlane(aod,kTRUE,cent_v0);//get the truth event plane
1736    if(gReactionPlane==999.) return;
1737  }
1738   
1739
1740    Double_t nooftrackstruth=0.0;//in case of pp this will give the multiplicity(for truth case) after the track loop(only for unidentified particles that pass  kinematic cuts)
1741
1742    //TObjArray* tracksMCtruth=0;
1743 TObjArray* tracksMCtruth=new TObjArray;//for truth MC particles with PID,here unidentified means any particle other than pion, kaon or proton(Basicaly Spundefined of AliHelperPID)******WARNING::different from data and reco MC
1744  tracksMCtruth->SetOwner(kTRUE);  //***********************************IMPORTANT!
1745
1746   eventno++;
1747
1748   //There is a small difference between zvtx and zVtxmc?????? 
1749   //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
1750   //cout<<"***********************************************zvtx="<<zvtx<<endl;
1751  
1752 //now process the truth particles(for both efficiency & correlation function)
1753 Int_t nMCTrack = fArrayMC->GetEntriesFast();
1754   
1755 for (Int_t iMC = 0; iMC < nMCTrack; iMC++) 
1756 {      //MC truth track loop starts
1757     
1758 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
1759     
1760     if(!partMC){
1761       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
1762       continue;
1763     }
1764
1765 //consider only charged particles
1766     if(partMC->Charge() == 0) continue;
1767
1768 //consider only primary particles; neglect all secondary particles including from weak decays
1769  if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
1770
1771
1772 //remove injected signals(primaries above <maxLabel>)
1773  if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
1774
1775 //remove duplicates
1776   Bool_t isduplicate=kFALSE;
1777  if (fRemoveDuplicates)
1778    { 
1779  for (Int_t j=iMC+1; j<nMCTrack; ++j) 
1780    {//2nd trutuh loop starts
1781 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
1782    if(!partMC2){
1783       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
1784       continue;
1785     }    
1786  if (partMC->GetLabel() == partMC2->GetLabel())
1787    {
1788 isduplicate=kTRUE;
1789  break;  
1790    }    
1791    }//2nd truth loop ends
1792    }
1793  if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
1794
1795 //give only kinematic cuts at the truth level  
1796  if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1797  if (partMC->Pt() < fminPt ||  partMC->Pt() > fmaxPt) continue;
1798
1799  if(!partMC) continue;//for safety
1800
1801  //To determine multiplicity in case of PP
1802  nooftrackstruth++;
1803  //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1804 //only physical primary(all/unidentified)  
1805 if(ffillhistQATruth)
1806     {
1807  MCtruthpt->Fill(partMC->Pt());
1808  MCtrutheta->Fill(partMC->Eta());
1809  MCtruthphi->Fill(partMC->Phi());
1810     }
1811  //get particle ID
1812 Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
1813 Int_t particletypeTruth=-999;
1814  if (TMath::Abs(pdgtruth)==211)
1815    {
1816  particletypeTruth=SpPion;
1817 if(ffillhistQATruth)
1818     {
1819  MCtruthpionpt->Fill(partMC->Pt());
1820  MCtruthpioneta->Fill(partMC->Eta());
1821  MCtruthpionphi->Fill(partMC->Phi());
1822     }
1823       }
1824  if (TMath::Abs(pdgtruth)==321)
1825    {
1826  particletypeTruth=SpKaon;
1827 if(ffillhistQATruth)
1828     {
1829  MCtruthkaonpt->Fill(partMC->Pt());
1830  MCtruthkaoneta->Fill(partMC->Eta());
1831  MCtruthkaonphi->Fill(partMC->Phi());
1832   }
1833     }
1834 if(TMath::Abs(pdgtruth)==2212)
1835   {
1836  particletypeTruth=SpProton;
1837 if(ffillhistQATruth)
1838     {
1839  MCtruthprotonpt->Fill(partMC->Pt());
1840  MCtruthprotoneta->Fill(partMC->Eta());
1841  MCtruthprotonphi->Fill(partMC->Phi());
1842     }
1843      }
1844  if(TMath::Abs(pdgtruth)!=211 && TMath::Abs(pdgtruth)!=321 && TMath::Abs(pdgtruth)!=2212)  particletypeTruth=unidentified;//*********************WARNING:: situation is different from reco MC and data case(here we don't have SpUndefined particles,because here unidentified=SpUndefined)
1845
1846     if(fRequestEventPlane){
1847       FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
1848     }
1849
1850  // -- Fill THnSparse for efficiency and contamination calculation
1851     if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) effcent=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicity i.e. multiplicity intregated out )**************Important
1852
1853  Double_t primmctruth[4] = {effcent, zVtxmc,partMC->Pt(), partMC->Eta()};
1854  if(ffillefficiency)
1855   {
1856     fTrackHistEfficiency[5]->Fill(primmctruth,0);//for all primary truth particles(4)
1857     if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[3]->Fill(primmctruth,0);//for  primary truth mesons(3)
1858     if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[4]->Fill(primmctruth,0);//for  primary truth kaons+protons(4)
1859     if (TMath::Abs(pdgtruth)==211)  fTrackHistEfficiency[0]->Fill(primmctruth,0);//for pions
1860     if (TMath::Abs(pdgtruth)==321)  fTrackHistEfficiency[1]->Fill(primmctruth,0);//for kaons
1861     if(TMath::Abs(pdgtruth)==2212)  fTrackHistEfficiency[2]->Fill(primmctruth,0);//for protons
1862   }
1863    
1864  Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1865 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1866   {
1867     Short_t chargeval=0;
1868     if(partMC->Charge()>0)   chargeval=1;
1869     if(partMC->Charge()<0)   chargeval=-1;
1870     if(chargeval==0) continue;
1871     const TBits *clustermap=0;
1872     const TBits *sharemap=0;
1873     LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth,clustermap,sharemap);
1874 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1875  copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1876  tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1877   }
1878   }//MC truth track loop ends
1879
1880 //*********************still in event loop
1881
1882    if (fRandomizeReactionPlane)//only for TRuth MC??
1883   {
1884     Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
1885     Double_t angle = TMath::TwoPi() * centralityDigits;
1886     AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
1887     ShiftTracks(tracksMCtruth, angle);  
1888   }
1889  
1890
1891  Float_t weghtval=1.0;
1892
1893 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1894   {
1895  //Fill Correlations for MC truth particles(same event)
1896 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1897   Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1898
1899 //start mixing
1900  AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
1901 if (pool2 && pool2->IsReady())
1902   {//start mixing only when pool->IsReady
1903 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1904   {//proceed only when no. of trigger particles >0 in current event
1905     Float_t nmix=(Float_t)pool2->GetCurrentNEvents();  
1906 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++) 
1907   { //pool event loop start
1908  TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1909   if(!bgTracks6) continue;
1910   Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1911   
1912    }// pool event loop ends mixing case
1913  }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1914 } //if pool->IsReady() condition ends mixing case
1915
1916  //still in main event loop
1917
1918  if(tracksMCtruth){
1919 if(pool2)  pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1920  }
1921   }
1922
1923  //still in main event loop
1924
1925 if(tracksMCtruth) delete tracksMCtruth;
1926
1927 //now deal with reco tracks
1928
1929
1930    Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
1931
1932 //detrmine the ref mult in case of Reco(not required if we get centrality info from AliCentrality)
1933  if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) cent_v0=refmultReco;
1934  effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1935
1936  if(fRequestEventPlane){
1937    gReactionPlane = GetEventPlane(aod,kFALSE,cent_v0);//get the reconstructed event plane
1938     if(gReactionPlane==999.) return;
1939  }
1940
1941
1942   TExMap *trackMap = new TExMap();
1943
1944 // --- track loop for mapping matrix
1945  if(fFilterBit==128)
1946    {
1947  for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++) 
1948 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation 
1949   AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
1950   if (!track) continue;
1951   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
1952   if(tracktype!=1) continue; 
1953
1954   if(!track) continue;//for safety
1955
1956    Int_t gid = track->GetID();
1957    trackMap->Add(gid,itrk);
1958  }//track looop ends
1959    }
1960
1961
1962   
1963    //TObjArray* tracksUNID=0;
1964    TObjArray* tracksUNID = new TObjArray;
1965    tracksUNID->SetOwner(kTRUE);
1966
1967    //TObjArray* tracksID=0;
1968    TObjArray* tracksID = new TObjArray;
1969    tracksID->SetOwner(kTRUE);
1970
1971
1972
1973
1974    Double_t trackscount=0.0;
1975 // loop over reconstructed tracks 
1976   for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++) 
1977 { // reconstructed track loop starts
1978   AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1979   if (!track) continue;
1980  //get the corresponding MC track at the truth level (doing reco matching)
1981   AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
1982   if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
1983
1984 //remove injected signals 
1985  if(fInjectedSignals)
1986    {
1987     AliAODMCParticle* mother = recomatched;
1988
1989       while (!mother->IsPhysicalPrimary())
1990       {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
1991         if (mother->GetMother() < 0)
1992         {
1993           mother = 0;
1994           break;
1995         }
1996           
1997    mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1998         if (!mother)
1999           break;
2000       }
2001  if (!mother)
2002     {
2003       Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
2004       continue;
2005     }
2006  if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
2007    }//remove injected signals
2008
2009  if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
2010         
2011   Bool_t isduplicate2=kFALSE;
2012 if (fRemoveDuplicates)
2013    {
2014   for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++) 
2015     {//2nd loop starts
2016  AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
2017  if (!track2) continue;
2018  AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
2019 if(!recomatched2) continue;
2020
2021 if (track->GetLabel() == track2->GetLabel())
2022    {
2023 isduplicate2=kTRUE;
2024  break;  
2025    }
2026     }//2nd loop ends
2027    }
2028  if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
2029      
2030   fHistQA[11]->Fill(track->GetTPCNcls());
2031   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2032
2033  if(tracktype==0) continue; 
2034  if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
2035 {
2036   if(!track) continue;//for safety
2037  //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
2038
2039   AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
2040
2041   if(fFilterBit==128){
2042 Int_t gid1 = track->GetID();
2043 //if(gid1>=0) PIDtrack = track;
2044  PIDtrack = aod->GetTrack(trackMap->GetValue(-1-gid1));
2045 if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
2046   }
2047
2048   trackscount++;
2049
2050 //check for eta , phi holes
2051  fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2052  fphiSpectraasso->Fill(track->Phi(),track->Pt());
2053
2054   Int_t particletypeMC=-9999;
2055
2056 //tag all particles as unidentified
2057  particletypeMC=unidentified;
2058
2059  Float_t effmatrix=1.;
2060
2061 // -- Fill THnSparse for efficiency calculation
2062  if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) effcent=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicity i.e. multiplicity intregated out )**************Important
2063  //NOTE:: this will be used for fillinfg THnSparse of efficiency & also to get the the track by track eff. factor on the fly(only in case of pp)
2064
2065  //Clone & Reduce track list(TObjArray) for unidentified particles
2066 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2067   {
2068      Short_t chargeval=0;
2069     if(track->Charge()>0)   chargeval=1;
2070     if(track->Charge()<0)   chargeval=-1;
2071     if(chargeval==0) continue;
2072  if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2073    effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);
2074  LRCParticlePID* copy = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2075    copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
2076    tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2077   }
2078
2079 //get the pdg code of the corresponding truth particle
2080  Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
2081
2082  Double_t allrecomatchedpid[4] = {effcent, zVtxmc,recomatched->Pt(), recomatched->Eta()};
2083  if(ffillefficiency) {
2084 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,2);//for allreco matched
2085  if(TMath::Abs(pdgCode)==211 ||  TMath::Abs(pdgCode)==321)   fTrackHistEfficiency[3]->Fill(allrecomatchedpid,2);//for mesons
2086  if(TMath::Abs(pdgCode)==321 ||  TMath::Abs(pdgCode)==2212)   fTrackHistEfficiency[4]->Fill(allrecomatchedpid,2);//for kaons+protons
2087  if(TMath::Abs(pdgCode)==211)  fTrackHistEfficiency[0]->Fill(allrecomatchedpid,2);//for pions  
2088  if(TMath::Abs(pdgCode)==321)  fTrackHistEfficiency[1]->Fill(allrecomatchedpid,2);//for kaons
2089  if(TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,2);//for protons
2090
2091  if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) {
2092  fTrackHistEfficiency[5]->Fill(allrecomatchedpid,1);//for primreco matched
2093  if(TMath::Abs(pdgCode)==211 ||  TMath::Abs(pdgCode)==321)   fTrackHistEfficiency[3]->Fill(allrecomatchedpid,1);//for mesons
2094  if(TMath::Abs(pdgCode)==321 ||  TMath::Abs(pdgCode)==2212)   fTrackHistEfficiency[4]->Fill(allrecomatchedpid,1);//for kaons+protons
2095  if( TMath::Abs(pdgCode)==211)  fTrackHistEfficiency[0]->Fill(allrecomatchedpid,1);//for pions  
2096  if( TMath::Abs(pdgCode)==321)  fTrackHistEfficiency[1]->Fill(allrecomatchedpid,1);//for kaons
2097  if( TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,1);//for protons
2098  }
2099  }
2100
2101  //now start the particle identification process:)
2102
2103 Float_t dEdx = PIDtrack->GetTPCsignal();
2104  fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2105
2106  if(HasTOFPID(PIDtrack))
2107 {
2108 Double_t beta = GetBeta(PIDtrack);
2109 fHistoTOFbeta->Fill(track->Pt(), beta);
2110  }
2111
2112 //do track identification(nsigma method)
2113   particletypeMC=GetParticle(PIDtrack,fFIllPIDQAHistos);//******************************problem is here
2114 switch(TMath::Abs(pdgCode)){
2115   case 2212:
2116     if(fFIllPIDQAHistos){
2117       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2118         if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2119         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
2120         h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
2121       }
2122     }
2123     break;
2124   case 321:
2125     if(fFIllPIDQAHistos){
2126       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2127         if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2128         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpKaon,ipid));
2129         h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
2130       }
2131     }
2132     break;
2133   case 211:
2134     if(fFIllPIDQAHistos){
2135       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2136         if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2137         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpPion,ipid));
2138         h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
2139       }
2140     }
2141     break;
2142   }
2143
2144
2145 //2-d TPCTOF map(for each Pt interval)
2146   if(HasTOFPID(PIDtrack)){
2147  fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2148  fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2149  fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]); 
2150   }
2151
2152  //Pt, Eta , Phi distribution of the reconstructed identified particles
2153 if(ffillhistQAReco)
2154     {
2155 if (particletypeMC==SpPion)
2156   {
2157     fPionPt->Fill(track->Pt());
2158     fPionEta->Fill(track->Eta());
2159     fPionPhi->Fill(track->Phi());
2160   }
2161 if (particletypeMC==SpKaon)
2162   {
2163     fKaonPt->Fill(track->Pt());
2164     fKaonEta->Fill(track->Eta());
2165     fKaonPhi->Fill(track->Phi());
2166   }
2167 if (particletypeMC==SpProton)
2168   {
2169     fProtonPt->Fill(track->Pt());
2170     fProtonEta->Fill(track->Eta());
2171     fProtonPhi->Fill(track->Phi());
2172   }
2173     }
2174
2175  //for misidentification fraction calculation(do it with tuneonPID)
2176  if(particletypeMC==SpPion )
2177    {
2178      if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
2179      if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
2180      if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
2181      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
2182    }
2183 if(particletypeMC==SpKaon )
2184    {
2185      if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
2186      if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
2187      if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
2188      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
2189    }
2190  if(particletypeMC==SpProton )
2191    {
2192      if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
2193      if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
2194      if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
2195      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
2196    }
2197  if(particletypeMC==SpUndefined )
2198    {
2199      if(TMath::Abs(pdgCode)==211) fUNIDcont->Fill(1.,track->Pt());
2200      if(TMath::Abs(pdgCode)==321) fUNIDcont->Fill(3.,track->Pt());
2201      if(TMath::Abs(pdgCode)==2212) fUNIDcont->Fill(5.,track->Pt());
2202      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fUNIDcont->Fill(7.,track->Pt());
2203    }
2204
2205  if(particletypeMC==SpUndefined) continue;
2206
2207     if(fRequestEventPlane){
2208       FillPIDEventPlane(cent_v0,particletypeMC,track->Phi(),gReactionPlane);
2209     }
2210
2211  //fill tracking efficiency
2212  if(ffillefficiency)
2213    {
2214  if(particletypeMC==SpPion || particletypeMC==SpKaon)
2215    {
2216      if(TMath::Abs(pdgCode)==211 ||  TMath::Abs(pdgCode)==321) {
2217        fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
2218  if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
2219      }
2220    }
2221  if(particletypeMC==SpKaon || particletypeMC==SpProton)
2222    {
2223      if(TMath::Abs(pdgCode)==321 ||  TMath::Abs(pdgCode)==2212) {
2224        fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
2225  if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
2226      }
2227    }
2228  if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211)  {
2229    fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
2230  if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
2231  } 
2232  if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
2233    fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
2234 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
2235  }
2236  if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
2237    fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
2238 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
2239  }
2240    }
2241
2242 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2243   {
2244     Short_t chargeval=0;
2245     if(track->Charge()>0)   chargeval=1;
2246     if(track->Charge()<0)   chargeval=-1;
2247     if(chargeval==0) continue;
2248 if (fapplyTrigefficiency || fapplyAssoefficiency)
2249     effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles 
2250  LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2251     copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2252     tracksID->Add(copy1);
2253   }
2254   }// if(tracktype==1) condition structure ands
2255
2256 }//reco track loop ends
2257
2258   //*************************************************************still in event loop
2259  
2260
2261 if(trackscount>0.0)
2262   { 
2263 //fill the centrality/multiplicity distribution of the selected events
2264  fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2265
2266  if (fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE) fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
2267
2268 //count selected events having centrality betn 0-100%
2269  fEventCounter->Fill(13);
2270
2271 //***************************************event no. counting
2272 Bool_t isbaryontrig=kFALSE;
2273 Bool_t ismesontrig=kFALSE;
2274 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2275
2276 if(tracksID && tracksID->GetEntriesFast()>0)
2277   {
2278 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2279     {  //trigger loop starts
2280       LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2281       if(!trig) continue;
2282       if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2283       Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2284       if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2285       if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2286     }//trig loop ends
2287  if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx); 
2288  if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2289   }
2290
2291  //same event delte-eta, delta-phi plot
2292 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2293   {//same event calculation starts
2294     if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2295     if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID)  Fillcorrelation(gReactionPlane,tracksUNID,tracksID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
2296   }
2297
2298 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2299   {//same event calculation starts
2300     if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID)  Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2301     if(ffilltrigIDassoID)   Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2302   }
2303
2304 //still in  main event loop
2305 //start mixing
2306  if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2307   AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2308 if (pool && pool->IsReady())
2309   {//start mixing only when pool->IsReady
2310     Float_t nmix1=(Float_t)pool->GetCurrentNEvents();  
2311  for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
2312   { //pool event loop start
2313  TObjArray* bgTracks = pool->GetEvent(jMix);
2314   if(!bgTracks) continue;
2315   if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2316     Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2317  if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2318    Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE 
2319    }// pool event loop ends mixing case
2320
2321 } //if pool->IsReady() condition ends mixing case
2322  if(tracksUNID) {
2323 if(pool)
2324   pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2325  }
2326  }//mixing with unidentified particles
2327
2328  if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2329   AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2330 if (pool1 && pool1->IsReady())
2331   {//start mixing only when pool->IsReady
2332   Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();  
2333 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++) 
2334   { //pool event loop start
2335  TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2336   if(!bgTracks2) continue;
2337 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2338   Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE  
2339 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2340   Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2341
2342    }// pool event loop ends mixing case
2343 } //if pool1->IsReady() condition ends mixing case
2344
2345 if(tracksID) {
2346 if(pool1) 
2347   pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2348  }
2349  }//mixing with identified particles
2350
2351   //no. of events analyzed
2352 fEventCounter->Fill(15);
2353   }
2354
2355 if(tracksUNID)  delete tracksUNID;
2356
2357 if(tracksID) delete tracksID;
2358
2359
2360 PostData(1, fOutput);
2361
2362 }
2363 //________________________________________________________________________
2364 void AliTwoParticlePIDCorr::doAODevent() 
2365 {
2366
2367   //get AOD
2368   AliVEvent *event = InputEvent();
2369   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
2370   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
2371   if (!aod) {
2372     AliError("Cannot get the AOD event");
2373     return;
2374   }
2375
2376 // count all events   
2377   fEventCounter->Fill(1);
2378
2379 if (!fPID) return;//this should be available with each event even if we don't do PID selection
2380
2381     fgPsi2v0a=999.;
2382     fgPsi2v0c=999.;
2383     fgPsi2tpc=999.;
2384     fgPsi3v0a=999.;
2385     fgPsi3v0c=999.;
2386     fgPsi3tpc=999.;
2387     gReactionPlane = 999.;
2388     
2389   Double_t cent_v0=   -999.;
2390   Double_t effcent=1.0;
2391   Float_t bSign = 0.;
2392   Double_t trackscount=0;//counts particles passed filterbit cuts and kinematic cuts used in this analysis
2393
2394
2395  bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
2396  Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
2397
2398
2399 // check event cuts and fill event histograms and return the centrality or reference multiplicity value
2400  if((cent_v0 = GetAcceptedEventMultiplicity(aod,kFALSE)) < 0){ 
2401     return;
2402   }
2403  effcent=cent_v0;//required for efficiency correction case********Extremely Important
2404   //get the event plane in case of PbPb
2405     if(fRequestEventPlane){
2406       gReactionPlane = GetEventPlane(aod,kFALSE,cent_v0);
2407     if(gReactionPlane==999.) return;
2408     }    
2409     
2410
2411 TExMap *trackMap = new TExMap();
2412 // --- track loop for mapping matrix
2413  if(fFilterBit==128)
2414    {
2415  for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++) 
2416 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation 
2417   AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
2418   if (!track) continue;
2419   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
2420   if(tracktype!=1) continue; 
2421
2422   if(!track) continue;//for safety
2423
2424    Int_t gid = track->GetID();
2425    trackMap->Add(gid,itrk);
2426  }//track looop ends
2427    }
2428
2429    TObjArray*  tracksUNID= new TObjArray;//track info before doing PID
2430    tracksUNID->SetOwner(kTRUE);  // IMPORTANT!
2431
2432    TObjArray* tracksID= new TObjArray;//only pions, kaons,protons i.e. after doing the PID selection
2433    tracksID->SetOwner(kTRUE);  // IMPORTANT!
2434  
2435     
2436     eventno++;
2437
2438     Bool_t fTrigPtmin1=kFALSE;
2439     Bool_t fTrigPtmin2=kFALSE;
2440     Bool_t fTrigPtJet=kFALSE;
2441
2442  for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++) 
2443 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation 
2444   AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2445   if (!track) continue;
2446   fHistQA[11]->Fill(track->GetTPCNcls());
2447   Int_t particletype=-9999;//required for PID filling
2448   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2449   if(tracktype!=1) continue; 
2450
2451   if(!track) continue;//for safety
2452
2453 AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
2454
2455   if(fFilterBit==128){
2456 Int_t gid1 = track->GetID();
2457 //if(gid1>=0) PIDtrack = track;
2458  PIDtrack = aod->GetTrack(trackMap->GetValue(-1-gid1));
2459 if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
2460   }
2461
2462 //check for eta , phi holes
2463  fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2464  fphiSpectraasso->Fill(track->Phi(),track->Pt());
2465
2466  trackscount++;
2467
2468  //if no applyefficiency , set the eff factor=1.0
2469  Float_t effmatrix=1.0;
2470
2471  //tag all particles as unidentified that passed filterbit & kinematic cuts 
2472  particletype=unidentified;
2473
2474  //To count the no. of tracks having an accepted track in a certain PT(e.g. Jet Pt) range
2475  if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
2476  if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
2477  if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
2478
2479
2480 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) effcent=15.0;//integrated over multiplicity [i.e each track has multiplicity 15.0](so put any fixed value for each track so that practically means there is only one bin in multiplicityi.e multiplicity intregated out )**************Important for efficiency related issues
2481
2482
2483  //to reduce memory consumption in pool
2484   if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig)) 
2485   {
2486  //Clone & Reduce track list(TObjArray) for unidentified particles
2487     Short_t chargeval=0;
2488     if(track->Charge()>0)   chargeval=1;
2489     if(track->Charge()<0)   chargeval=-1;
2490     if(chargeval==0) continue;
2491  if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2492    effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);
2493  LRCParticlePID* copy = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2494   copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2495   tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2496   }
2497
2498 //now start the particle identificaion process:) 
2499
2500 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
2501
2502   Float_t dEdx = PIDtrack->GetTPCsignal();
2503   fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2504
2505   //fill beta vs Pt plots only for tracks having proper TOF response(much less tracks compared to the no. that pass the filterbit & kinematic cuts)
2506  if(HasTOFPID(PIDtrack))
2507 {
2508   Double_t beta = GetBeta(PIDtrack);
2509   fHistoTOFbeta->Fill(track->Pt(), beta);
2510  }
2511   
2512
2513 //track identification(using nsigma method)
2514      particletype=GetParticle(PIDtrack,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
2515
2516 //2-d TPCTOF map(for each Pt interval)
2517   if(HasTOFPID(PIDtrack)){
2518  fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2519  fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2520  fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]); 
2521   }
2522
2523 //ignore the Spundefined particles as they also contain pion, kaon, proton outside the nsigma cut(also if tracks don't have proper TOF PID in a certain Pt interval) & these tracks are actually counted when we do the the efficiency correction, so considering them as unidentified particles & doing the efficiency correction(i.e defining unidentified=pion+Kaon+proton+SpUndefined is right only without efficiency correction) for them will be two times wrong!!!!! 
2524   if (particletype==SpUndefined) continue;//this condition creating a modulated structure in delphi projection in mixed event case(only when we are dealing with identified particles i.e. tracksID)!!!!!!!!!!!
2525
2526     if(fRequestEventPlane){
2527       FillPIDEventPlane(cent_v0,particletype,track->Phi(),gReactionPlane);
2528     }
2529
2530
2531  //Pt, Eta , Phi distribution of the reconstructed identified particles
2532 if(ffillhistQAReco)
2533     {
2534 if (particletype==SpPion)
2535   {
2536     fPionPt->Fill(track->Pt());
2537     fPionEta->Fill(track->Eta());
2538     fPionPhi->Fill(track->Phi());
2539   }
2540 if (particletype==SpKaon)
2541   {
2542     fKaonPt->Fill(track->Pt());
2543     fKaonEta->Fill(track->Eta());
2544     fKaonPhi->Fill(track->Phi());
2545   }
2546 if (particletype==SpProton)
2547   {
2548     fProtonPt->Fill(track->Pt());
2549     fProtonEta->Fill(track->Eta());
2550     fProtonPhi->Fill(track->Phi());
2551   }
2552     }
2553  
2554  if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig)) 
2555   {
2556     Short_t chargeval=0;
2557     if(track->Charge()>0)   chargeval=1;
2558     if(track->Charge()<0)   chargeval=-1;
2559     if(chargeval==0) continue;
2560 if (fapplyTrigefficiency || fapplyAssoefficiency)
2561   effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles; Bool_t mesoneffrequired=kFALSE
2562  LRCParticlePID* copy1 = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2563     copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2564     tracksID->Add(copy1);
2565   }
2566 } //track loop ends but still in event loop
2567
2568 if(trackscount<1.0){
2569   if(tracksUNID) delete tracksUNID;
2570   if(tracksID) delete tracksID;
2571   return;
2572  }
2573
2574  if (fTrigPtmin1) fhistJetTrigestimate->Fill(cent_v0,0.0);
2575  if (fTrigPtmin2) fhistJetTrigestimate->Fill(cent_v0,2.0);
2576  if (fTrigPtJet) fhistJetTrigestimate->Fill(cent_v0,4.0);
2577
2578  Float_t weightval=1.0;
2579
2580
2581   
2582 //fill the centrality/multiplicity distribution of the selected events
2583  fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2584
2585 if(fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE) fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
2586
2587 //count selected events having centrality betn 0-100%
2588  fEventCounter->Fill(13);
2589
2590 //***************************************event no. counting
2591 Bool_t isbaryontrig=kFALSE;
2592 Bool_t ismesontrig=kFALSE;
2593 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2594
2595 if(tracksID && tracksID->GetEntriesFast()>0)
2596   {
2597 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2598     {  //trigger loop starts
2599       LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2600       if(!trig) continue;
2601       if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2602       Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2603       if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2604       if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2605     }//trig loop ends
2606  if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx); 
2607  if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2608   }
2609
2610
2611 //same event delta-eta-deltaphi plot 
2612
2613 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2614   {//same event calculation starts
2615     if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2616     if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID)  Fillcorrelation(gReactionPlane,tracksUNID,tracksID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
2617   }
2618
2619 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2620   {//same event calculation starts
2621     if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID)  Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2622     if(ffilltrigIDassoID)   Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2623   }
2624
2625 //still in  main event loop
2626
2627
2628 //start mixing
2629  if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2630 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2631 if (pool && pool->IsReady())
2632   {//start mixing only when pool->IsReady
2633   Float_t nmix1=(Float_t)pool->GetCurrentNEvents();  
2634  for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
2635   { //pool event loop start
2636  TObjArray* bgTracks = pool->GetEvent(jMix);
2637   if(!bgTracks) continue;
2638   if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2639     Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2640  if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2641    Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE 
2642    }// pool event loop ends mixing case
2643 } //if pool->IsReady() condition ends mixing case
2644  if(tracksUNID) {
2645 if(pool)
2646   pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2647  }
2648  }//mixing with unidentified particles
2649
2650
2651  if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2652  AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2653 if (pool1 && pool1->IsReady())
2654   {//start mixing only when pool->IsReady
2655   Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();  
2656 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++) 
2657   { //pool event loop start
2658  TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2659   if(!bgTracks2) continue;
2660 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2661   Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE  
2662 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2663   Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2664
2665    }// pool event loop ends mixing case
2666 } //if pool1->IsReady() condition ends mixing case
2667
2668 if(tracksID) {
2669 if(pool1) 
2670   pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2671  }
2672  }//mixing with identified particles
2673
2674
2675   //no. of events analyzed
2676 fEventCounter->Fill(15);
2677  
2678 //still in main event loop
2679
2680
2681 if(tracksUNID)  delete tracksUNID;
2682
2683 if(tracksID) delete tracksID;
2684
2685
2686 PostData(1, fOutput);
2687
2688 } // *************************event loop ends******************************************//_______________________________________________________________________
2689 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
2690 {
2691   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
2692   
2693   TObjArray* tracksClone = new TObjArray;
2694   tracksClone->SetOwner(kTRUE);
2695   
2696   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2697   {
2698     LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
2699     LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval(),particle->GetTPCPadMap(),particle->GetTPCSharedMap());
2700     copy100->SetUniqueID(particle->GetUniqueID());
2701     tracksClone->Add(copy100);
2702   }
2703   
2704   return tracksClone;
2705 }
2706
2707 //--------------------------------------------------------------------------------
2708 void AliTwoParticlePIDCorr::Fillcorrelation(Float_t ReactionPlane,TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t weight,Bool_t firstTime,Float_t bSign,Bool_t fPtOrder,Bool_t twoTrackEfficiencyCut,Bool_t mixcase,TString fillup)
2709 {
2710
2711   //before calling this function check that either trackstrig & tracksasso are available 
2712
2713  // Eta() is extremely time consuming, therefore cache it for the inner loop here:
2714   TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
2715   TArrayF eta(input->GetEntriesFast());
2716   for (Int_t i=0; i<input->GetEntriesFast(); i++)
2717     eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
2718
2719   //if(trackstrig)
2720     Int_t jmax=trackstrig->GetEntriesFast();
2721   if(tracksasso)
2722      jmax=tracksasso->GetEntriesFast();
2723
2724 // identify K, Lambda candidates and flag those particles
2725     // a TObject bit is used for this
2726 const UInt_t kResonanceDaughterFlag = 1 << 14;
2727     if (fRejectResonanceDaughters > 0)
2728     {
2729       Double_t resonanceMass = -1;
2730       Double_t massDaughter1 = -1;
2731       Double_t massDaughter2 = -1;
2732       const Double_t interval = 0.02;
2733  switch (fRejectResonanceDaughters)
2734       {
2735         case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
2736         case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
2737         case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
2738         default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
2739       }      
2740
2741 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2742         trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2743  for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
2744           tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2745
2746  for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2747       {
2748       LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2749 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
2750         {
2751         LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2752
2753  // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
2754 if (trig->IsEqual(asso)) continue;
2755
2756 if (trig->Charge() * asso->Charge() > 0) continue;
2757
2758  Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2759      
2760 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
2761           {
2762             mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2763
2764             if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
2765             {
2766               trig->SetBit(kResonanceDaughterFlag);
2767               asso->SetBit(kResonanceDaughterFlag);
2768               
2769 //            Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
2770             }
2771           }
2772         }
2773       }
2774     }
2775
2776       //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
2777
2778     Float_t TriggerPtMin=fminPtTrig;
2779     Float_t TriggerPtMax=fmaxPtTrig;
2780     Int_t HighestPtTriggerIndx=-99999;
2781     TH1* triggerWeighting = 0;
2782
2783 if(fSelectHighestPtTrig || fWeightPerEvent)//**************add this data member to the constructor
2784       {
2785 if (fWeightPerEvent)
2786     {
2787       TAxis* axis=0;
2788    if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID) axis = fTHnTrigcount->GetGrid(0)->GetGrid()->GetAxis(2);                                          
2789   if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH)    axis = fTHnTrigcountMCTruthPrim->GetGrid(0)->GetGrid()->GetAxis(2);
2790       triggerWeighting = new TH1F("triggerWeighting", "", axis->GetNbins(), axis->GetXbins()->GetArray());
2791     }
2792 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2793     {  //trigger loop starts(highest Pt trigger selection)
2794       LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2795       if(!trig) continue;
2796       Float_t trigpt=trig->Pt();
2797     //to avoid overflow qnd underflow
2798       if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2799       Int_t particlepidtrig=trig->getparticle();
2800       if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2801
2802       Float_t trigeta=trig->Eta();
2803
2804       // some optimization
2805  if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2806         continue;
2807
2808 if (fOnlyOneEtaSide != 0)
2809       {
2810         if (fOnlyOneEtaSide * trigeta < 0)
2811           continue;
2812       }
2813   if (fTriggerSelectCharge != 0)
2814         if (trig->Charge() * fTriggerSelectCharge < 0)
2815           continue;
2816         
2817       if (fRejectResonanceDaughters > 0)
2818         if (trig->TestBit(kResonanceDaughterFlag)) continue;
2819
2820       if(fSelectHighestPtTrig){
2821  if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
2822           {         
2823           HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
2824           TriggerPtMin=trigpt;
2825           }
2826       }
2827
2828 if (fWeightPerEvent)  triggerWeighting->Fill(trigpt);
2829
2830     }//trigger loop ends(highest Pt trigger selection)
2831
2832       }//******************(fSelectHighestPtTrig || fWeightPerEvent) condition ends
2833
2834
2835  //two particle correlation filling
2836 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2837     {  //trigger loop starts
2838       LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2839       if(!trig) continue;
2840       Float_t trigpt=trig->Pt();
2841     //to avoid overflow qnd underflow
2842       if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2843       Int_t particlepidtrig=trig->getparticle();
2844       if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2845
2846       Float_t trigeta=trig->Eta();
2847
2848       // some optimization
2849  if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2850         continue;
2851
2852 if (fOnlyOneEtaSide != 0)
2853       {
2854         if (fOnlyOneEtaSide * trigeta < 0)
2855           continue;
2856       }
2857   if (fTriggerSelectCharge != 0)
2858         if (trig->Charge() * fTriggerSelectCharge < 0)
2859           continue;
2860         
2861       if (fRejectResonanceDaughters > 0)
2862         if (trig->TestBit(kResonanceDaughterFlag)) continue;
2863
2864       if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
2865         if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
2866       }
2867
2868       Float_t trigphi=trig->Phi();
2869       Float_t trackefftrig=1.0;
2870       if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
2871
2872     // Event plane (determine psi bin)
2873     Double_t gPsiMinusPhi    =   0.;
2874     Double_t gPsiMinusPhiBin = -10.;
2875 if(fRequestEventPlane){
2876     gPsiMinusPhi   = TMath::Abs(trigphi - ReactionPlane);
2877     //in-plane(Note thet event plane angle has to be defined within 0 to 180 degree(do not use this if event ), otherwise the definition of in plane and out plane particles is wrong)
2878     if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
2879       (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
2880        ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
2881       gPsiMinusPhiBin = 0.0;
2882     /*
2883  if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
2884        ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
2885       gPsiMinusPhiBin = 0.0;
2886     */
2887     //intermediate
2888     else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
2889             ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
2890             ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
2891             ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
2892       gPsiMinusPhiBin = 1.0;
2893     //out of plane
2894     else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
2895             ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
2896       gPsiMinusPhiBin = 2.0;
2897     //everything else
2898     else 
2899       gPsiMinusPhiBin = 3.0;
2900     
2901     fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
2902  }
2903
2904       //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
2905         Double_t* trigval;
2906         Int_t dim=3;
2907         Int_t eventplaneAxis=0;
2908         if(fRequestEventPlane) eventplaneAxis=1;
2909         if(fcontainPIDtrig && SetChargeAxis==0) dim=4+eventplaneAxis;
2910         if(!fcontainPIDtrig && SetChargeAxis==2) dim=4+eventplaneAxis;
2911         if(fcontainPIDtrig && SetChargeAxis==2) dim=5+eventplaneAxis;
2912         trigval= new Double_t[dim];
2913       trigval[0] = cent;
2914       trigval[1] = vtx;
2915       trigval[2] = trigpt;
2916
2917       if(fRequestEventPlane){
2918       trigval[3] = gPsiMinusPhiBin;
2919       if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = particlepidtrig;
2920       if(!fcontainPIDtrig && SetChargeAxis==2) trigval[4] = trig->Charge();
2921       if(fcontainPIDtrig && SetChargeAxis==2) {
2922       trigval[4] = particlepidtrig;
2923       trigval[5] = trig->Charge();
2924        }
2925       }
2926
2927   if(!fRequestEventPlane){
2928       if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = particlepidtrig;
2929       if(!fcontainPIDtrig && SetChargeAxis==2) trigval[3] = trig->Charge();
2930       if(fcontainPIDtrig && SetChargeAxis==2) {
2931       trigval[3] = particlepidtrig;
2932       trigval[4] = trig->Charge();
2933        }
2934       }
2935
2936  
2937
2938         if (fWeightPerEvent)
2939         {
2940           // leads effectively to a filling of one entry per filled trigger particle pT bin
2941           Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(trigval[2]);
2942 //        Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
2943           trackefftrig *= triggerWeighting->GetBinContent(weightBin);
2944         }
2945
2946
2947       //trigger particle counting for both same and mixed event case;;;;;step=0->same event case;;;;;step=1->mixed event case;;;;;;
2948 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
2949       if(fillup=="trigassoUNID" ) {
2950 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig); 
2951 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig); 
2952       }
2953     }
2954  if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
2955    if(fillup=="trigassoUNID" )  
2956      {
2957 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig); 
2958 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig); 
2959      }
2960     }
2961 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
2962   if(fillup=="trigUNIDassoID")  
2963     {
2964 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig); 
2965 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig); 
2966     }
2967     }
2968  //ensure that trigIDassoID , trigassoUNID, trigIDassoUNID & trigUNIDassoID  case FillCorrelation called only once in the event loop for same event correlation function calculation, otherwise there will be multiple counting of pion, kaon,proton,unidentified
2969 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
2970   if(fillup=="trigIDassoID")  
2971     {
2972 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig); 
2973 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig); 
2974     }
2975     }
2976  if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
2977    if(fillup=="trigIDassoUNID" ) 
2978      {
2979 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig); 
2980 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig); 
2981      } 
2982     }
2983 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
2984   if(fillup=="trigIDassoID")  
2985     {
2986 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig); 
2987 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig); 
2988     }
2989     }
2990
2991  if(fillup=="trigIDassoIDMCTRUTH") { //In truth MC case "Unidentified" means any particle other than pion,kaon or proton and no efficiency correction(default value 1.0)************************be careful!!!! 
2992 if(mixcase==kFALSE)   fTHnTrigcountMCTruthPrim->Fill(trigval,0,1.0/trackefftrig); 
2993 if(mixcase==kTRUE && firstTime)   fTHnTrigcountMCTruthPrim->Fill(trigval,1,1.0/trackefftrig); 
2994   }
2995
2996     //asso loop starts within trigger loop
2997    for(Int_t j=0;j<jmax;j++)
2998              {
2999     LRCParticlePID *asso=0;
3000     if(!tracksasso)
3001     asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
3002     else
3003     asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
3004
3005     if(!asso) continue;
3006
3007     //to avoid overflow and underflow
3008  if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
3009
3010  //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
3011
3012   if(!tracksasso && j==i) continue;
3013
3014    // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event,i.e. both Trig and asso TObjArray belongs to the same Pt range but say Trig is Unidentified but asso is identified then the serial no. wise particles are not same and and j==i doesn't aplly)
3015    // if (tracksasso && trig->IsEqual(asso))  continue;
3016
3017   if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
3018           
3019  if (fPtOrder)
3020  if (asso->Pt() >= trig->Pt()) continue;
3021
3022   Int_t particlepidasso=asso->getparticle(); 
3023   if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
3024             
3025
3026 if (fAssociatedSelectCharge != 0)
3027 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
3028             
3029  if (fSelectCharge > 0)
3030         {
3031           // skip like sign
3032           if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
3033             continue;
3034             
3035           // skip unlike sign
3036           if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
3037             continue;
3038         }
3039
3040 if (fEtaOrdering)
3041         {
3042           if (trigeta < 0 && asso->Eta() < trigeta)
3043             continue;
3044           if (trigeta > 0 && asso->Eta() > trigeta)
3045             continue;
3046         }
3047
3048 if (fRejectResonanceDaughters > 0)
3049           if (asso->TestBit(kResonanceDaughterFlag))
3050           {
3051 //          Printf("Skipped j=%d", j);
3052             continue;
3053           }
3054
3055         // conversions
3056         if (fCutConversions && asso->Charge() * trig->Charge() < 0)
3057         {
3058           Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3059           
3060           if (mass < 0.1)
3061           {
3062             mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3063             
3064             fControlConvResoncances->Fill(0.0, mass);
3065
3066             if (mass < 0.04*0.04) 
3067               continue;
3068           }
3069         }
3070
3071         // K0s
3072         if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3073         {
3074           Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
3075           
3076           const Float_t kK0smass = 0.4976;
3077           
3078           if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
3079           {
3080             mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
3081             
3082             fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
3083
3084             if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
3085               continue;
3086           }
3087         }
3088         
3089         // Lambda
3090         if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3091         {
3092           Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
3093           Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3094           
3095           const Float_t kLambdaMass = 1.115;
3096
3097           if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
3098           {
3099             mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
3100
3101             fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
3102             
3103             if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3104               continue;
3105           }
3106           if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
3107           {
3108             mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3109
3110             fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
3111
3112             if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3113               continue;
3114           }
3115         }
3116
3117         if (twoTrackEfficiencyCut)
3118         {
3119           // the variables & cuthave been developed by the HBT group 
3120           // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
3121           Float_t phi1 = trig->Phi();
3122           Float_t pt1 = trig->Pt();
3123           Float_t charge1 = trig->Charge();
3124           Float_t phi2 = asso->Phi();
3125           Float_t pt2 = asso->Pt();
3126           Float_t charge2 = asso->Charge();
3127
3128           Float_t deta= trigeta - eta[j]; 
3129     
3130  // optimization
3131           if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
3132           {
3133
3134   // check first boundaries to see if is worth to loop and find the minimum
3135             Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMinRadius, bSign);
3136             Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMaxRadius, bSign);
3137
3138  const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
3139
3140             Float_t dphistarminabs = 1e5;
3141             Float_t dphistarmin = 1e5;
3142
3143  if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
3144             {
3145               for (Double_t rad=0.8; rad<2.51; rad+=0.01) 
3146               {
3147                 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
3148
3149                 Float_t dphistarabs = TMath::Abs(dphistar);
3150
3151         if (dphistarabs < dphistarminabs)
3152                 {
3153                   dphistarmin = dphistar;
3154                   dphistarminabs = dphistarabs;
3155                 }
3156               }
3157               if(mixcase==kFALSE)  fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3158               if(mixcase==kTRUE)  fTwoTrackDistancePtmix[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3159
3160 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
3161               {
3162 //              Printf("Removed track pair %d %d with %f %f %f %f %f %f %f %f %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);
3163                 continue;
3164               }
3165              if(mixcase==kFALSE) fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3166              if(mixcase==kTRUE) fTwoTrackDistancePtmix[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3167
3168             }
3169           }
3170         }
3171
3172         //pair sharedfraction cut(only between the trig and asso track)
3173    if(fillup!="trigIDassoIDMCTRUTH")//******************************************NOT for TRUTH MC particles
3174   {
3175     if(fSharedfraction_Pair_cut>=0){
3176         Bool_t passsharedfractionpaircut=CalculateSharedFraction(trig->GetTPCPadMap(),asso->GetTPCPadMap(),trig->GetTPCSharedMap(),asso->GetTPCSharedMap());
3177         if(!passsharedfractionpaircut) continue;
3178     }
3179   }
3180         Float_t weightperevent=weight;
3181         Float_t trackeffasso=1.0;
3182         if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
3183         //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
3184         Float_t deleta=trigeta-eta[j];
3185         Float_t delphi=PhiRange(trigphi-asso->Phi()); 
3186
3187         Float_t delpt=trigpt-asso->Pt();
3188         //fill it with/without two track efficiency cut    
3189         if(mixcase==kFALSE)  fTwoTrackDistancePtdip->Fill(deleta, delphi, TMath::Abs(delpt));//for same event
3190         if(mixcase==kTRUE)  fTwoTrackDistancePtdipmix->Fill(deleta, delphi, TMath::Abs(delpt));//for mixed event
3191
3192  //here get the two particle efficiency correction factor
3193         Float_t effweight=trackefftrig*trackeffasso*weightperevent;
3194         // if(mixcase==kFALSE)  cout<<"*******************effweight="<<effweight<<endl;
3195         Double_t* vars;
3196         Int_t dimused=kTrackVariablesPair+eventplaneAxis;
3197         vars= new Double_t[dimused];
3198         vars[0]=cent;
3199         vars[1]=vtx;
3200         vars[2]=trigpt;
3201         vars[3]=asso->Pt();
3202         vars[4]=deleta;
3203         vars[5]=delphi;
3204
3205         Int_t dimension=6;
3206         if(fRequestEventPlane) 
3207         {
3208        vars[6]=gPsiMinusPhiBin;
3209        dimension=7;
3210         }
3211
3212 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
3213         vars[dimension]=trig->Charge();
3214         vars[dimension+1]=asso->Charge();
3215  }
3216 if(fcontainPIDtrig && !fcontainPIDasso){
3217         vars[dimension]=particlepidtrig;
3218 if(SetChargeAxis==2){
3219         vars[dimension+1]=trig->Charge();
3220         vars[dimension+2]=asso->Charge();
3221  }
3222         }
3223 if(!fcontainPIDtrig && fcontainPIDasso){
3224         vars[dimension]=particlepidasso;
3225 if(SetChargeAxis==2){
3226         vars[dimension+1]=trig->Charge();
3227         vars[dimension+2]=asso->Charge();
3228    }
3229  }
3230  if(fcontainPIDtrig && fcontainPIDasso){
3231         vars[dimension]=particlepidtrig;
3232         vars[dimension+1]=particlepidasso;
3233 if(SetChargeAxis==2){
3234         vars[dimension+2]=trig->Charge();
3235         vars[dimension+3]=asso->Charge();
3236    }
3237  }
3238
3239         if (fWeightPerEvent)
3240         {
3241           Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[2]);
3242 //        Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3243          effweight *= triggerWeighting->GetBinContent(weightBin);
3244         }
3245     
3246
3247         //Fill Correlation ThnSparses
3248     if(fillup=="trigassoUNID")
3249       {
3250     if(mixcase==kFALSE)  fTHnCorrUNID->Fill(vars,0,1.0/effweight);
3251     if(mixcase==kTRUE)   fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
3252       }
3253     if(fillup=="trigIDassoID")
3254       {
3255         if(mixcase==kFALSE)  fTHnCorrID->Fill(vars,0,1.0/effweight);
3256         if(mixcase==kTRUE)  fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
3257       }
3258     if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
3259       {//in this case effweight should be 1 always 
3260         if(mixcase==kFALSE)  fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight); 
3261         if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
3262     }   
3263     if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
3264       {
3265         if(mixcase==kFALSE)  fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
3266         if(mixcase==kTRUE)   fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
3267        }
3268         
3269 delete[] vars;
3270    }//asso loop ends 
3271 delete[] trigval;
3272  }//trigger loop ends 
3273
3274  if (triggerWeighting)
3275     {
3276       delete triggerWeighting;
3277       triggerWeighting = 0;
3278     }
3279 }
3280
3281 //------------------------------------------------------------------------------------------------
3282 Bool_t AliTwoParticlePIDCorr:: CalculateSharedFraction(const TBits *triggerPadMap,const TBits *assocPadMap,const TBits *triggerShareMap,const TBits *assocShareMap)
3283 {//source code-AliFemtoShareQualityPairCut.cxx
3284 Double_t nofhits=0;
3285 Double_t nofsharedhits=0;
3286
3287 for(UInt_t imap=0;imap< (triggerPadMap->GetNbits() );imap++)
3288 {
3289 //if they are in same pad
3290 //cout<<triggerPadMap->TestBitNumber(imap)<<"    "<< assocPadMap->TestBitNumber(imap)<<endl;
3291 if (triggerPadMap->TestBitNumber(imap) &&
3292       assocPadMap->TestBitNumber(imap))
3293 {
3294 //if they share
3295 //cout<<triggerShareMap->TestBitNumber(imap)<<"   "<<assocShareMap->TestBitNumber(imap)<<endl;
3296 if (triggerShareMap->TestBitNumber(imap) &&
3297       assocShareMap->TestBitNumber(imap))
3298 {
3299   //cout<<triggerShareMap->TestBitNumber(imap)<<"   "<<assocShareMap->TestBitNumber(imap)<<endl;
3300 nofhits+=2;
3301 nofsharedhits+=2;
3302 }
3303
3304
3305
3306 //not shared
3307     else {
3308      
3309       nofhits+=2;
3310     }
3311
3312
3313 }
3314 //different pad
3315
3316 //cout<< (triggerPadMap->TestBitNumber(imap) || assocPadMap->TestBitNumber(imap))<<endl;
3317 else if (triggerPadMap->TestBitNumber(imap) ||
3318       assocPadMap->TestBitNumber(imap)) {
3319     // One track has a hit, the other does not
3320    
3321     nofhits++;
3322     //cout<<"No hits :"<<nofhits<<endl;
3323    
3324       }
3325
3326
3327
3328 }
3329
3330 Double_t SharedFraction=0.0;
3331 if(nofhits>0) SharedFraction=(nofsharedhits/nofhits);
3332
3333 //cout<<"Fraction shared hits :"<<SharedFraction<<endl;
3334
3335 if(SharedFraction>fSharedfraction_Pair_cut) return kFALSE;
3336
3337 return kTRUE;
3338
3339 }
3340
3341 //________________________________________________________________________________________________
3342 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
3343 {
3344   //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
3345  Int_t effVars[4];
3346  Float_t effvalue=1.; 
3347
3348   if(parpid==unidentified)
3349             {
3350             effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
3351             effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx); 
3352             effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt()); 
3353             effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta()); 
3354             effvalue=effcorection[5]->GetBinContent(effVars);
3355             }
3356 if(parpid==SpPion || parpid==SpKaon)
3357             {
3358               if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3359                 {
3360             effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
3361             effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx); 
3362             effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt()); 
3363             effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
3364             effvalue=effcorection[3]->GetBinContent(effVars);
3365                 }
3366               else{
3367  if(parpid==SpPion)
3368             {
3369             effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
3370             effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx); 
3371             effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt()); 
3372             effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta()); 
3373             effvalue=effcorection[0]->GetBinContent(effVars);
3374             }
3375             
3376  if(parpid==SpKaon)
3377             {
3378             effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
3379             effVars[1] =  effcorection[1]->GetAxis(1)->FindBin(evzvtx); 
3380             effVars[2] =  effcorection[1]->GetAxis(2)->FindBin(track->Pt()); 
3381             effVars[3] =  effcorection[1]->GetAxis(3)->FindBin(track->Eta()); 
3382             effvalue=effcorection[1]->GetBinContent(effVars);
3383             }
3384               }
3385             }   
3386              
3387  if(parpid==SpProton)
3388             {
3389             effVars[0] =  effcorection[2]->GetAxis(0)->FindBin(cent);
3390             effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx); 
3391             effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt()); 
3392             effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta()); 
3393             effvalue=effcorection[2]->GetBinContent(effVars);
3394             }
3395
3396  if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3397                 {
3398   if(parpid==SpProton || parpid==SpKaon)
3399             {
3400             effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
3401             effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx); 
3402             effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt()); 
3403             effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
3404             effvalue=effcorection[4]->GetBinContent(effVars);
3405             }
3406                 }           
3407             //    Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
3408      if(effvalue==0.) effvalue=1.;
3409
3410      return effvalue; 
3411
3412 }
3413 //---------------------------------------------------------------------------------
3414
3415 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield, Bool_t fill)
3416 {  
3417  
3418   if(!track) return 0;
3419   Bool_t trackOK = track->TestFilterBit(fFilterBit);
3420   if(!trackOK) return 0;
3421   if (fTrackStatus != 0 && !CheckTrack(track)) return 0;
3422   //select only primary traks(for data & reco MC tracks) 
3423   if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
3424   if(track->Charge()==0) return 0;
3425   if (fill) fHistQA[12]->Fill(track->GetTPCNcls());  
3426   Float_t dxy, dz;                
3427   dxy = track->DCA();
3428   dz = track->ZAtDCA();
3429   if (fill) fHistQA[6]->Fill(dxy);
3430   if (fill) fHistQA[7]->Fill(dz);
3431   Float_t chi2ndf = track->Chi2perNDF();
3432   if (fill) fHistQA[13]->Fill(chi2ndf);  
3433   // Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
3434   Float_t nCrossedRowsTPC = track->GetTPCNCrossedRows();
3435   if (fill) fHistQA[14]->Fill(nCrossedRowsTPC); 
3436   //Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
3437   if (track->GetTPCNclsF()>0) {
3438    Float_t  ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
3439    if (fill) fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
3440     }
3441 //accepted tracks  
3442      Float_t pt=track->Pt();
3443      if(pt< fminPt || pt> fmaxPt) return 0;
3444      if(TMath::Abs(track->Eta())> fmaxeta) return 0;
3445      if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
3446      //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required???
3447 // DCA XY
3448         if (fdcacut && fDCAXYCut)
3449         {
3450           if (!vertex)
3451             return 0;
3452           
3453           Double_t pos[2];
3454           Double_t covar[3];
3455           AliAODTrack* clone =(AliAODTrack*) track->Clone();
3456           Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
3457           delete clone;
3458           if (!success)
3459             return 0;
3460
3461 //        Printf("%f", ((AliAODTrack*)part)->DCA());
3462 //        Printf("%f", pos[0]);
3463           if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
3464             return 0;
3465         }
3466
3467         if (fSharedClusterCut >= 0)
3468         {
3469           Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
3470           if (frac > fSharedClusterCut)
3471             return 0;
3472         }
3473
3474    // Rejects tracks with shared clusters after filling a control histogram
3475    // This overload is used for primaries
3476
3477      // Get the shared maps
3478       const TBits sharedMap = track->GetTPCSharedMap();
3479      // Fill a control histogram
3480       fPriHistShare->Fill(sharedMap.CountBits());
3481
3482     // Reject shared clusters
3483        if (fSharedTPCmapCut >= 0)
3484         {     
3485       if((sharedMap.CountBits()) >= 1)  return 0;// Bad track, has too many shared clusters!
3486         }
3487
3488      if (fill) fHistQA[8]->Fill(pt);
3489      if (fill) fHistQA[9]->Fill(track->Eta());
3490      if (fill) fHistQA[10]->Fill(track->Phi());
3491      return 1;
3492   }
3493   //________________________________________________________________________________
3494 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos) 
3495 {
3496 //This function is called within the func GetParticle() for accepted tracks only i.e.after call of Classifytrack() & for those tracks which have proper TPC PID response . combined nsigma(circular) cut only for particles having pt upto  4.0 Gev/c and beyond that use the asymmetric nsigma cut around pion's mean position in TPC ( while filling the  TObjArray for trig & asso )
3497 Float_t pt=track->Pt();
3498
3499 //plot the separation power
3500
3501 Double_t bethe[AliPID::kSPECIES]={0.};
3502 Double_t sigma_TPC[AliPID::kSPECIES]={0.}; 
3503
3504  Double_t Pi_Ka_sep[NSigmaPIDType+1]={0.};
3505  Double_t Pi_Pr_sep[NSigmaPIDType+1]={0.};
3506  Double_t Ka_Pr_sep[NSigmaPIDType+1]={0.};
3507
3508
3509     Double_t ptpc = track->GetTPCmomentum();
3510     Int_t dEdxN = track->GetTPCsignalN();
3511  for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
3512        bethe[ipart] = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
3513       //Double_t diff = dEdx - bethe;
3514        sigma_TPC[ipart] = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
3515       //nSigma[ipart] = diff / sigma;
3516     }
3517  Pi_Ka_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kKaon])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kKaon])/2.0);
3518  Pi_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kProton])/2.0);
3519  Ka_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kKaon]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kKaon]+sigma_TPC[AliPID::kProton])/2.0);
3520
3521
3522 Double_t sigma_TOF[AliPID::kSPECIES]={0.}; 
3523
3524 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3525    {
3526  Double_t timei[AliPID::kSPECIES];
3527  track->GetIntegratedTimes(timei);
3528  for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {  sigma_TOF[ipart]= fPID->GetTOFResponse().GetExpectedSigma(track->P(), timei[ipart], AliPID::ParticleMass(ipart));}
3529  Pi_Ka_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kKaon])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kKaon])/2.0);
3530  Pi_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kProton])/2.0);
3531  Ka_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kKaon]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kKaon]+sigma_TOF[AliPID::kProton])/2.0);
3532
3533   Pi_Ka_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Ka_sep[NSigmaTPC]*Pi_Ka_sep[NSigmaTPC]+Pi_Ka_sep[NSigmaTOF]*Pi_Ka_sep[NSigmaTOF]);
3534   Pi_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Pr_sep[NSigmaTPC]*Pi_Pr_sep[NSigmaTPC]+Pi_Pr_sep[NSigmaTOF]*Pi_Pr_sep[NSigmaTOF]);
3535   Ka_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Ka_Pr_sep[NSigmaTPC]*Ka_Pr_sep[NSigmaTPC]+Ka_Pr_sep[NSigmaTOF]*Ka_Pr_sep[NSigmaTOF]);
3536    }
3537
3538
3539 //fill separation power histograms
3540  for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3541    if(ipid==0){
3542         TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3543         h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3544         TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3545         h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3546         TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3547         h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3548    }
3549    if(HasTOFPID(track) && pt>fPtTOFPIDmin && ipid!=0){
3550        TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3551         h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3552         TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3553         h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3554         TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3555         h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3556    }
3557  }
3558
3559
3560 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
3561 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
3562 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
3563 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
3564
3565 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
3566 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
3567
3568  if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3569    {
3570
3571 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
3572 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
3573 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
3574 //---combined
3575 nsigmaTPCTOFkPion   = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
3576 nsigmaTPCTOFkKaon   = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
3577 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
3578
3579
3580    }
3581 else{
3582     // --- combined
3583     // if TOF is missing and below fPtTOFPID only the TPC information is used
3584     nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
3585     nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
3586     nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
3587
3588   }
3589
3590 //set data member fnsigmas
3591   fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
3592   fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
3593   fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
3594
3595   //for all tracks below fPtTOFPIDmin  and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
3596   fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
3597   fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
3598   fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
3599
3600  //for all tracks below fPtTOFPIDmin  and also for tracks above fPtTOFPIDmin without proper TOF response these TPCTOF nsigma values will be TMath::Abs(TPC only nsigma)
3601   fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
3602   fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
3603   fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
3604
3605  if(FIllQAHistos){
3606     //Fill NSigma SeparationPlot
3607     for(Int_t ipart=0;ipart<NSpecies;ipart++){
3608       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3609         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3610         TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
3611         h->Fill(track->Pt(),fnsigmas[ipart][ipid]);
3612       }
3613     }
3614   }
3615
3616 }
3617 //----------------------------------------------------------------------------
3618 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos) 
3619 {
3620   //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3621 if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && (!HasTOFPID(track)) )return SpUndefined;//so any track having Pt>0.6 withot having proper TOF response will be defined as SpUndefined
3622 //get the identity of the particle with the minimum Nsigma
3623   Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3624   switch (fPIDType){
3625   case NSigmaTPC:
3626     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3627     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPC])  ;
3628     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPC])  ;
3629     break;
3630   case NSigmaTOF:
3631     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3632     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTOF])  ;
3633     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTOF])  ;
3634     break;
3635   case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3636     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3637     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
3638     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
3639     break;
3640   case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3641     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3642     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
3643     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
3644     break;
3645   }
3646
3647
3648 if(fdiffPIDcutvalues){
3649   if(track->Pt()<=4) fNSigmaPID=fPIDCutval1;
3650   if(track->Pt()>4 && track->Pt()<=6) fNSigmaPID=fPIDCutval2;
3651   if(track->Pt()>6 && track->Pt()<=8) fNSigmaPID=fPIDCutval3;
3652   if(track->Pt()>8) fNSigmaPID=fPIDCutval4;
3653   }
3654
3655  // guess the particle based on the smaller nsigma (within fNSigmaPID)
3656   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
3657
3658   if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
3659     if((fHighPtKaonNSigmaPID>0) && (track->Pt()>fHighPtKaonSigma) && (nsigmaKaon > fHighPtKaonNSigmaPID)) return SpUndefined;//different nsigma cut for kaons above a particular Pt range(within the TPC-TOF PID range)
3660 if(FillQAHistos){
3661       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3662         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3663         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpKaon,ipid));
3664         h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
3665       }
3666     }
3667  return SpKaon;
3668   }
3669   if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
3670  if(FillQAHistos){
3671       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3672         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3673         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpPion,ipid));
3674         h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
3675       }
3676     }
3677 return SpPion;
3678   }
3679   if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
3680 if(FillQAHistos){
3681       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3682         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3683         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpProton,ipid));
3684         h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
3685       }
3686  }
3687 return SpProton;
3688   }
3689
3690 // else, return undefined
3691   return SpUndefined;
3692   
3693  
3694 }
3695
3696 //------------------------------------------------------------------------------------------
3697 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){ 
3698   //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3699
3700   //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
3701   //fill DC histos
3702   for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
3703   
3704   Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
3705   
3706   
3707   if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
3708   
3709   Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3710   switch (fPIDType) {
3711   case NSigmaTPC:
3712     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3713     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPC])  ;
3714     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPC])  ;
3715     break;
3716   case NSigmaTOF:
3717     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3718     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTOF])  ;
3719     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTOF])  ;
3720     break;
3721   case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3722     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3723     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
3724     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
3725     break;
3726   case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3727     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3728     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
3729     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
3730     break;
3731   }
3732
3733   // Actually the tracks in the overlapping region(in TPC-TOF nSigma plane) will be ignored
3734
3735   if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
3736   if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
3737   if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
3738      
3739   
3740
3741 if(FIllQAHistos){
3742     //fill NSigma distr for double counting
3743     for(Int_t ipart=0;ipart<NSpecies;ipart++){
3744       if(fHasDoubleCounting[ipart]){//this may be kTRUE only for particles having Pt<=4.0 GeV/C, so this histo contains all the particles having Pt<=4.0 GeV/C in the nsigma overlapping region in TPC/TPC-TOF plane 
3745         for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3746           if((ipid!=NSigmaTPC) && (!HasTOFPID(trk)))continue;//not filling TOF and combined if no TOF PID
3747           TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
3748           h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
3749         }
3750       }
3751     }
3752   }
3753  
3754  
3755   return fHasDoubleCounting;
3756 }
3757
3758 //////////////////////////////////////////////////////////////////////////////////////////////////
3759
3760 Bool_t* AliTwoParticlePIDCorr::GetAllCompatibleIdentitiesNSigma(AliAODTrack * trk,Bool_t FIllQAHistos){ 
3761  //mainly intended to check the probability of the PID of the tracks which are in the overlapping nSigma regions and near about the middle position from the   mean position of two ID particle
3762   Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
3763   IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
3764   return IDs;
3765   
3766 }
3767 //////////////////////////////////////////////////////////////////////////////////////////////////
3768
3769 UInt_t AliTwoParticlePIDCorr::CalcPIDCombined(AliAODTrack *track, Int_t detMask, Double_t* prob) const{
3770   //
3771   // Bayesian PID calculation
3772   //
3773   for(Int_t i=0;i<AliPID::kSPECIES;i++)
3774     {
3775       prob[i]=0.;
3776     }
3777   fPIDCombined->SetDetectorMask(detMask);
3778   
3779   return fPIDCombined->ComputeProbabilities((AliAODTrack*)track, fPID, prob);
3780 }
3781
3782 //////////////////////////////////////////////////////////////////////////////////////////////////
3783
3784 Int_t AliTwoParticlePIDCorr::GetIDBayes(AliAODTrack * trk, Bool_t FIllQAHistos){ 
3785   
3786   Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
3787
3788
3789   //Filling of Probability histos
3790         Double_t probTPC[AliPID::kSPECIES]={0.};
3791         Double_t probTOF[AliPID::kSPECIES]={0.};
3792         Double_t probTPCTOF[AliPID::kSPECIES]={0.};
3793
3794         UInt_t detUsedTPC = 0;
3795         UInt_t detUsedTOF = 0;
3796         UInt_t detUsedTPCTOF = 0;
3797
3798  //get the TPC probability
3799           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
3800           detUsedTPC = fPIDCombined->ComputeProbabilities(trk, fPID, probTPC);
3801 if(detUsedTPC == AliPIDResponse::kDetTPC)
3802   {
3803 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3804
3805         TH2F *h=GetHistogram2D(Form("probBayes_TPC_%d",ipart));
3806         if(ipart==0)    h->Fill(trk->Pt(),probTPC[AliPID::kPion]);
3807         if(ipart==1)    h->Fill(trk->Pt(),probTPC[AliPID::kKaon]);
3808         if(ipart==2)    h->Fill(trk->Pt(),probTPC[AliPID::kProton]);
3809  }
3810   }
3811
3812   //get the TOF probability
3813           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
3814           detUsedTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTOF);
3815 if(detUsedTOF == AliPIDResponse::kDetTOF)
3816   {
3817 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3818         TH2F *h=GetHistogram2D(Form("probBayes_TOF_%d",ipart));
3819         if(ipart==0)    h->Fill(trk->Pt(),probTOF[AliPID::kPion]);
3820         if(ipart==1)    h->Fill(trk->Pt(),probTOF[AliPID::kKaon]);
3821         if(ipart==2)    h->Fill(trk->Pt(),probTOF[AliPID::kProton]);
3822  }
3823   }
3824
3825  //get the TPC-TOF probability
3826           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
3827           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTPCTOF);
3828 if(detUsedTPCTOF == (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))
3829   {
3830 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3831         TH2F *h=GetHistogram2D(Form("probBayes_TPCTOF_%d",ipart));
3832         if(ipart==0)    h->Fill(trk->Pt(),probTPCTOF[AliPID::kPion]);
3833         if(ipart==1)    h->Fill(trk->Pt(),probTPCTOF[AliPID::kKaon]);
3834         if(ipart==2)    h->Fill(trk->Pt(),probTPCTOF[AliPID::kProton]); 
3835 }
3836   }
3837
3838   
3839   Double_t probBayes[AliPID::kSPECIES];
3840   
3841   UInt_t detUsed= 0;
3842   if(HasTOFPID(trk) && trk->Pt()>fPtTOFPIDmin){//use TOF information
3843     detUsed = CalcPIDCombined(trk, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
3844     if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return SpUndefined;//check that TPC and TOF are used
3845   }else{
3846     detUsed = CalcPIDCombined(trk,AliPIDResponse::kDetTPC, probBayes);
3847     if(detUsed != AliPIDResponse::kDetTPC)return SpUndefined;//check that TPC is used
3848   }
3849   
3850   //the probability has to be normalized to one, we check it
3851   Double_t sump=0.;
3852   for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
3853   if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
3854     AliFatal("Bayesian probability not normalized to one");
3855   }
3856
3857   if(fdiffPIDcutvalues){
3858   if(trk->Pt()<=4) fBayesCut=fPIDCutval1;
3859   if(trk->Pt()>4 && trk->Pt()<=6) fBayesCut=fPIDCutval2;
3860   if(trk->Pt()>6 && trk->Pt()<=8) fBayesCut=fPIDCutval3;
3861   if(trk->Pt()>8) fBayesCut=fPIDCutval4;
3862   }
3863
3864   
3865   //probabilities are normalized to one, if the cut is above .5 there is no problem
3866   if(probBayes[AliPID::kPion]>fBayesCut && IDs[SpPion]==1){
3867     TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpPion));
3868     h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
3869     return SpPion;
3870   }
3871   else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[SpKaon]==1){
3872     TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpKaon));
3873     h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
3874     return SpKaon;
3875   }
3876   else if(probBayes[AliPID::kProton]>fBayesCut && IDs[SpProton]==1){
3877     TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpProton));
3878     h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
3879     return SpProton;
3880   }
3881   else{
3882     return SpUndefined;
3883   }
3884 }
3885
3886
3887 //////////////////////////////////////////////////////////////////////////////////////////////////
3888 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){ 
3889   //return the specie according to the minimum nsigma value
3890   //no double counting, this has to be evaluated using CheckDoubleCounting()
3891   
3892   Int_t ID=SpUndefined;
3893
3894   CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
3895
3896
3897  //Do PID
3898   if(fPIDType==Bayes){//use bayesianPID
3899     
3900     if(!fPIDCombined) {
3901       AliFatal("PIDCombined object has to be set in the steering macro");
3902     }
3903     
3904     ID = GetIDBayes(trk,FIllQAHistos);
3905     
3906   }else{ //use nsigma PID  
3907
3908    ID=FindMinNSigma(trk,FIllQAHistos);
3909 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
3910       Bool_t *HasDC;
3911       HasDC=GetDoubleCounting(trk,FIllQAHistos);
3912       for(Int_t ipart=0;ipart<NSpecies;ipart++){
3913         if(HasDC[ipart]==kTRUE)  ID = SpUndefined;
3914       }
3915     }
3916   }
3917 //Fill PID signal plot
3918   if(ID != SpUndefined){
3919     for(Int_t idet=0;idet<fNDetectors;idet++){
3920       TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
3921       if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3922       if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3923       if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3924     }
3925   }
3926   //Fill PID signal plot without cuts
3927   for(Int_t idet=0;idet<fNDetectors;idet++){
3928     TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
3929     if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3930     if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3931     if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3932   }
3933   
3934   return ID;
3935 }
3936
3937 //-------------------------------------------------------------------------------------
3938 Bool_t
3939 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
3940 {  
3941   // check PID signal 
3942    AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
3943    if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
3944    //ULong_t status=track->GetStatus();
3945    //if  (!( (status & AliAODTrack::kTPCpid  ) == AliAODTrack::kTPCpid  )) return kFALSE;//remove light nuclei
3946    //if (track->GetTPCsignal() <= 0.) return kFALSE;
3947    // if(track->GetTPCsignalN() < 60) return kFALSE;//tracks with TPCsignalN< 60 have questionable dEdx,cutting on TPCsignalN > 70 or > 60 shouldn't make too much difference in statistics,also  it is IMO safe to use TPC also for MIPs.
3948    
3949   return kTRUE;  
3950 }
3951 //___________________________________________________________
3952
3953 Bool_t
3954 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
3955 {
3956   // check TOF matched track 
3957   //ULong_t status=track->GetStatus();
3958   //if  (!( (status & AliAODTrack::kITSin  ) == AliAODTrack::kITSin  )) return kFALSE;
3959  AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
3960  if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
3961   if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
3962  //if(!((status & AliAODTrack::kTOFpid  ) == AliAODTrack::kTOFpid  )) return kFALSE;
3963  //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
3964  // if (probMis > 0.01) return kFALSE;
3965 if(fRemoveTracksT0Fill)
3966     {
3967 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
3968       if (startTimeMask < 0)return kFALSE; 
3969     }
3970   return kTRUE;
3971 }
3972
3973 //________________________________________________________________________
3974 Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
3975 {
3976   //it is called only when TOF PID is available
3977 //TOF beta calculation
3978   Double_t tofTime=track->GetTOFsignal();
3979   
3980   Double_t c=TMath::C()*1.E-9;// m/ns
3981   Float_t startTime = fPID->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
3982   Double_t length= fPID->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
3983   tofTime -= startTime;      // subtract startTime to the signal
3984   Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector         
3985   tof=tof*c;
3986   return length/tof;
3987
3988
3989   /*
3990   Double_t p = track->P();
3991   Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
3992   Double_t timei[5];
3993   track->GetIntegratedTimes(timei);
3994   return timei[0]/time;
3995   */
3996 }
3997 //------------------------------------------------------------------------------------------------------
3998
3999 Float_t AliTwoParticlePIDCorr::GetInvMassSquared(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
4000 {
4001   // calculate inv mass squared
4002   // same can be achieved, but with more computing time with
4003   /*TLorentzVector photon, p1, p2;
4004   p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
4005   p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
4006   photon = p1+p2;
4007   photon.M()*/
4008   
4009   Float_t tantheta1 = 1e10;
4010   
4011   if (eta1 < -1e-10 || eta1 > 1e-10)
4012     tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
4013   
4014   Float_t tantheta2 = 1e10;
4015   if (eta2 < -1e-10 || eta2 > 1e-10)
4016     tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
4017   
4018   Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4019   Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4020   
4021   Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1 - phi2) + 1.0 / tantheta1 / tantheta2 ) ) );
4022   
4023   return mass2;
4024 }
4025 //---------------------------------------------------------------------------------
4026
4027 Float_t AliTwoParticlePIDCorr::GetInvMassSquaredCheap(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
4028 {
4029   // calculate inv mass squared approximately
4030   
4031   Float_t tantheta1 = 1e10;
4032   
4033   if (eta1 < -1e-10 || eta1 > 1e-10)
4034   {
4035     Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
4036     tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4037   }
4038   
4039   Float_t tantheta2 = 1e10;
4040   if (eta2 < -1e-10 || eta2 > 1e-10)
4041   {
4042     Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
4043     tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4044   }
4045   
4046   Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4047   Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4048   
4049   // fold onto 0...pi
4050   Float_t deltaPhi = TMath::Abs(phi1 - phi2);
4051   while (deltaPhi > TMath::TwoPi())
4052     deltaPhi -= TMath::TwoPi();
4053   if (deltaPhi > TMath::Pi())
4054     deltaPhi = TMath::TwoPi() - deltaPhi;
4055   
4056   Float_t cosDeltaPhi = 0;
4057   if (deltaPhi < TMath::Pi()/3)
4058     cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
4059   else if (deltaPhi < 2*TMath::Pi()/3)
4060     cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
4061   else
4062     cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
4063   
4064   Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
4065   
4066 //   Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
4067   
4068   return mass2;
4069 }
4070 //--------------------------------------------------------------------------------
4071 Float_t  AliTwoParticlePIDCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign)
4072
4073   //
4074   // calculates dphistar
4075   //
4076   
4077   Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
4078   
4079   static const Double_t kPi = TMath::Pi();
4080   
4081   // circularity
4082 //   if (dphistar > 2 * kPi)
4083 //     dphistar -= 2 * kPi;
4084 //   if (dphistar < -2 * kPi)
4085 //     dphistar += 2 * kPi;
4086   
4087   if (dphistar > kPi)
4088     dphistar = kPi * 2 - dphistar;
4089   if (dphistar < -kPi)
4090     dphistar = -kPi * 2 - dphistar;
4091   if (dphistar > kPi) // might look funny but is needed
4092     dphistar = kPi * 2 - dphistar;
4093   
4094   return dphistar;
4095 }
4096 //_________________________________________________________________________
4097 /*
4098 void AliTwoParticlePIDCorr ::DefineEventPool()
4099 {
4100 Int_t MaxNofEvents=1000;
4101 const Int_t NofVrtxBins=10+(1+10)*2;
4102 Double_t ZvrtxBins[NofVrtxBins+1]={ -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10, 
4103                                        90,  92,  94,  96,  98, 100, 102, 104, 106, 108, 110, 
4104                                       190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210 
4105
4106 //default values are for centrality
4107 Int_t  NofCentBins=15;
4108 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
4109
4110  if(fCentralityMethod.EndsWith("_MANUAL"))
4111    {
4112  Int_t  NofCentBins=9;
4113  CentralityBins[NofCentBins+1]={0.,9.,14.,19.,26.,34.,44.,58.,80.,500.};//Is This binning is fine for pp, or we don't require them....
4114    }
4115 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
4116
4117
4118
4119   
4120 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
4121
4122 //if(!fPoolMgr) return kFALSE;
4123 //return kTRUE;
4124
4125 }
4126 */
4127 //------------------------------------------------------------------------
4128 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
4129 {
4130   // This method is a copy from AliUEHist::GetBinning
4131   // takes the binning from <configuration> identified by <tag>
4132   // configuration syntax example:
4133   // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
4134   // phi: .....
4135   //
4136   // returns bin edges which have to be deleted by the caller
4137   
4138   TString config(configuration);
4139   TObjArray* lines = config.Tokenize("\n");
4140   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
4141   {
4142     TString line(lines->At(i)->GetName());
4143     if (line.BeginsWith(TString(tag) + ":"))
4144     {
4145       line.Remove(0, strlen(tag) + 1);
4146       line.ReplaceAll(" ", "");
4147       TObjArray* binning = line.Tokenize(",");
4148       Double_t* bins = new Double_t[binning->GetEntriesFast()];
4149       for (Int_t j=0; j<binning->GetEntriesFast(); j++)
4150         bins[j] = TString(binning->At(j)->GetName()).Atof();
4151       
4152       nBins = binning->GetEntriesFast() - 1;
4153
4154       delete binning;
4155       delete lines;
4156       return bins;
4157     }
4158   }
4159   
4160   delete lines;
4161   AliFatal(Form("Tag %s not found in %s", tag, configuration));
4162   return 0;
4163 }
4164
4165 //____________________________________________________________________
4166
4167 Bool_t AliTwoParticlePIDCorr::CheckTrack(AliAODTrack * part)
4168 {
4169   // check if the track status flags are set
4170   
4171   UInt_t status=((AliAODTrack*)part)->GetStatus();
4172   if ((status & fTrackStatus) == fTrackStatus)
4173     return kTRUE;
4174   return kFALSE;
4175 }
4176 //________________________________________________________________________
4177
4178 Bool_t AliTwoParticlePIDCorr::AcceptEventCentralityWeight(Double_t centrality)
4179 {
4180   // rejects "randomly" events such that the centrality gets flat
4181   // uses fCentralityWeights histogram
4182
4183   // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
4184   
4185   Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
4186   Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
4187   
4188   Bool_t result = kFALSE;
4189   if (centralityDigits < weight) 
4190     result = kTRUE;
4191   
4192   AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
4193   
4194   return result;
4195 }
4196
4197 //____________________________________________________________________
4198 void AliTwoParticlePIDCorr::ShiftTracks(TObjArray* tracks, Double_t angle)
4199 {
4200   // shifts the phi angle of all tracks by angle
4201   // 0 <= angle <= 2pi
4202   
4203   for (Int_t i=0; i<tracks->GetEntriesFast(); ++i) 
4204   {
4205    LRCParticlePID *part=(LRCParticlePID*)(tracks->UncheckedAt(i));
4206
4207     Double_t newAngle = part->Phi() + angle; 
4208     if (newAngle >= TMath::TwoPi())
4209       newAngle -= TMath::TwoPi();
4210     
4211     part->SetPhi(newAngle);
4212   }
4213 }
4214
4215
4216 //________________________________________________________________________
4217 void  AliTwoParticlePIDCorr::SetVZEROCalibrationFile(const char* filename,const char* lhcPeriod) {
4218   //Function to setup the VZERO gain equalization
4219     //============Get the equilization map============//
4220   TFile *calibrationFile = TFile::Open(filename);
4221   if((!calibrationFile)||(!calibrationFile->IsOpen())) {
4222     Printf("No calibration file found!!!");
4223     return;
4224   }
4225
4226   TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
4227   if(!list) {
4228     Printf("Calibration TList not found!!!");
4229     return;
4230   }
4231
4232   fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
4233   if(!fHistVZEROAGainEqualizationMap) {
4234     Printf("VZERO-A calibration object not found!!!");
4235     return;
4236   }
4237   fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
4238   if(!fHistVZEROCGainEqualizationMap) {
4239     Printf("VZERO-C calibration object not found!!!");
4240     return;
4241   }
4242
4243   fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
4244   if(!fHistVZEROChannelGainEqualizationMap) {
4245     Printf("VZERO channel calibration object not found!!!");
4246     return;
4247   }
4248 }
4249
4250 //________________________________________________________________________
4251 Double_t AliTwoParticlePIDCorr::GetChannelEqualizationFactor(Int_t run,Int_t channel) {
4252   //
4253   if(!fHistVZEROAGainEqualizationMap) return 1.0;
4254
4255   for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
4256     Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4257     if(gRunNumber == run)
4258       return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
4259   }
4260
4261   return 1.0;
4262 }
4263
4264 //________________________________________________________________________
4265 Double_t AliTwoParticlePIDCorr::GetEqualizationFactor(Int_t run, const char* side) {
4266   //
4267   if(!fHistVZEROAGainEqualizationMap) return 1.0;
4268
4269   TString gVZEROSide = side;
4270   for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
4271     Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4272     //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
4273     if(gRunNumber == run) {
4274       if(gVZEROSide == "A") 
4275         return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
4276       else if(gVZEROSide == "C") 
4277         return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
4278     }
4279   }
4280
4281   return 1.0;
4282 }
4283 //________________________________________________________________________
4284 Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliAODEvent *event){
4285   //Function that returns the reference multiplicity from AODs (data or reco MC, Not for Truth)
4286   //Different ref. mult. implemented: V0M, V0A, V0C, TPC
4287   Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
4288   Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
4289
4290   AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
4291   if(!header) {
4292     Printf("ERROR: AOD header not available");
4293     return -999;
4294   }
4295   Int_t gRunNumber = header->GetRunNumber();
4296  Float_t bSign1=header->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
4297
4298
4299  for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) 
4300 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation 
4301   AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
4302   if (!track) continue;
4303   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
4304   if(tracktype!=1) continue; 
4305
4306   if(!track) continue;//for safety
4307
4308     gRefMultiplicityTPC += 1.0;
4309
4310  }//track looop ends
4311
4312  if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
4313   //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
4314   for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
4315     fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
4316     
4317     if(iChannel < 32) 
4318       gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
4319     else if(iChannel >= 32) 
4320       gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
4321   }//loop over PMTs
4322   
4323   //Equalization of gain
4324   Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
4325   if(gFactorA != 0)
4326     gRefMultiplicityVZEROA /= gFactorA;
4327   Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
4328   if(gFactorC != 0)
4329     gRefMultiplicityVZEROC /= gFactorC;
4330   if((gFactorA != 0)&&(gFactorC != 0)) 
4331     gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
4332
4333       
4334   //EQVZERO vs TPC multiplicity
4335   fHistEQVZEROvsTPCmultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
4336   fHistEQVZEROAvsTPCmultiplicity->Fill(gRefMultiplicityVZEROA,gRefMultiplicityTPC);
4337   fHistEQVZEROCvsTPCmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityTPC);
4338
4339   //EQVZERO vs VZERO multiplicity
4340   fHistVZEROCvsEQVZEROCmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),gRefMultiplicityVZEROC);
4341   fHistVZEROAvsEQVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0A(),gRefMultiplicityVZEROA);
4342
4343   //VZEROC vs VZEROA multiplicity
4344   fHistVZEROCvsVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),event->GetVZEROData()->GetMTotV0A());
4345
4346   //EQVZEROC vs EQVZEROA multiplicity
4347   fHistEQVZEROCvsEQVZEROAmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityVZEROA);
4348  }
4349
4350  if(fCentralityMethod == "TRACKS_MANUAL") {
4351     gRefMultiplicity = gRefMultiplicityTPC;
4352     fHistRefmult->Fill(3.,gRefMultiplicityTPC);
4353  }
4354  else if(fCentralityMethod == "V0M_MANUAL"){
4355     gRefMultiplicity = gRefMultiplicityVZERO;
4356     fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
4357
4358  }
4359  else if(fCentralityMethod == "V0A_MANUAL"){
4360     gRefMultiplicity = gRefMultiplicityVZEROA;
4361     fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
4362
4363  }
4364  else if(fCentralityMethod == "V0C_MANUAL"){
4365     gRefMultiplicity = gRefMultiplicityVZEROC;
4366     fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
4367  }
4368  else {
4369 gRefMultiplicity = gRefMultiplicityTPC;
4370  } 
4371   return gRefMultiplicity;
4372 }
4373
4374 //-------------------------------------------------------------------------------------------------------
4375 Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliAODEvent *event, Bool_t truth){
4376
4377   if(!event) return -1;
4378   // get centrality object and check quality
4379   Double_t cent_v0=-1;
4380   Double_t nooftrackstruth=0;
4381   Bool_t shift_to_TRACKS_MANUAL=kFALSE;//in case of wrong setting automatic shift to Tracks_Manual method
4382
4383 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")//for PbPb, pPb, pp7TeV(still to be introduced)//data or RecoMC and also for TRUTH
4384     {
4385       /*
4386 if(fSampleType=="pp_7" && fPPVsMultUtils)
4387 {//for pp 7 TeV case only using Alianalysisutils class
4388         if(fAnalysisUtils) cent_v0 = fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,fCentralityMethod);
4389         else cent_v0 = -1;
4390   fHistCentStats->Fill(0.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0A"));
4391   fHistCentStats->Fill(1.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0C"));
4392   fHistCentStats->Fill(2.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0M"));   
4393   fHistCentStats->Fill(3.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0AEq"));//only available for LHC10d at present (Quantile info)
4394   fHistCentStats->Fill(4.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0CEq"));//only available for LHC10d at present (Quantile info)
4395   fHistCentStats->Fill(5.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0MEq"));//only available for LHC10d at present (Quantile info)
4396       }
4397       */
4398   if(fSampleType=="pPb" || fSampleType=="PbPb")
4399   {
4400   AliCentrality *centralityObj=0;
4401   AliAODHeader *header = (AliAODHeader*) event->GetHeader();
4402   if(!header) return -1;
4403   centralityObj = header->GetCentralityP();
4404   // if (centrality->GetQuality() != 0) return ;
4405   if(centralityObj){
4406   fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0A"));
4407   fHistCentStats->Fill(1.,centralityObj->GetCentralityPercentile("V0C"));
4408   fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0M"));
4409   fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
4410   fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
4411   fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
4412
4413   fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
4414   fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA")); 
4415
4416    cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
4417     }
4418   else cent_v0= -1;    
4419   }
4420  else shift_to_TRACKS_MANUAL=kTRUE;    
4421
4422     }//centralitymethod condition
4423
4424  else if(fCentralityMethod=="V0M_MANUAL" || fCentralityMethod=="V0A_MANUAL" || fCentralityMethod=="V0C_MANUAL" || fCentralityMethod=="TRACKS_MANUAL" || shift_to_TRACKS_MANUAL)//data or RecoMc and also for TRUTH
4425    {
4426      if(!truth){//for data or RecoMC
4427     cent_v0 = GetReferenceMultiplicityVZEROFromAOD(event);
4428    }//for data or RecoMC
4429
4430     if(truth && (fAnalysisType == "MCAOD")){//condition for TRUTH case
4431 //check for TClonesArray(truth track MC information)
4432 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
4433   if (!fArrayMC) {
4434     //AliFatal("Error: MC particles branch not found!\n");
4435     return -1;
4436   }
4437 //now process the truth particles(for both efficiency & correlation function)
4438 Int_t nMCTrack = fArrayMC->GetEntriesFast();
4439   
4440 for (Int_t iMC = 0; iMC < nMCTrack; iMC++) 
4441 {//MC truth track loop starts
4442     
4443 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
4444     
4445     if(!partMC){
4446       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
4447       continue;
4448     }
4449
4450 //consider only charged particles
4451     if(partMC->Charge() == 0) continue;
4452
4453 //consider only primary particles; neglect all secondary particles including from weak decays
4454  if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
4455
4456
4457 //remove injected signals(primaries above <maxLabel>)
4458  if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
4459
4460 //remove duplicates
4461   Bool_t isduplicate=kFALSE;
4462  if (fRemoveDuplicates)
4463    { 
4464  for (Int_t j=iMC+1; j<nMCTrack; ++j) 
4465    {//2nd trutuh loop starts
4466 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
4467    if(!partMC2){
4468       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
4469       continue;
4470     }    
4471  if (partMC->GetLabel() == partMC2->GetLabel())
4472    {
4473 isduplicate=kTRUE;
4474  break;  
4475    }    
4476    }//2nd truth loop ends
4477    }
4478  if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
4479
4480
4481       if (fCentralityMethod=="V0M_MANUAL") {
4482         if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8)    continue;
4483         if (partMC->Eta() < -3.7 || partMC->Eta() > -1.7) continue;
4484 }
4485       else if (fCentralityMethod=="V0A_MANUAL") {
4486         if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8)  continue;}
4487       else if (fCentralityMethod=="V0C_MANUAL") {
4488         if(partMC->Eta() > -1.7 || partMC->Eta() < -3.7)  continue;}
4489       else if (fCentralityMethod=="TRACKS_MANUAL") {
4490         if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4491         if (partMC->Pt() < fminPt ||  partMC->Pt() > fmaxPt) continue;
4492            }
4493       else{//basically returns the tracks manual case
4494 //give only kinematic cuts at the truth level  
4495        if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4496        if (partMC->Pt() < fminPt ||  partMC->Pt() > fmaxPt) continue;
4497       }
4498
4499  //To determine multiplicity in case of PP
4500  nooftrackstruth+= 1;;
4501
4502  }//truth track loop ends
4503  cent_v0=nooftrackstruth;
4504
4505     }//condition for TRUTH case
4506
4507    }//end of MANUAL method
4508
4509  else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC
4510     {
4511     AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4512     if (!header)
4513     return -1;
4514     
4515       AliGenEventHeader* eventHeader = header->GetCocktailHeader(0);  // get first MC header from either ESD/AOD (including cocktail header if available)
4516       if (!eventHeader)
4517       {
4518         // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing 
4519         // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
4520         AliError("Event header not found. Skipping this event.");
4521         return -1;
4522       }
4523       
4524       AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4525      
4526       
4527      if (collGeometry)   cent_v0 = collGeometry->ImpactParameter();
4528       else cent_v0=-1.;
4529     }//end of Impact parameter method
4530
4531 //else return -1
4532  else cent_v0=-1.;
4533
4534  return cent_v0;
4535 }
4536 //-----------------------------------------------------------------------------------------
4537 Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliAODEvent *aod,Bool_t truth){
4538   //do the event selection(zvtx, pileup, centrality/multiplicity cut) and then return the value of the centrality of that event
4539   if(!aod) return -1;
4540
4541   Float_t gRefMultiplicity = -1.;
4542
4543   // check first event in chunk (is not needed for new reconstructions)
4544   if(fCheckFirstEventInChunk){
4545     AliAnalysisUtils ut;
4546     if(ut.IsFirstEventInChunk(aod)) 
4547       return -1.;
4548   }
4549
4550  if(frejectPileUp){
4551     AliAnalysisUtils ut;
4552     ut.SetUseMVPlpSelection(kTRUE);
4553     ut.SetUseOutOfBunchPileUp(kTRUE);
4554     if(ut.IsPileUpEvent(aod))
4555       return -1.;
4556   }
4557
4558 //count events after pileup selection
4559    fEventCounter->Fill(3);
4560
4561   //vertex selection(is it fine for PP?)
4562  if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return; 
4563   trkVtx = aod->GetPrimaryVertex();
4564   if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;
4565   TString vtxTtl = trkVtx->GetTitle();
4566   if (!vtxTtl.Contains("VertexerTracks")) return -1;
4567    zvtx = trkVtx->GetZ();
4568   const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
4569   if (!spdVtx || spdVtx->GetNContributors()<=0) return -1;
4570   TString vtxTyp = spdVtx->GetTitle();
4571   Double_t cov[6]={0};
4572   spdVtx->GetCovarianceMatrix(cov);
4573   Double_t zRes = TMath::Sqrt(cov[5]);
4574   if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return -1;
4575    if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return -1;
4576   }
4577   else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
4578         Int_t nVertex = aod->GetNumberOfVertices();
4579         if( nVertex > 0 ) { 
4580      trkVtx = aod->GetPrimaryVertex();
4581                 Int_t nTracksPrim = trkVtx->GetNContributors();
4582                  zvtx = trkVtx->GetZ();
4583                 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
4584                 // Reject TPC only vertex
4585                 TString name(trkVtx->GetName());
4586                 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return -1;
4587
4588                 // Select a quality vertex by number of tracks?
4589                 if( nTracksPrim < fnTracksVertex ) {
4590                   //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
4591                         return -1;
4592                         }
4593                 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
4594                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
4595                 //  return kFALSE;
4596                 //      if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
4597         }
4598         else return -1;
4599
4600   }
4601  else if(fVertextype==0){//default case
4602   trkVtx = aod->GetPrimaryVertex();
4603   if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;//proper number of contributors
4604   zvtx = trkVtx->GetZ();
4605   Double32_t fCov[6];
4606   trkVtx->GetCovarianceMatrix(fCov);
4607   if(fCov[5] == 0) return -1;//proper vertex resolution
4608   }
4609   else {
4610    AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
4611    return -1;//as there is no proper sample type
4612   }
4613
4614 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
4615
4616 //count events having a proper vertex
4617    fEventCounter->Fill(5);
4618
4619  if (TMath::Abs(zvtx) > fzvtxcut) return -1;
4620
4621 //count events after vertex cut
4622   fEventCounter->Fill(7);
4623
4624
4625  //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4626   
4627  fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
4628
4629  //get the centrality or multiplicity
4630  if(truth)  {gRefMultiplicity = GetRefMultiOrCentrality(aod,kTRUE);}//kTRUE-->for Truth case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
4631
4632  else {gRefMultiplicity = GetRefMultiOrCentrality(aod,kFALSE);}//kFALSE-->for data and RecoMc case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
4633
4634   if(gRefMultiplicity<0) return -1;
4635
4636  // take events only within the  multiplicity class mentioned in the custom binning
4637   if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
4638
4639 //count events having proper centrality/ref multiplicity
4640   fEventCounter->Fill(9);
4641
4642
4643 // centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
4644  if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
4645   {
4646     AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
4647     return -1;
4648   }
4649
4650 //count events after rejection due to centrality weighting
4651   fEventCounter->Fill(11);
4652
4653   return gRefMultiplicity;
4654
4655 }
4656 //--------------------------------------------------------------------------------------------------------
4657 Float_t AliTwoParticlePIDCorr::GetEventPlane(AliAODEvent *event,Bool_t truth, Double_t v0Centr)
4658 {
4659   // Get the event plane
4660  //reset Q vector info  
4661
4662     Int_t run = event->GetRunNumber();
4663
4664     if(run != fRun){
4665         // Load the calibrations run dependent
4666       if(! fIsAfter2011) OpenInfoCalbration(run);
4667       fRun=run;
4668     }
4669
4670
4671   Int_t iC = -1;  
4672 if (v0Centr < 80){ // analysis only for 0-80% centrality classes
4673  // centrality bins
4674     if(v0Centr < 5) iC = 0;
4675     else if(v0Centr < 10) iC = 1;
4676     else if(v0Centr < 20) iC = 2;
4677     else if(v0Centr < 30) iC = 3;
4678     else if(v0Centr < 40) iC = 4;
4679     else if(v0Centr < 50) iC = 5;
4680     else if(v0Centr < 60) iC = 6;
4681     else if(v0Centr < 70) iC = 7;
4682     else iC = 8;
4683
4684
4685      Int_t iCcal = iC;
4686
4687  //reset Q vector info   
4688     Double_t Qxa2 = 0, Qya2 = 0;
4689     Double_t Qxc2 = 0, Qyc2 = 0;
4690     Double_t Qxa3 = 0, Qya3 = 0;
4691     Double_t Qxc3 = 0, Qyc3 = 0;
4692
4693
4694   //MC: from reaction plane
4695  if(truth)
4696 {
4697     AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4698     if (header){
4699       evplaneMC = header->GetReactionPlaneAngle();//[0, 360]
4700         //make it within [-pi/2,pi/2] to make it general
4701         if(evplaneMC > TMath::Pi()/2 && evplaneMC <=  TMath::Pi()*3/2) evplaneMC-=TMath::Pi(); 
4702         else if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
4703          fHistEventPlaneTruth->Fill(iC,evplaneMC);
4704         /*
4705         AliGenEventHeader* eventHeader = header->GetCocktailHeader(0);  // get first MC header from either ESD/AOD (including cocktail header if available)
4706       if (eventHeader)
4707       {
4708               
4709         AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);     
4710       
4711         if (collGeometry){//get the reaction plane from MC header   
4712           gReactionPlane = collGeometry->ReactionPlaneAngle();//[0,180]
4713  }
4714       }
4715         */   
4716      //taken from vnv0 code(get the TPC, V0A, V0C event plane using truth tracks)
4717        TClonesArray *mcArray = NULL;
4718         mcArray = (TClonesArray*)event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
4719         if(mcArray){
4720           Float_t QxMCv2[3] = {0,0,0};
4721           Float_t QyMCv2[3] = {0,0,0};
4722           Float_t QxMCv3[3] = {0,0,0};
4723           Float_t QyMCv3[3] = {0,0,0};
4724           Float_t EvPlaneMCV2[3] = {0,0,0};
4725           Float_t EvPlaneMCV3[3] = {0,0,0};
4726           Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
4727           Float_t etaMax[3] = {4.88,-1.8,0.8};
4728
4729           // analysis on MC tracks
4730           Int_t nMCtrack = mcArray->GetEntries() ;
4731
4732           // EP computation with MC tracks
4733           for(Int_t iT=0;iT < nMCtrack;iT++){
4734             AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
4735             if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
4736             
4737             Float_t eta = mctr->Eta();
4738   for(Int_t iD=0;iD<3;iD++){
4739               if(eta > etaMin[iD] && eta < etaMax[iD]){
4740                 Float_t phi = mctr->Phi();
4741                 QxMCv2[iD] += TMath::Cos(2*phi);
4742                 QyMCv2[iD] += TMath::Sin(2*phi);
4743                 QxMCv3[iD] += TMath::Cos(3*phi);
4744                 QyMCv3[iD] += TMath::Sin(3*phi);
4745               }
4746             }
4747           }
4748
4749             EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
4750             EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
4751             EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
4752             fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
4753             fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
4754             fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
4755             fgPsi2v0aMC = EvPlaneMCV2[0];
4756             fgPsi2v0cMC = EvPlaneMCV2[1];
4757             fgPsi2tpcMC = EvPlaneMCV2[2];
4758           
4759
4760             EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
4761             EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
4762             EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
4763             fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
4764             fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
4765             fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
4766             fgPsi3v0aMC = EvPlaneMCV3[0];
4767             fgPsi3v0cMC = EvPlaneMCV3[1];
4768             fgPsi3tpcMC = EvPlaneMCV3[2];
4769           
4770         }    
4771     }
4772
4773     }
4774  else{
4775     Int_t nAODTracks = event->GetNumberOfTracks();
4776
4777 // TPC EP needed for resolution studies (TPC subevent)
4778    //AliEventplane * ep = (fAOD->GetHeader())->GetEventplaneP();
4779    //Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
4780     Double_t Qx2 = 0, Qy2 = 0;
4781     Double_t Qx3 = 0, Qy3 = 0;
4782
4783     for(Int_t iT = 0; iT < nAODTracks; iT++) {
4784       
4785       AliAODTrack* aodTrack = event->GetTrack(iT);
4786       
4787       if (!aodTrack){
4788         continue;
4789       }
4790       
4791       Bool_t trkFlag = aodTrack->TestFilterBit(1);
4792
4793       if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag) 
4794         continue;
4795         
4796       Double_t b[2] = {-99., -99.};
4797       Double_t bCov[3] = {-99., -99., -99.};
4798
4799
4800       AliAODTrack param(*aodTrack);
4801       if (!param.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov)){
4802         continue;
4803       }
4804             
4805       if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
4806         continue;
4807       
4808       Qx2 += TMath::Cos(2*aodTrack->Phi()); 
4809       Qy2 += TMath::Sin(2*aodTrack->Phi());
4810       Qx3 += TMath::Cos(3*aodTrack->Phi()); 
4811       Qy3 += TMath::Sin(3*aodTrack->Phi());
4812       
4813     }
4814     
4815    Float_t evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
4816    Float_t evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
4817
4818     fgPsi2tpc = evPlAng2;
4819     fgPsi3tpc = evPlAng3;
4820
4821      fPhiRPTPC->Fill(iC,evPlAng2);
4822      fPhiRPTPCv3->Fill(iC,evPlAng3);
4823
4824
4825
4826 //V0 info    
4827     AliAODVZERO* aodV0 = event->GetVZEROData();
4828
4829     for (Int_t iv0 = 0; iv0 < 64; iv0++) {
4830       Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
4831       Float_t multv0 = aodV0->GetMultiplicity(iv0);
4832
4833       if(! fIsAfter2011){
4834         if(fAnalysisType == "AOD"){//not for reco MC tracks, only for real data
4835           if (iv0 < 32){ // V0C
4836             Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4837             Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4838             Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4839             Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4840           } else {       // V0A
4841             Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4842             Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4843             Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4844             Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4845           }
4846         }
4847         else{
4848           if (iv0 < 32){ // V0C
4849             Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4850             Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4851             Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4852             Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4853           } else {       // V0A
4854             Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4855             Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4856             Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4857             Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4858           }
4859         }
4860       }
4861     }
4862    //grab for each centrality the proper histo with the Qx and Qy to do the recentering
4863     Double_t Qxamean2 = fMeanQ[iCcal][1][0];
4864     Double_t Qxarms2  = fWidthQ[iCcal][1][0];
4865     Double_t Qyamean2 = fMeanQ[iCcal][1][1];
4866     Double_t Qyarms2  = fWidthQ[iCcal][1][1];
4867     Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
4868     Double_t Qxarms3  = fWidthQv3[iCcal][1][0];
4869     Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
4870     Double_t Qyarms3  = fWidthQv3[iCcal][1][1];
4871     
4872     Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
4873     Double_t Qxcrms2  = fWidthQ[iCcal][0][0];
4874     Double_t Qycmean2 = fMeanQ[iCcal][0][1];
4875     Double_t Qycrms2  = fWidthQ[iCcal][0][1];   
4876     Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
4877     Double_t Qxcrms3  = fWidthQv3[iCcal][0][0];
4878     Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
4879     Double_t Qycrms3  = fWidthQv3[iCcal][0][1]; 
4880     
4881     Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
4882     Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
4883     Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
4884     Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
4885     Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
4886     Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
4887     Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
4888     Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
4889     /*
4890     //to calculate 2nd order event plane with v0M
4891  Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
4892     /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
4893   Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
4894     /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
4895
4896   //here the calculated event plane is within -Pi to +Pi(delete it , no use here , only for definition)
4897   Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
4898   Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
4899   Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
4900
4901     */
4902
4903     Float_t evPlAngV0ACor2=999.;
4904     Float_t evPlAngV0CCor2=999.;
4905     Float_t evPlAngV0ACor3=999.;
4906     Float_t evPlAngV0CCor3=999.;
4907
4908    if(! fIsAfter2011){
4909       if(fAnalysisType == "AOD"){
4910         evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
4911         evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
4912         evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
4913         evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
4914       }
4915       else{
4916         evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
4917         evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
4918         evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
4919         evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
4920       }
4921     }
4922     else{
4923       AliEventplane *ep =  event->GetEventplane();
4924       evPlAngV0ACor2 = ep->GetEventplane("V0A", event, 2);
4925       evPlAngV0CCor2 = ep->GetEventplane("V0C", event, 2);
4926       evPlAngV0ACor3 = ep->GetEventplane("V0A", event, 3);
4927       evPlAngV0CCor3 = ep->GetEventplane("V0C", event, 3);
4928     }
4929
4930     fgPsi2v0a = evPlAngV0ACor2;
4931     fgPsi2v0c = evPlAngV0CCor2;
4932     fgPsi3v0a = evPlAngV0ACor3;
4933     fgPsi3v0c = evPlAngV0CCor3;
4934
4935  // Fill EP distribution histograms evPlAng
4936     
4937      fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
4938      fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
4939     
4940      fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
4941      fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
4942
4943     // Fill histograms needed for resolution evaluation
4944     fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
4945     fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
4946     fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
4947     
4948     fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
4949     fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
4950     fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
4951
4952
4953     /*   
4954  Float_t gVZEROEventPlane    = -10.;
4955   Float_t gReactionPlane      = -10.;
4956   Double_t qxTot = 0.0, qyTot = 0.0;
4957
4958     AliEventplane *ep = event->GetEventplane();
4959     if(ep){ 
4960       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
4961       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
4962       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
4963       gReactionPlane = gVZEROEventPlane;
4964     }
4965     */
4966   }//AOD,ESD,ESDMC
4967  //return gReactionPlane;
4968
4969  //make the final 2nd order event plane within 0 to Pi
4970      //using data and reco tracks only
4971       if(fgPsi2v0a!=999. && fgPsi2v0a < 0.) fgPsi2v0a += TMath::Pi();
4972       if(fgPsi2v0c!=999. && fgPsi2v0c < 0.) fgPsi2v0c += TMath::Pi();
4973       if(fgPsi2tpc!=999. && fgPsi2tpc < 0.) fgPsi2tpc += TMath::Pi();
4974       //using truth tracks only
4975       if(evplaneMC!=999. && evplaneMC < 0.) evplaneMC += TMath::Pi();
4976       if(fgPsi2v0aMC!=999. && fgPsi2v0aMC < 0.) fgPsi2v0aMC += TMath::Pi();
4977       if(fgPsi2v0cMC!=999. && fgPsi2v0cMC < 0.) fgPsi2v0cMC += TMath::Pi();
4978       if(fgPsi2tpcMC!=999. && fgPsi2tpcMC < 0.) fgPsi2tpcMC += TMath::Pi();
4979       //for the time being leave the 3rd order event planes within -pi/3 t0 +pi/3
4980
4981       if(truth){//for truth MC
4982         if(fV2 && fEPdet=="header")gReactionPlane=evplaneMC;
4983         if(fV2 && fEPdet=="V0A")gReactionPlane=fgPsi2v0aMC;
4984         if(fV2 && fEPdet=="V0C")gReactionPlane=fgPsi2v0cMC;
4985         if(fV2 && fEPdet=="TPC")gReactionPlane=fgPsi2tpcMC;
4986
4987         if(fV3 && fEPdet=="V0A")gReactionPlane=fgPsi3v0aMC;
4988         if(fV3 && fEPdet=="V0C")gReactionPlane=fgPsi3v0cMC;
4989         if(fV3 && fEPdet=="TPC")gReactionPlane=fgPsi3tpcMC;
4990 }
4991       else{//for data and recoMC
4992         if(fV2 && fEPdet=="V0A")gReactionPlane=fgPsi2v0a;
4993         if(fV2 && fEPdet=="V0C")gReactionPlane=fgPsi2v0c;
4994         if(fV2 && fEPdet=="TPC")gReactionPlane=fgPsi2tpc;
4995
4996         if(fV3 && fEPdet=="V0A")gReactionPlane=fgPsi3v0a;
4997         if(fV3 && fEPdet=="V0C")gReactionPlane=fgPsi3v0c;
4998         if(fV3 && fEPdet=="TPC")gReactionPlane=fgPsi3tpc;
4999
5000       }
5001      
5002  } //centrality cut condition
5003
5004 return gReactionPlane;
5005 }
5006 //------------------------------------------------------------------------------------------------------------------
5007 void AliTwoParticlePIDCorr::OpenInfoCalbration(Int_t run){
5008     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
5009     TFile *foadb = TFile::Open(oadbfilename.Data());
5010
5011     if(!foadb){
5012         printf("OADB file %s cannot be opened\n",oadbfilename.Data());
5013         return;
5014     }
5015
5016     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
5017     if(!cont){
5018         printf("OADB object hMultV0BefCorr is not available in the file\n");
5019         return; 
5020     }
5021
5022     if(!(cont->GetObject(run))){
5023         printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
5024         run = 137366;
5025     }
5026     fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
5027
5028     TF1 *fpol0 = new TF1("fpol0","pol0"); 
5029     fMultV0->Fit(fpol0,"","",0,31);
5030     fV0Cpol = fpol0->GetParameter(0);
5031     fMultV0->Fit(fpol0,"","",32,64);
5032     fV0Apol = fpol0->GetParameter(0);
5033
5034     for(Int_t iside=0;iside<2;iside++){
5035         for(Int_t icoord=0;icoord<2;icoord++){
5036             for(Int_t i=0;i  < 9;i++){
5037                 char namecont[100];
5038                 if(iside==0 && icoord==0)
5039                   snprintf(namecont,100,"hQxc2_%i",i);
5040                 else if(iside==1 && icoord==0)
5041                   snprintf(namecont,100,"hQxa2_%i",i);
5042                 else if(iside==0 && icoord==1)
5043                   snprintf(namecont,100,"hQyc2_%i",i);
5044                 else if(iside==1 && icoord==1)
5045                   snprintf(namecont,100,"hQya2_%i",i);
5046
5047                 cont = (AliOADBContainer*) foadb->Get(namecont);
5048                 if(!cont){
5049                     printf("OADB object %s is not available in the file\n",namecont);
5050                     return;     
5051                 }
5052                 
5053                 if(!(cont->GetObject(run))){
5054                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5055                     run = 137366;
5056                 }
5057                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5058                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5059
5060                 //for v3
5061                 if(iside==0 && icoord==0)
5062                   snprintf(namecont,100,"hQxc3_%i",i);
5063                 else if(iside==1 && icoord==0)
5064                   snprintf(namecont,100,"hQxa3_%i",i);
5065                 else if(iside==0 && icoord==1)
5066                   snprintf(namecont,100,"hQyc3_%i",i);
5067                 else if(iside==1 && icoord==1)
5068                   snprintf(namecont,100,"hQya3_%i",i);
5069
5070                 cont = (AliOADBContainer*) foadb->Get(namecont);
5071                 if(!cont){
5072                     printf("OADB object %s is not available in the file\n",namecont);
5073                     return;     
5074                 }
5075                 
5076                 if(!(cont->GetObject(run))){
5077                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5078                     run = 137366;
5079                 }
5080                 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5081                 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5082
5083             }
5084         }
5085     }
5086 }
5087 //____________________________________________________________________
5088 void AliTwoParticlePIDCorr::FillPIDEventPlane(Double_t centrality,Int_t par,Float_t trigphi,Float_t fReactionPlane) 
5089 {
5090
5091  // Event plane (determine psi bin)
5092     Double_t gPsiMinusPhi    =   0.;
5093     Double_t gPsiMinusPhiBin = -10.;
5094 if(fRequestEventPlane){
5095     gPsiMinusPhi   = TMath::Abs(trigphi - fReactionPlane);
5096     //in-plane
5097     if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
5098       (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
5099        ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
5100       gPsiMinusPhiBin = 0.0;
5101     //intermediate
5102     else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
5103             ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
5104             ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
5105             ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
5106       gPsiMinusPhiBin = 1.0;
5107     //out of plane
5108     else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
5109             ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
5110       gPsiMinusPhiBin = 2.0;
5111     //everything else
5112     else 
5113       gPsiMinusPhiBin = 3.0;
5114
5115     fEventPlanePID->Fill(centrality,gPsiMinusPhiBin,(Float_t)par); 
5116  }
5117 }
5118
5119 //----------------------------------------------------------
5120 void AliTwoParticlePIDCorr::Terminate(Option_t *) 
5121 {
5122   // Draw result to screen, or perform fitting, normalizations
5123   // Called once at the end of the query
5124   fOutput = dynamic_cast<TList*> (GetOutputData(1));
5125   if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
5126   
5127   
5128 }
5129 //------------------------------------------------------------------