#include "AliAODInputHandler.h"
#include "AliAODMCParticle.h"
#include "AliAODTracklets.h"
+#include "AliAnalysisUtils.h"
#include "AliThreePionRadii.h"
fAODcase(kTRUE),
fPbPbcase(kTRUE),
fGenerateSignal(kFALSE),
+ fGeneratorOnly(kFALSE),
fPdensityPairCut(kTRUE),
fRMax(11),
fFilterBit(7),
fEDbin(0),
fMbins(fCentBins),
fMultLimit(0),
+ fKt3bins(1),
+ fV0Mbinning(kFALSE),
fCentBinLowLimit(0),
fCentBinHighLimit(1),
+ fTriggerType(0),
fEventCounter(0),
fEventsToMix(0),
fZvertexBins(0),
fAODcase(kTRUE),
fPbPbcase(kTRUE),
fGenerateSignal(kFALSE),
+ fGeneratorOnly(kFALSE),
fPdensityPairCut(kTRUE),
fRMax(11),
fFilterBit(7),
fEDbin(0),
fMbins(fCentBins),
fMultLimit(0),
+ fKt3bins(1),
+ fV0Mbinning(kFALSE),
fCentBinLowLimit(0),
fCentBinHighLimit(1),
+ fTriggerType(0),
fEventCounter(0),
fEventsToMix(0),
fZvertexBins(0),
fAODcase(obj.fAODcase),
fPbPbcase(obj.fPbPbcase),
fGenerateSignal(obj.fGenerateSignal),
+ fGeneratorOnly(obj.fGeneratorOnly),
fPdensityPairCut(obj.fPdensityPairCut),
fRMax(obj.fRMax),
fFilterBit(obj.fFilterBit),
fEDbin(obj.fEDbin),
fMbins(obj.fMbins),
fMultLimit(obj.fMultLimit),
+ fKt3bins(obj.fKt3bins),
+ fV0Mbinning(obj.fV0Mbinning),
fCentBinLowLimit(obj.fCentBinLowLimit),
fCentBinHighLimit(obj.fCentBinHighLimit),
+ fTriggerType(obj.fTriggerType),
fEventCounter(obj.fEventCounter),
fEventsToMix(obj.fEventsToMix),
fZvertexBins(obj.fZvertexBins),
fAODcase = obj.fAODcase;
fPbPbcase = obj.fPbPbcase;
fGenerateSignal = obj.fGenerateSignal;
+ fGeneratorOnly = obj.fGeneratorOnly;
fPdensityPairCut = obj.fPdensityPairCut;
fRMax = obj.fRMax;
fFilterBit = obj.fFilterBit;
fEDbin = obj.fEDbin;
fMbins = obj.fMbins;
fMultLimit = obj.fMultLimit;
+ fKt3bins = obj.fKt3bins;
+ fV0Mbinning = obj.fV0Mbinning;
fCentBinLowLimit = obj.fCentBinLowLimit;
fCentBinHighLimit = obj.fCentBinHighLimit;
+ fTriggerType = obj.fTriggerType;
fEventCounter = obj.fEventCounter;
fEventsToMix = obj.fEventsToMix;
fZvertexBins = obj.fZvertexBins;
fMultLimits[20]=2000;
- if(fPbPbcase && fCentBinLowLimit < 10) {// PbPb 0-50%
+ if(fPbPbcase && fCentBinLowLimit < 6) {// PbPb 0-30%, was 0-50%
fMultLimit=kMultLimitPbPb;
fMbins=fCentBins;
fQcut[0]=0.1;//pi-pi, pi-k, pi-p
fQcut[1]=0.1;//k-k
fQcut[2]=0.6;//the rest
- fNormQcutLow[0] = 0.15;//0.15
- fNormQcutHigh[0] = 0.175;//0.175
+ fNormQcutLow[0] = 0.15;// was 0.15
+ fNormQcutHigh[0] = 0.175;// was 0.175
fNormQcutLow[1] = 1.34;//1.34
fNormQcutHigh[1] = 1.4;//1.4
fNormQcutLow[2] = 1.1;//1.1
fQupperBound = fQcut[0];
fQbins = kQbins;
//
- fDampStart = 0.5;// was 0.3
+ fDampStart = 0.5;
fDampStep = 0.02;
- }else if(fPbPbcase && fCentBinLowLimit >= 10) {// PbPb 50-100%
+ }else if(fPbPbcase && fCentBinLowLimit >= 6) {// PbPb 30-100%, was 50-100%
fMultLimit=kMultLimitPbPb;
fMbins=fCentBins;
fQcut[0]=0.2;//pi-pi, pi-k, pi-p
fQcut[1]=0.2;//k-k
fQcut[2]=1.2;//the rest
- fNormQcutLow[0] = 0.3;//0.15
- fNormQcutHigh[0] = 0.35;//0.175
+ fNormQcutLow[0] = 0.3;// was 0.3
+ fNormQcutHigh[0] = 0.35;// was 0.35
fNormQcutLow[1] = 1.34;//1.34
fNormQcutHigh[1] = 1.4;//1.4
fNormQcutLow[2] = 1.1;//1.1
fQupperBound = fQcut[0];
fQbins = 2*kQbins;
//
- fDampStart = 0.5;// was 0.3
+ fDampStart = 0.5;
fDampStep = 0.02;
}else {// pp or pPb
fMultLimit=kMultLimitPP;
fQcut[0]=2.0;// 0.4
fQcut[1]=2.0;
fQcut[2]=2.0;
- fNormQcutLow[0] = 1.0;
- fNormQcutHigh[0] = 1.2;// 1.5
+ fNormQcutLow[0] = 1.0;// was 1.0
+ fNormQcutHigh[0] = 1.2;// was 1.2
fNormQcutLow[1] = 1.0;
fNormQcutHigh[1] = 1.2;
fNormQcutLow[2] = 1.0;
//
fQlimitC2 = 2.0;
fQbinsC2 = 200;
- fQupperBound = 0.4;
+ fQupperBound = 0.5;// was 0.4
fQbins = kQbinsPP;
//
- fDampStart = 0.5;// was 0.3
+ fDampStart = 0.5;
fDampStep = 0.02;
}
- fQLowerCut = 0.005;// was 0.005
+ fQLowerCut = 0.005;
fKupperBound = 1.0;
//
TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
fOutputList->Add(fEvents2);
+ TH1F *fMultDist0 = new TH1F("fMultDist0","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
+ fMultDist0->GetXaxis()->SetTitle("Multiplicity");
+ fOutputList->Add(fMultDist0);
+
TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
fMultDist1->GetXaxis()->SetTitle("Multiplicity");
fOutputList->Add(fMultDist1);
fMultDist3->GetXaxis()->SetTitle("Multiplicity");
fOutputList->Add(fMultDist3);
+ TH1F *fMultDist4 = new TH1F("fMultDist4","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
+ fMultDist4->GetXaxis()->SetTitle("Multiplicity");
+ fOutputList->Add(fMultDist4);
+
TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
fOutputList->Add(fPtEtaDist);
if(fMCcase) fOutputList->Add(fPrimaryMCPionPairs);
TH3D *fAllMCPionPairs = new TH3D("fAllMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
if(fMCcase) fOutputList->Add(fAllMCPionPairs);
-
+ //
+ TH3D *fMuonContamSmearedNum2 = new TH3D("fMuonContamSmearedNum2","",2,-0.5,1.5, 20,0,1, 100,0,0.5);
+ if(fMCcase) fOutputList->Add(fMuonContamSmearedNum2);
+ TH3D *fMuonContamSmearedDen2 = new TH3D("fMuonContamSmearedDen2","",2,-0.5,1.5, 20,0,1, 100,0,0.5);
+ if(fMCcase) fOutputList->Add(fMuonContamSmearedDen2);
+ TH3D *fMuonContamIdealNum2 = new TH3D("fMuonContamIdealNum2","",2,-0.5,1.5, 20,0,1, 100,0,0.5);
+ if(fMCcase) fOutputList->Add(fMuonContamIdealNum2);
+ TH3D *fMuonContamIdealDen2 = new TH3D("fMuonContamIdealDen2","",2,-0.5,1.5, 20,0,1, 100,0,0.5);
+ if(fMCcase) fOutputList->Add(fMuonContamIdealDen2);
+ //
+ TH3D *fMuonContamSmearedNum3 = new TH3D("fMuonContamSmearedNum3","",2,-0.5,1.5, 2,-0.5,1.5, 50,0,0.5);
+ if(fMCcase) fOutputList->Add(fMuonContamSmearedNum3);
+ TH3D *fMuonContamSmearedDen3 = new TH3D("fMuonContamSmearedDen3","",2,-0.5,1.5, 2,-0.5,1.5, 50,0,0.5);
+ if(fMCcase) fOutputList->Add(fMuonContamSmearedDen3);
+ TH3D *fMuonContamIdealNum3 = new TH3D("fMuonContamIdealNum3","",2,-0.5,1.5, 2,-0.5,1.5, 50,0,0.5);
+ if(fMCcase) fOutputList->Add(fMuonContamIdealNum3);
+ TH3D *fMuonContamIdealDen3 = new TH3D("fMuonContamIdealDen3","",2,-0.5,1.5, 2,-0.5,1.5, 50,0,0.5);
+ if(fMCcase) fOutputList->Add(fMuonContamIdealDen3);
+ //
+ TH1D *fMuonParents = new TH1D("fMuonParents","",500,0.5,500.5);
+ if(fMCcase) fOutputList->Add(fMuonParents);
+ TH1D *fSecondaryMuonParents = new TH1D("fSecondaryMuonParents","",500,0.5,500.5);
+ if(fMCcase) fOutputList->Add(fSecondaryMuonParents);
+ TH3D *fMuonPionDeltaQinv = new TH3D("fMuonPionDeltaQinv","",2,-0.5,1.5, 20,0,1, 100,-0.2,0.2);
+ if(fMCcase) fOutputList->Add(fMuonPionDeltaQinv);
+ TH1D *fPionCandidates = new TH1D("fPionCandidates","",500,0.5,500.5);
+ if(fMCcase) fOutputList->Add(fPionCandidates);
+ //
+ TH3D *fMuonPionK2 = new TH3D("fMuonPionK2","",2,-0.5,1.5, 20,0,1, 100,0,0.5);
+ TH3D *fPionPionK2 = new TH3D("fPionPionK2","",2,-0.5,1.5, 20,0,1, 100,0,0.5);
+ TH3D *fMuonPionK3 = new TH3D("fMuonPionK3","",2,-0.5,1.5, 2,-0.5,1.5, 50,0,0.5);
+ TH3D *fPionPionK3 = new TH3D("fPionPionK3","",2,-0.5,1.5, 2,-0.5,1.5, 50,0,0.5);
+ if(fMCcase) fOutputList->Add(fMuonPionK2);
+ if(fMCcase) fOutputList->Add(fPionPionK2);
+ if(fMCcase) fOutputList->Add(fMuonPionK3);
+ if(fMCcase) fOutputList->Add(fPionPionK3);
+ //
TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
fOutputList->Add(fAvgMult);
TH2D *fAvgMultHisto2D = new TH2D("fAvgMultHisto2D","",fMbins,.5,fMbins+.5, 1000,0.5,2000.5);
fOutputList->Add(fAvgMultHisto2D);
-
+ TH2D *fAvgMultHisto2DV0C = new TH2D("fAvgMultHisto2DV0C","",fMbins,.5,fMbins+.5, 1000,0.5,2000.5);
+ fOutputList->Add(fAvgMultHisto2DV0C);
+ TH2D *fAvgMultHisto2DV0AplusC = new TH2D("fAvgMultHisto2DV0AplusC","",fMbins,.5,fMbins+.5, 1000,0.5,2000.5);
+ fOutputList->Add(fAvgMultHisto2DV0AplusC);
TH2D *fTrackChi2NDF = new TH2D("fTrackChi2NDF","",20,0,100, 100,0,10);
fOutputList->Add(fTrackChi2NDF);
fOutputList->Add(fTrackTPCncls);
- TH3D *fTPNRejects1 = new TH3D("fTPNRejects1","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
- fOutputList->Add(fTPNRejects1);
- TH3D *fTPNRejects2 = new TH3D("fTPNRejects2","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
- fOutputList->Add(fTPNRejects2);
- TH3D *fTPNRejects3 = new TH3D("fTPNRejects3","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
- fOutputList->Add(fTPNRejects3);
- TH3D *fTPNRejects4 = new TH3D("fTPNRejects4","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
- fOutputList->Add(fTPNRejects4);
- TH3D *fTPNRejects5 = new TH3D("fTPNRejects5","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
- fOutputList->Add(fTPNRejects5);
-
TH3D *fKt3DistTerm1 = new TH3D("fKt3DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
TH3D *fKt3DistTerm5 = new TH3D("fKt3DistTerm5","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
fOutputList->Add(fMCWeight3DTerm4MCden);
TH2D *fdNchdEtaResponse = new TH2D("fdNchdEtaResponse","",15,0,15, 15,0,15);
- TH2D *fNpionTrueDist = new TH2D("fNpionTrueDist","",fMbins,.5,fMbins+.5, 1000,0.5,2000.5);
- TH2D *fNchTrueDist = new TH2D("fNchTrueDist","",fMbins,.5,fMbins+.5, 1000,0.5,2000.5);
+ TH2D *fNpionTrueDist = new TH2D("fNpionTrueDist","",fMbins,.5,fMbins+.5, 3000,0.5,3000.5);
+ TH2D *fNchTrueDist = new TH2D("fNchTrueDist","",fMbins,.5,fMbins+.5, 3000,0.5,3000.5);// default Nch mapping
+ TH2D *fNchTrueDistCMS = new TH2D("fNchTrueDistCMS","",fMbins,.5,fMbins+.5, 3000,0.5,3000.5);// default Nch mapping
+ TH2D *fNchTrueDistFullPt = new TH2D("fNchTrueDistFullPt","",fMbins,.5,fMbins+.5, 3000,0.5,3000.5);// full Pt Nch range mapping
+ TH2D *fNchTrueDistPubMethod = new TH2D("fNchTrueDistPubMethod","",fMbins,.5,fMbins+.5, 3000,0.5,3000.5);// Published pp Nch mapping
+ Float_t PubBins[9]={1.,12.,17.,23.,29.,35.,42.,52.,152.};
+ TProfile *fAvgNchTrueDistvsPubMethodBin = new TProfile("fAvgNchTrueDistvsPubMethodBin","",8,PubBins,"");
+ TProfile *fAvgRecRate = new TProfile("fAvgRecRate","",3000,0.5,3000.5, 0,3000, "");
if(fMCcase) fOutputList->Add(fdNchdEtaResponse);
if(fMCcase) fOutputList->Add(fNpionTrueDist);
if(fMCcase) fOutputList->Add(fNchTrueDist);
+ if(fMCcase) fOutputList->Add(fNchTrueDistCMS);
+ if(fMCcase) fOutputList->Add(fNchTrueDistFullPt);
+ if(fMCcase) fOutputList->Add(fNchTrueDistPubMethod);
+ if(fMCcase) fOutputList->Add(fAvgRecRate);
+ if(fMCcase) fOutputList->Add(fAvgNchTrueDistvsPubMethodBin);
+
TH2D *fdCentVsNchdEta = new TH2D("fdCentVsNchdEta","",fMbins,.5,fMbins+.5, 15,0,15);
if(fPbPbcase) fOutputList->Add(fdCentVsNchdEta);
+ TH1D *fV0TotSignal = new TH1D("fV0TotSignal","",3000, 0,30000);
+ if(fV0Mbinning) fOutputList->Add(fV0TotSignal);
+ TH2D *fMultBinVsCent = new TH2D("fMultBinVsCent","",fMbins,.5,fMbins+.5, 100,0,100);
+ fOutputList->Add(fMultBinVsCent);
+
+ TH1D *fExtendedQ3Histo_term1 = new TH1D("fExtendedQ3Histo_term1","",50,0,0.5);
+ TH1D *fExtendedQ3Histo_term2 = new TH1D("fExtendedQ3Histo_term2","",50,0,0.5);
+ TH1D *fExtendedQ3Histo_term5 = new TH1D("fExtendedQ3Histo_term5","",50,0,0.5);
+ fOutputList->Add(fExtendedQ3Histo_term1);
+ fOutputList->Add(fExtendedQ3Histo_term2);
+ fOutputList->Add(fExtendedQ3Histo_term5);
+
if(fPdensityPairCut){
for(Int_t mb=0; mb<fMbins; mb++){
if((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit)) continue;
for(Int_t edB=0; edB<fEDbins; edB++){
+ if(edB >= fKt3bins) continue;
+
for(Int_t c1=0; c1<2; c1++){
for(Int_t c2=0; c2<2; c2++){
for(Int_t sc=0; sc<kSCLimit2; sc++){
Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2 = new TH2D(nameEx2->Data(),"Two Particle Distribution",20,0,1, fQbinsC2,0,fQlimitC2);
fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2);
+ //
+ TString *nameMeanKt=new TString(nameEx2->Data());
+ nameMeanKt->Append("_MeanKt");
+ Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMeanKt = new TH1D(nameMeanKt->Data(),"Two Particle Distribution",200,0,1);
+ fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMeanKt);
+ //
TString *nameEx2QW=new TString(nameEx2->Data());
nameEx2QW->Append("_QW");
Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2QW = new TH2D(nameEx2QW->Data(),"Two Particle Distribution",20,0,1, fQbinsC2,0,fQlimitC2);
fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
nameEx2PIDpurityNum->Append("_PIDpurityNum");
- Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum = new TH2D(nameEx2PIDpurityNum->Data(),"Two Particle Distribution",20,0,1, fQbinsC2,0,fQlimitC2);
+ Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum = new TH3D(nameEx2PIDpurityNum->Data(),"Two Particle Distribution",16,0.5,16.5, 20,0,1, fQbinsC2,0,fQlimitC2);
fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
}
Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = new TH3D(name3DQ->Data(),"", fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3);
//
+ TString *nameMeanKt=new TString(namePC3->Data());
+ nameMeanKt->Append("_MeanKt");
+ Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fMeanKt = new TH1D(nameMeanKt->Data(),"", 200,0,1);
+ fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fMeanKt);
if(sc==0 && fMCcase==kTRUE){
TString *name3DMomResIdeal=new TString(namePC3->Data());
-
+ TH1D *frstar4VectDist = new TH1D("frstar4VectDist","",10000,0,100);
+ fOutputList->Add(frstar4VectDist);
TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
TH1D *fQDist = new TH1D("fQDist","",200,-.2,.2);
fOutputList->Add(fQDist);
-
-
+
////////////////////////////////////
///////////////////////////////////
}
//________________________________________________________________________
-void AliThreePionRadii::Exec(Option_t *)
+void AliThreePionRadii::UserExec(Option_t *)
{
// Main loop
// Called for each event
Bool_t isSelected3 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
if(!isSelected1 && !isSelected2 && !isSelected3 && !fMCcase) {return;}
}
+ }else{// pp and pPb
+ Bool_t isSelected[4]={kFALSE};
+ isSelected[0] = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
+ isSelected[1] = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);
+ isSelected[2] = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
+ isSelected[3] = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kHighMult);
+ if(!isSelected[fTriggerType] && !fMCcase) return;
}
-
+
+
///////////////////////////////////////////////////////////
const AliAODVertex *primaryVertexAOD;
AliCentrality *centrality;// for AODs and ESDs
TClonesArray *mcArray = 0x0;
- Int_t mcdNch=0;
- Int_t mcdNpion=0;
+ Int_t mcNch=0, mcNchCMS=0, mcNchFullPt=0, mcNchPubMethod=0;
+ Int_t mcNpion=0;
if(fMCcase){
if(fAODcase){
mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
if(!mcParticle) continue;
- if(fabs(mcParticle->Eta())>0.8) continue;
if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
- if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
- if(!mcParticle->IsPrimary()) continue;
+ //if(!mcParticle->IsPrimary()) continue;// superfluous when IsPhysicalPrimary() is used
if(!mcParticle->IsPhysicalPrimary()) continue;
- mcdNch++;
- if(abs(mcParticle->GetPdgCode())==211) mcdNpion++;
+ //
+ if(fabs(mcParticle->Eta())<=0.5) mcNchPubMethod++;// Published pp binning
+ if(fabs(mcParticle->Eta())<=0.8) mcNchFullPt++;// My binning in full Pt range
+
+ if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
+
+ //
+ if(mcParticle->P() < 1.0) {
+ if(fabs(mcParticle->Eta())<=0.8) {
+ mcNch++;// My binning in my pt range
+ if(abs(mcParticle->GetPdgCode())==211) mcNpion++;
+ }
+ }
+ // p-Pb CMS boost counting
+ Double_t newPz = mcParticle->Pz()*cosh(0.465) - mcParticle->E()*sinh(0.465);
+ Double_t newP = sqrt(pow(mcParticle->Pt(),2) + pow(newPz,2));
+ if(newP < 1.0){
+ Double_t newEta = 0.5 * log( (newP+newPz) / (newP-newPz));
+ if(TMath::Abs(newEta)<=0.8) {
+ mcNchCMS++;
+ }
+ }
}
-
}
}// fMCcase
UInt_t status=0;
Int_t positiveTracks=0, negativeTracks=0;
Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
-
+ Int_t FBTracks=0, AODTracks=0;
+
Double_t vertex[3]={0};
Int_t zbin=0;
Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
/////////////////////////////////////////////////
-
Float_t centralityPercentile=0;
//Float_t cStep=5.0, cStart=0;
//cout<<"AOD multiplicity = "<<fAOD->GetNumberOfTracks()<<endl;
}
-
-
+ ((TH1F*)fOutputList->FindObject("fMultDist0"))->Fill(fAOD->GetNumberOfTracks());
+
+ // Pile-up rejection
+ AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
+ if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
+ else AnaUtil->SetUseMVPlpSelection(kFALSE);
+ Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD);
+ if(pileUpCase) return;
////////////////////////////////
// Vertexing
if(fabs(vertex[2]) > 10) {/*cout<<"Zvertex Out of Range. Skip Event"<<endl;*/ return;} // Z-Vertex Cut
((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
- if(fAOD->IsPileupFromSPD()) {/*cout<<"PileUpEvent. Skip Event"<<endl;*/ return;} // Reject Pile-up events
+ for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
+ AliAODTrack* aodtrack = fAOD->GetTrack(i);
+ if (!aodtrack) continue;
+ AODTracks++;
+ if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
+ FBTracks++;
+ }
+ ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(FBTracks);
+
+ //if(fAOD->IsPileupFromSPD()) {/*cout<<"PileUpEvent. Skip Event"<<endl;*/ return;} // Old Pile-up cut
if(primaryVertexAOD->GetNContributors() < 1) {/*cout<<"Bad Vertex. Skip Event"<<endl;*/ return;}
- ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
+ ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(FBTracks);
fBfield = fAOD->GetMagneticField();
if(tracklets->GetTheta(trackletN) > 1.0904 && tracklets->GetTheta(trackletN) < 2.0512) trackletMult++;// |eta|<0.5 tracklets
}
-
/////////////////////////////
// Create Shuffled index list
Int_t randomIndex[fAOD->GetNumberOfTracks()];
AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
if (!aodtrack) continue;
if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
-
+
status=aodtrack->GetStatus();
-
- if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
+ if(!aodtrack->TestFilterBit(BIT(7))) continue;// AOD filterBit cut
+ if(aodtrack->GetTPCNcls() < 70) continue;// TPC nCluster cut
+
+ // FilterBit Overlap Check
+ if(fFilterBit != 7){
+ Bool_t goodTrackOtherFB = kFALSE;
+ if(fMCcase && fAOD->GetRunNumber()<=126437) goodTrackOtherFB=kTRUE;// FB7 to FB5 mapping in 10f6a MC does not work
+
+ for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
+ AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
+ if(!aodtrack2) continue;
+ if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
+
+ if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
+
+ }
+ if(!goodTrackOtherFB) continue;
+ }
+
if(aodtrack->Pt() < 0.16) continue;
if(fabs(aodtrack->Eta()) > 0.8) continue;
-
+
Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
if(!goodMomentum) continue;
aodtrack->GetXYZ( fTempStruct[myTracks].fX);
if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
- if(fTempStruct[myTracks].fCharge==+1) {
- ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
- ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
- }else {
- ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
- ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
- }
-
- ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
- ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
-
+
// PID section
fTempStruct[myTracks].fElectron = kFALSE;
Float_t signalTPC=0, signalTOF=0;
Double_t integratedTimesTOF[10]={0};
-
- if(fFilterBit != 7 || (fMCcase && !fPbPbcase)) {
- nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron));
- nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon));
- nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion));
- nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon));
- nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton));
+ Bool_t DoPIDWorkAround=kTRUE;
+ //if(fFilterBit == 7) DoPIDWorkAround=kTRUE;
+ if(fMCcase && !fPbPbcase) DoPIDWorkAround=kFALSE;
+ if(DoPIDWorkAround==kFALSE && fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion)) < 900) {
+ nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron);
+ nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon);
+ nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion);
+ nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon);
+ nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton);
//
- nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron));
- nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon));
- nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion));
- nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon));
- nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton));
+ nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron);
+ nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon);
+ nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion);
+ nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon);
+ nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton);
signalTPC = aodtrack->GetTPCsignal();
if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
fTempStruct[myTracks].fTOFhit = kTRUE;
signalTOF = aodtrack->GetTOFsignal();
aodtrack->GetIntegratedTimes(integratedTimesTOF);
}else fTempStruct[myTracks].fTOFhit = kFALSE;
-
+
}else {// FilterBit 7 PID workaround
-
+
for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
if (!aodTrack2) continue;
UInt_t status2=aodTrack2->GetStatus();
- nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron));
- nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon));
- nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion));
- nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon));
- nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton));
+ nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron);
+ nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon);
+ nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion);
+ nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon);
+ nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton);
//
- nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron));
- nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon));
- nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion));
- nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon));
- nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton));
+ nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron);
+ nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon);
+ nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion);
+ nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon);
+ nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton);
signalTPC = aodTrack2->GetTPCsignal();
if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
}// aodTrack2
}// FilterBit 7 PID workaround
-
+ //cout<<nSigmaTPC[2]<<endl;
///////////////////
((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
if(fTempStruct[myTracks].fTOFhit) {
// Use TOF if good hit and above threshold
if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
- if(nSigmaTOF[0]<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
- if(nSigmaTOF[2]<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
- if(nSigmaTOF[3]<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
- if(nSigmaTOF[4]<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
+ if(fabs(nSigmaTOF[0])<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
+ if(fabs(nSigmaTOF[2])<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
+ if(fabs(nSigmaTOF[3])<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
+ if(fabs(nSigmaTOF[4])<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
}else {// TPC info instead
- if(nSigmaTPC[0]<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
- if(nSigmaTPC[2]<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
- if(nSigmaTPC[3]<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
- if(nSigmaTPC[4]<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
+ if(fabs(nSigmaTPC[0])<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
+ if(fabs(nSigmaTPC[2])<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
+ if(fabs(nSigmaTPC[3])<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
+ if(fabs(nSigmaTPC[4])<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
}
if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons
if(!fTempStruct[myTracks].fPion) continue;// only pions
-
+
+
+
+ if(fTempStruct[myTracks].fCharge==+1) {
+ ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
+ ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
+ }else {
+ ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
+ ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
+ }
+
+ ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
+ ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
+
if(fTempStruct[myTracks].fPion) {// pions
myTracks++;
+ if(fMCcase){// muon mothers
+ AliAODMCParticle *tempMCTrack=(AliAODMCParticle*)mcArray->At(abs(aodtrack->GetLabel()));
+ if(abs(tempMCTrack->GetPdgCode())==13 && tempMCTrack->GetMother()>0){// muons
+ AliAODMCParticle *parent=(AliAODMCParticle*)mcArray->At(tempMCTrack->GetMother());
+ if(parent->IsPhysicalPrimary()){
+ ((TH1D*)fOutputList->FindObject("fMuonParents"))->Fill(abs(parent->GetPdgCode()));
+ }else ((TH1D*)fOutputList->FindObject("fSecondaryMuonParents"))->Fill(abs(parent->GetPdgCode()));
+ }
+ ((TH1D*)fOutputList->FindObject("fPionCandidates"))->Fill(abs(tempMCTrack->GetPdgCode()));
+ }
}
}else {// ESD tracks
cout<<"ESDs not supported currently"<<endl;
return;
}
+ // Generator info only
+ if(fMCcase && fGeneratorOnly){
+ myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
+ for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
+ if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
+ if(myTracks >= 1300) continue;// additional cut to limit high mult events which exceed pair # limits
+
+ AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
+ if(!mcParticle) continue;
+ if(fabs(mcParticle->Eta())>0.8) continue;
+ if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
+ if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
+ if(!mcParticle->IsPrimary()) continue;
+ if(!mcParticle->IsPhysicalPrimary()) continue;
+ if(abs(mcParticle->GetPdgCode())!=211) continue;
+
+ fTempStruct[myTracks].fP[0] = mcParticle->Px();
+ fTempStruct[myTracks].fP[1] = mcParticle->Py();
+ fTempStruct[myTracks].fP[2] = mcParticle->Pz();
+ fTempStruct[myTracks].fX[0] = 0.; fTempStruct[myTracks].fX[1] = 0.; fTempStruct[myTracks].fX[2] = 0.;
+
+ fTempStruct[myTracks].fId = myTracks;// use my track counter
+ fTempStruct[myTracks].fLabel = mctrackN;
+ fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
+ if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
+ fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
+ fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
+ fTempStruct[myTracks].fEta = mcParticle->Eta();
+ fTempStruct[myTracks].fCharge = int(mcParticle->Charge()/3.);
+ fTempStruct[myTracks].fDCAXY = 0.;
+ fTempStruct[myTracks].fDCAZ = 0.;
+ fTempStruct[myTracks].fDCA = 0.;
+ fTempStruct[myTracks].fPion = kTRUE;
+ fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
+ fTempStruct[myTracks].fKey = 1;
+
+ myTracks++;
+ pionCount++;
+ }
+ }
+
+
if(myTracks >= 1) {
- ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
+ ((TH1F*)fOutputList->FindObject("fMultDist4"))->Fill(myTracks);
}
fMbin=-1;
for(Int_t i=0; i<fCentBins; i++){
if( pionCount >= fMultLimits[i] && pionCount < fMultLimits[i+1]) {fMbin = fCentBins-i-1; break;}
-
- //if(fPbPbcase){
- //if(centralityPercentile >= 5*i && centralityPercentile < 5*(i+1)) {fMbin=i; break;}
- //}else{
- //if( (pow(trackletMult,1/3.) >= fMultLimits[i]) && (pow(trackletMult,1/3.) < fMultLimits[i+1])) {fMbin=i; break;}
- //}
- //if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}// Mbin 0 has 0-5 pions
}
-
-
-
- if(fMbin==-1) {cout<<pionCount<<" Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
- if(fMbin < fCentBinLowLimit || fMbin > fCentBinHighLimit) {cout<<"Mult out of range"<<endl; return;}
-
+
fFSIindex=0;
if(fPbPbcase){
if(fMbin==0) fFSIindex = 0;//0-5%
else fFSIindex = 8;//90-100%
}else fFSIindex = 9;// pp and pPb
- //////////////////////////////////////////////////
- fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
- //////////////////////////////////////////////////
-
- //cout<<fMbin<<" "<<pionCount<<endl;
+ if(fMCcase){// FSI binning for MC
+ if(fRMax>=10) fFSIindex = 0;
+ else if(fRMax>=9) fFSIindex = 1;
+ else if(fRMax>=8) fFSIindex = 2;
+ else if(fRMax>=7) fFSIindex = 3;
+ else if(fRMax>=6) fFSIindex = 4;
+ else if(fRMax>=5) fFSIindex = 5;
+ else if(fRMax>=4) fFSIindex = 6;
+ else if(fRMax>=3) fFSIindex = 7;
+ else if(fRMax>=2) fFSIindex = 8;
+ else fFSIindex = 9;
+ }
+
+ if(fV0Mbinning){
+ Bool_t useV0=kFALSE;
+ if(fPbPbcase) useV0=kTRUE;
+ if(!fPbPbcase && fAOD->GetRunNumber() >= 195344 && fAOD->GetRunNumber() <= 195677) useV0=kTRUE;
+ if(useV0){
+ AliAODVZERO *vZero = fAOD->GetVZEROData();
+ Float_t vZeroAmp = vZero->GetMTotV0A();
+ centrality = fAOD->GetCentrality();
+ centralityPercentile = centrality->GetCentralityPercentile("V0M");
+ for(Int_t i=0; i<fCentBins; i++){
+ if(vZeroAmp/4.4 >= fMultLimits[i] && vZeroAmp/4.4 < fMultLimits[i+1]) {fMbin = fCentBins-i-1; break;}
+ }
+ ((TH1D*)fOutputList->FindObject("fV0TotSignal"))->Fill(vZeroAmp);
+ //cout<<centralityPercentile<<" "<<vZeroAmp<<" "<<fMbin<<endl;
+ //
+ Int_t fMbinV0C=-1;
+ vZeroAmp = vZero->GetMTotV0C();
+ for(Int_t i=0; i<fCentBins; i++){
+ if(vZeroAmp/4.4 >= fMultLimits[i] && vZeroAmp/4.4 < fMultLimits[i+1]) {fMbinV0C = fCentBins-i-1; break;}
+ }
+ //
+ Int_t fMbinV0AplusC=-1;
+ vZeroAmp = vZero->GetMTotV0A() + vZero->GetMTotV0C();
+ for(Int_t i=0; i<fCentBins; i++){
+ if(vZeroAmp/4.4 >= fMultLimits[i] && vZeroAmp/4.4 < fMultLimits[i+1]) {fMbinV0AplusC = fCentBins-i-1; break;}
+ }
+ ((TH2D*)fOutputList->FindObject("fAvgMultHisto2DV0C"))->Fill(fMbinV0C+1., pionCount);
+ ((TH2D*)fOutputList->FindObject("fAvgMultHisto2DV0AplusC"))->Fill(fMbinV0AplusC+1., pionCount);
+ }
+ }
+ if(fMbin==-1) {cout<<pionCount<<" Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
+
((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
((TH2D*)fOutputList->FindObject("fAvgMultHisto2D"))->Fill(fMbin+1., pionCount);
if(fMCcase){
- ((TH2D*)fOutputList->FindObject("fdNchdEtaResponse"))->Fill(pow(trackletMult,1/3.), pow(mcdNch,1/3.));
- ((TH2D*)fOutputList->FindObject("fNpionTrueDist"))->Fill(fMbin+1., mcdNpion);
- ((TH2D*)fOutputList->FindObject("fNchTrueDist"))->Fill(fMbin+1., mcdNch);
+ ((TH2D*)fOutputList->FindObject("fdNchdEtaResponse"))->Fill(pow(trackletMult,1/3.), pow(mcNch,1/3.));
+ ((TH2D*)fOutputList->FindObject("fNpionTrueDist"))->Fill(fMbin+1., mcNpion);
+ ((TH2D*)fOutputList->FindObject("fNchTrueDist"))->Fill(fMbin+1., mcNch);// default Nch mapping
+ ((TH2D*)fOutputList->FindObject("fNchTrueDistCMS"))->Fill(fMbin+1., mcNchCMS);// p-Pb CMS counting
+ ((TH2D*)fOutputList->FindObject("fNchTrueDistFullPt"))->Fill(fMbin+1., mcNchFullPt);// full Pt Nch range mapping
+ ((TH2D*)fOutputList->FindObject("fNchTrueDistPubMethod"))->Fill(fMbin+1., mcNchPubMethod);// Published pp Method Nch mapping
+ ((TProfile*)fOutputList->FindObject("fAvgNchTrueDistvsPubMethodBin"))->Fill(mcNchPubMethod, mcNch);// Mapping of Published bins to default Nch bins
+ ((TProfile*)fOutputList->FindObject("fAvgRecRate"))->Fill(mcNpion, pionCount);
}
if(fPbPbcase){
((TH2D*)fOutputList->FindObject("fdCentVsNchdEta"))->Fill(fMbin+1, pow(trackletMult,1/3.));
+ centrality = fAOD->GetCentrality();
+ centralityPercentile = centrality->GetCentralityPercentile("V0M");
+ ((TH2D*)fOutputList->FindObject("fMultBinVsCent"))->Fill(fMbin+1, centralityPercentile);
}
+
+ // Mult cut
+ if(fMbin < fCentBinLowLimit || fMbin > fCentBinHighLimit) {cout<<"Mult out of range"<<endl; return;}
- //cout<<trackletMult<<" "<<mcdNchdEta<<endl;
+ //////////////////////////////////////////////////
+ fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
+ //////////////////////////////////////////////////
+
+
+
+ //return;// un-comment for a run to calculate Nrec to Nch Mapping
+ // to test the eta dependence of radii
+ /*Int_t firstTrackCount=myTracks;
+ Int_t newTrackCount=0;
+ myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
+ for(Int_t newTracks=0; newTracks<firstTrackCount; newTracks++){
+
+ if(fTempStruct[newTracks].fEta > -0.4) continue;
+
+ fTempStruct[newTrackCount]=fTempStruct[newTracks];
+
+ newTrackCount++;
+ pionCount++;
+ }
+ myTracks=newTrackCount;// re-assign main counter
+ */
////////////////////////////////////
// Add event to buffer if > 0 tracks
(fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
(fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
(fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
- }
+ }
}
}
-
-
+
Float_t qinv12=0, qinv13=0, qinv23=0;
Float_t qout=0, qside=0, qlong=0;
}
+
+ /////////////////////////////////////////////////////////////////
+ // extended range Q3 baseline
+ /*for(Int_t iter=0; iter<3; iter++){
+ for (Int_t i=0; i<myTracks; i++) {
+
+ Int_t en2=0;
+ if(iter==2) en2=1;
+ Int_t start2=i+1;
+ if(en2!=0) start2=0;
+ // 2nd particle
+ for (Int_t j=start2; j<(fEvt+en2)->fNtracks; j++) {
+ if((fEvt)->fTracks[i].fCharge != (fEvt+en2)->fTracks[j].fCharge) continue;
+ key1 = (fEvt)->fTracks[i].fKey;
+ key2 = (fEvt+en2)->fTracks[j].fKey;
+ Short_t fillIndex2 = FillIndex2part(key1+key2);
+ pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
+ pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
+ pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
+ pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
+ qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
+
+ if(qinv12>0.5) continue;
+ if(!AcceptPair(&((fEvt)->fTracks[i]), &((fEvt+en2)->fTracks[j]))) continue;
+
+ Int_t en3=0;
+ if(iter==1) en3=1;
+ if(iter==2) en3=2;
+ Int_t start3=j+1;
+ if(iter>0) start3=0;
+ // 3nd particle
+ for (Int_t k=start3; k<(fEvt+en3)->fNtracks; k++) {
+ if((fEvt)->fTracks[i].fCharge != (fEvt+en3)->fTracks[k].fCharge) continue;
+ pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
+ pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
+ pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
+ pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
+ qinv13 = GetQinv(fillIndex2, pVect1, pVect3);
+ if(qinv13>0.5) continue;
+ qinv23 = GetQinv(fillIndex2, pVect2, pVect3);
+ if(qinv23>0.5) continue;
+ if(!AcceptPair(&((fEvt)->fTracks[i]), &((fEvt+en3)->fTracks[k]))) continue;
+ if(!AcceptPair(&((fEvt+en2)->fTracks[j]), &((fEvt+en3)->fTracks[k]))) continue;
+
+ q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
+
+ if(iter==0) ((TH1D*)fOutputList->FindObject("fExtendedQ3Histo_term1"))->Fill(q3);
+ if(iter==1) ((TH1D*)fOutputList->FindObject("fExtendedQ3Histo_term2"))->Fill(q3);
+ if(iter==2) ((TH1D*)fOutputList->FindObject("fExtendedQ3Histo_term5"))->Fill(q3);
+
+ }
+ }
+ }
+ }
+ */
///////////////////////////////////////////////////////
// Start the pairing process
// P11 pairing
// 1st Particle
-
+
for (Int_t i=0; i<myTracks; i++) {
Int_t en2=0;
// Pair Splitting/Merging cut
if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
- if(ch1 == ch2){
+ if(ch1 == ch2 && !fGeneratorOnly){
if(!AcceptPair(&((fEvt)->fTracks[i]), &((fEvt+en2)->fTracks[j]))) {
fPairSplitCut[0][i]->AddAt('1',j);
((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMCqinv->Fill(qinv12MC, MCWeight(ch1,ch2,rForQW,10,qinv12MC));// was 4,5
Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(ch1,ch2,rForQW,10,qinv12MC));// was 4,5
- // pion purity
+ // pion purity
Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityDen->Fill(transK12, qinv12);
- if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){// Pions
- Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityNum->Fill(transK12, qinv12);
- }
-
+ Int_t SCNumber = 1;
+
+ if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==22) SCNumber=1;// e-e
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==24) SCNumber=2;// e-mu
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==222) SCNumber=3;// e-pi
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==332) SCNumber=4;// e-k
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==2223) SCNumber=5;// e-p
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==26) SCNumber=6;// mu-mu
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==224) SCNumber=7;// mu-pi
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==334) SCNumber=8;// mu-k
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==2225) SCNumber=9;// mu-p
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==422) SCNumber=10;// pi-pi
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==532) SCNumber=11;// pi-k
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==2423) SCNumber=12;// pi-p
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==642) SCNumber=13;// k-k
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==2533) SCNumber=14;// k-p
+ else if((abs(mcParticle1->GetPdgCode()) + abs(mcParticle2->GetPdgCode()))==4424) SCNumber=15;// p-p
+ else SCNumber=16;
+
+ Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityNum->Fill(SCNumber, transK12, qinv12);
+
+
+ ///////////////////////
+ // muon contamination
+ if(qinv12 < fQcut[0] && ((fEvt)->fTracks[i].fLabel != (fEvt+en2)->fTracks[j].fLabel)){
+ if(abs(mcParticle1->GetPdgCode())==13 || abs(mcParticle2->GetPdgCode())==13){// muon check
+ Float_t Pparent1[4]={pVect1MC[0],pVect1MC[1],pVect1MC[2],pVect1MC[3]};
+ Float_t Pparent2[4]={pVect2MC[0],pVect2MC[1],pVect2MC[2],pVect2MC[3]};
+ Bool_t pionParent1=kFALSE, pionParent2=kFALSE;
+ if(abs(mcParticle1->GetPdgCode())==13) {
+ AliAODMCParticle *parent1=(AliAODMCParticle*)mcArray->At(mcParticle1->GetMother());
+ if(abs(parent1->GetPdgCode())==211) {
+ pionParent1=kTRUE;
+ Pparent1[1] = parent1->Px(); Pparent1[2] = parent1->Py(); Pparent1[3] = parent1->Pz();
+ Pparent1[0] = sqrt(pow(Pparent1[1],2)+pow(Pparent1[2],2)+pow(Pparent1[3],2)+pow(fTrueMassPi,2));
+ }
+ }
+ //
+ if(abs(mcParticle2->GetPdgCode())==13) {
+ AliAODMCParticle *parent2=(AliAODMCParticle*)mcArray->At(mcParticle2->GetMother());
+ if(abs(parent2->GetPdgCode())==211) {
+ pionParent2=kTRUE;
+ Pparent2[1] = parent2->Px(); Pparent2[2] = parent2->Py(); Pparent2[3] = parent2->Pz();
+ Pparent2[0] = sqrt(pow(Pparent2[1],2)+pow(Pparent2[2],2)+pow(Pparent2[3],2)+pow(fTrueMassPi,2));
+ }
+ }
+
+ Float_t parentQinv12 = GetQinv(0, Pparent1, Pparent2);
+ Float_t WInput = 1.0, muonPionK12=1.0, pionPionK12=1.0;
+ if(parentQinv12 > 0.001 && parentQinv12 < 0.3) {
+ WInput = MCWeight(ch1,ch2, fRMax, 25, parentQinv12);
+ muonPionK12 = FSICorrelation2(ch1, ch2, qinv12MC);
+ pionPionK12 = FSICorrelation2(ch1, ch2, parentQinv12);
+ }
+ Int_t ChComb=0;
+ if(ch1 != ch2) ChComb=1;
+ if(pionParent1 || pionParent2){
+ ((TH3D*)fOutputList->FindObject("fMuonContamSmearedNum2"))->Fill(ChComb, transK12, qinv12MC, WInput);
+ ((TH3D*)fOutputList->FindObject("fMuonContamSmearedDen2"))->Fill(ChComb, transK12, qinv12MC);
+ ((TH3D*)fOutputList->FindObject("fMuonContamIdealNum2"))->Fill(ChComb, transK12, parentQinv12, WInput);
+ ((TH3D*)fOutputList->FindObject("fMuonContamIdealDen2"))->Fill(ChComb, transK12, parentQinv12);
+ ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(ChComb, transK12, qinv12MC-parentQinv12);
+ ((TH3D*)fOutputList->FindObject("fMuonPionK2"))->Fill(ChComb, transK12, qinv12MC, muonPionK12);
+ ((TH3D*)fOutputList->FindObject("fPionPionK2"))->Fill(ChComb, transK12, parentQinv12, pionPionK12);
+ }
+ ////////////////////////////////////
+ // 3rd particle
+ Int_t en3=0;
+ for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {
+ pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
+ pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
+ pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
+ pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
+ //
+ qinv13 = GetQinv(0, pVect1, pVect3);
+ qinv23 = GetQinv(0, pVect2, pVect3);
+ if(qinv13 > fQcut[0] || qinv23 > fQcut[0]) continue;
+
+ if(qinv13 < fQLowerCut || qinv23 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
+ if(ch1 == ch3 && !fGeneratorOnly){
+ if(!AcceptPair(&((fEvt)->fTracks[i]), &((fEvt+en3)->fTracks[k]))) {
+ continue;
+ }
+ }
+ if(ch2 == ch3 && !fGeneratorOnly){
+ if(!AcceptPair(&((fEvt+en2)->fTracks[j]), &((fEvt+en3)->fTracks[k]))) {
+ continue;
+ }
+ }
+
+ if((fEvt+en3)->fTracks[k].fLabel == (fEvt+en2)->fTracks[j].fLabel) continue;
+ if((fEvt+en3)->fTracks[k].fLabel == (fEvt)->fTracks[i].fLabel) continue;
+
+ if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
+ AliAODMCParticle *mcParticle3 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en3)->fTracks[k].fLabel));
+
+ ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
+ pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
+ pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
+ pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
+ pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
+ qinv13MC = GetQinv(0, pVect1MC, pVect3MC);
+ qinv23MC = GetQinv(0, pVect2MC, pVect3MC);
+
+ q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
+ transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
+ Int_t K3index=0;
+ if(transK3>0.3) K3index=1;
+
+ Float_t Pparent3[4]={pVect3MC[0],pVect3MC[1],pVect3MC[2],pVect3MC[3]};
+ Bool_t pionParent3=kFALSE;
+ if(abs(mcParticle3->GetPdgCode())==13){// muon check
+ AliAODMCParticle *parent3=(AliAODMCParticle*)mcArray->At(mcParticle3->GetMother());
+ if(abs(parent3->GetPdgCode())==211) {
+ pionParent3=kTRUE;
+ Pparent3[1] = parent3->Px(); Pparent3[2] = parent3->Py(); Pparent3[3] = parent3->Pz();
+ Pparent3[0] = sqrt(pow(Pparent3[1],2)+pow(Pparent3[2],2)+pow(Pparent3[3],2)+pow(fTrueMassPi,2));
+ }
+ }
+
+ Float_t parentQinv13 = GetQinv(0, Pparent1, Pparent3);
+ Float_t parentQinv23 = GetQinv(0, Pparent2, Pparent3);
+ Float_t parentQ3 = sqrt(pow(parentQinv12,2) + pow(parentQinv13,2) + pow(parentQinv23,2));
+ if(parentQ3 >= 0.5) continue;
+ if(parentQinv12 < 0.001) continue;
+ if(parentQinv13 < 0.001) continue;
+ if(parentQinv23 < 0.001) continue;
+
+ if(!pionParent1 && !pionParent2 && !pionParent3) continue;// want at least one pion-->muon
+
+ Int_t ChCombtriplet=0;
+ if(ch1!=ch2 || ch1!=ch3 || ch2!=ch3) ChCombtriplet=1;
+ Float_t WInput3=1.0;
+ Float_t muonPionK13=1.0, muonPionK23=1.0;
+ Float_t pionPionK13=1.0, pionPionK23=1.0;
+ muonPionK13 = FSICorrelation2(ch1, ch3, qinv13MC);
+ pionPionK13 = FSICorrelation2(ch1, ch3, parentQinv13);
+ muonPionK23 = FSICorrelation2(ch2, ch3, qinv23MC);
+ pionPionK23 = FSICorrelation2(ch2, ch3, parentQinv23);
+ if(ChCombtriplet==0) WInput3 = MCWeight3D(kTRUE, 1, 25, parentQinv12, parentQinv13, parentQinv23);
+ else{
+ if(ch1==ch2) WInput3 = MCWeight3D(kFALSE, 1, 25, parentQinv12, parentQinv13, parentQinv23);
+ else if(ch1==ch3) WInput3 = MCWeight3D(kFALSE, 1, 25, parentQinv13, parentQinv12, parentQinv23);
+ else WInput3 = MCWeight3D(kFALSE, 1, 25, parentQinv23, parentQinv12, parentQinv13);
+ }
+ if(WInput3>0 && WInput3<10.) {
+ ((TH3D*)fOutputList->FindObject("fMuonContamSmearedNum3"))->Fill(ChCombtriplet, K3index, q3MC, WInput3);
+ ((TH3D*)fOutputList->FindObject("fMuonContamSmearedDen3"))->Fill(ChCombtriplet, K3index, q3MC);
+ ((TH3D*)fOutputList->FindObject("fMuonContamIdealNum3"))->Fill(ChCombtriplet, K3index, parentQ3, WInput3);
+ ((TH3D*)fOutputList->FindObject("fMuonContamIdealDen3"))->Fill(ChCombtriplet, K3index, parentQ3);
+ ((TH3D*)fOutputList->FindObject("fMuonPionK3"))->Fill(ChCombtriplet, K3index, q3MC, muonPionK12*muonPionK13*muonPionK23);
+ ((TH3D*)fOutputList->FindObject("fPionPionK3"))->Fill(ChCombtriplet, K3index, parentQ3, pionPionK12*pionPionK13*pionPionK23);
+ }
+ }//label check of 3
+ }// 3rd particle
+ }// muon code check of 1 and 2
+ }// qinv12 cut
}// label check 2
}// MC case
Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fAvgP->Fill(transK12, qinv12, (fEvt)->fTracks[i].fMom);
Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fAvgP->Fill(transK12, qinv12, (fEvt+en2)->fTracks[j].fMom);
-
+ //
+ if(qinv12<fQupperBound) Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMeanKt->Fill(transK12);
//////////////////////////////////////////
}// j particle
}// i particle
-
+
//////////////////////////////////////////////
// P12 pairing
//if(transK12 <= 0.35) fEDbin=0;
//else fEDbin=1;
-
-
+
///////////////////////////////
ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fSmeared->Fill(denIndex, qinv12);
}
}
-
-
+
/////////////////////////////////////////////////////
if(qinv12 <= fQcut[qCutBin]) {// 3-particle MRC
q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
+ if(qinv12 < fQLowerCut) continue;
+ if(qinv13 < fQLowerCut) continue;
+ if(qinv23 < fQLowerCut) continue;
+ if(ch1 == ch2){
+ if(!AcceptPair(&((fEvt)->fTracks[i]), &((fEvt+en2)->fTracks[j]))) continue;
+ }
+ if(ch1 == ch3){
+ if(!AcceptPair(&((fEvt)->fTracks[i]), &((fEvt+en3)->fTracks[k]))) continue;
+ }
+ if(ch2 == ch3){
+ if(!AcceptPair(&((fEvt+en2)->fTracks[j]), &((fEvt+en3)->fTracks[k]))) continue;
+ }
//
// The below call to SetFillBins3 will work for all 3-particle terms since all are for fully mixed events. part is set to 1, but only matters for terms 2-4.
if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
- if(ch1 == ch2){
+ if(ch1 == ch2 && !fGeneratorOnly){
if(!AcceptPair(&((fEvt)->fTracks[i]), &((fEvt+en2)->fTracks[j]))) {
fPairSplitCut[1][i]->AddAt('1',j);
continue;
// 2-particle term
Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
-
+ if(qinv12<fQupperBound) Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMeanKt->Fill(transK12);
//////////////////////////////////////////
ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
- if(ch1 == ch2){
+ if(ch1 == ch2 && !fGeneratorOnly){
if(!AcceptPair(&((fEvt)->fTracks[i]), &((fEvt+en2)->fTracks[j]))) {
fPairSplitCut[2][i]->AddAt('1',j);
continue;
ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
- if(ch1 == ch2){
+ if(ch1 == ch2 && !fGeneratorOnly){
if(!AcceptPair(&((fEvt+en1)->fTracks[i]), &((fEvt+en2)->fTracks[j]))) {
fPairSplitCut[3][i]->AddAt('1',j);
continue;
// if all three pair cuts are the same then the case (config=2 && term=2) never reaches here.
-
+ Float_t Kt12 = sqrt( pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
+ Float_t Kt13 = sqrt( pow(pVect1[1]+pVect3[1],2) + pow(pVect1[2]+pVect3[2],2))/2.;
+ Float_t Kt23 = sqrt( pow(pVect2[1]+pVect3[1],2) + pow(pVect2[2]+pVect3[2],2))/2.;
q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
- if(fEDbins>1){
+ if(fKt3bins==1) fEDbin=0;
+ else if(fKt3bins==2){
if(transK3<0.3) fEDbin=0;
else fEDbin=1;
+ }else{// fKt3bins==3 is the limit set by fEDbins
+ if(transK3<0.25) fEDbin=0;
+ else if(transK3<0.35) fEDbin=1;
+ else fEDbin=2;
}
+
firstQ=0; secondQ=0; thirdQ=0;
////
Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTermsQ3->Fill(q3, WInput);
Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTerms3->Fill(firstQ, secondQ, thirdQ, WInput);
+ //
+ Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fMeanKt->Fill(Kt12);
+ Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fMeanKt->Fill(Kt13);
+ Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fMeanKt->Fill(Kt23);
////
//
if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
////
Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fTermsQ3->Fill(q3, WInput);
Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fTerms3->Fill(firstQ, secondQ, thirdQ, WInput);
-
+ //
+ Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fMeanKt->Fill(Kt12);
+ Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fMeanKt->Fill(Kt13);
+ Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fMeanKt->Fill(Kt23);
}
}
ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, 5, firstQ, secondQ, thirdQ);
Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTermsQ3->Fill(q3);
Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTerms3->Fill(firstQ, secondQ, thirdQ);
+ //
+ Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fMeanKt->Fill(Kt12);
+ Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fMeanKt->Fill(Kt13);
+ Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fMeanKt->Fill(Kt23);
}