877f679dbfc06fba4c39c2620ace3c42e7b5d61a
[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   }