#include "AliAODInputHandler.h"
#include "AliAODMCParticle.h"
#include "AliAODTracklets.h"
+#include "AliAnalysisUtils.h"
#include "AliThreePionRadii.h"
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
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
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;// was 0.4
+ fQupperBound = 0.5;// was 0.4
fQbins = kQbinsPP;
//
fDampStart = 0.5;
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);
TH2D *fdNchdEtaResponse = new TH2D("fdNchdEtaResponse","",15,0,15, 15,0,15);
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.};
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);
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);
}
-
+ 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
}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::kAnyINT);
+ 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;
TClonesArray *mcArray = 0x0;
- Int_t mcNch=0, mcNchFullPt=0, mcNchPubMethod=0;
+ Int_t mcNch=0, mcNchCMS=0, mcNchFullPt=0, mcNchPubMethod=0;
Int_t mcNpion=0;
if(fMCcase){
if(fAODcase){
AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
if(!mcParticle) continue;
if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
- if(!mcParticle->IsPrimary()) continue;
+ //if(!mcParticle->IsPrimary()) continue;// superfluous when IsPhysicalPrimary() is used
if(!mcParticle->IsPhysicalPrimary()) continue;
//
- if(fabs(mcParticle->Eta())>0.8) continue;
if(fabs(mcParticle->Eta())<=0.5) mcNchPubMethod++;// Published pp binning
- mcNchFullPt++;// My binning in full Pt range
+ if(fabs(mcParticle->Eta())<=0.8) mcNchFullPt++;// My binning in full Pt range
+
if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
- mcNch++;// My binning in my pt range
- if(abs(mcParticle->GetPdgCode())==211) mcNpion++;
+
+ //
+ 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
}
- //cout<<fAOD->GetFiredTriggerClasses()<<endl;
/////////////////////////////
// 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(7))) continue;// AOD filterBit cut
+ if(aodtrack->GetTPCNcls() < 70) continue;// TPC nCluster cut
// FilterBit Overlap Check
if(fFilterBit != 7){
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(fFilterBit == 7) DoPIDWorkAround=kTRUE;
if(fMCcase && !fPbPbcase) DoPIDWorkAround=kFALSE;
if(DoPIDWorkAround==kFALSE && fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion)) < 900) {
- 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));
+ 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++) {
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
// 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) {
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;
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);
}
else fFSIindex = 8;//90-100%
}else fFSIindex = 9;// pp and pPb
-
+ 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;
((TH2D*)fOutputList->FindObject("fAvgMultHisto2DV0AplusC"))->Fill(fMbinV0AplusC+1., pionCount);
}
}
-
+
if(fMbin==-1) {cout<<pionCount<<" Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
((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
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;
// Start the pairing process
// P11 pairing
// 1st Particle
-
+
for (Int_t i=0; i<myTracks; i++) {
Int_t en2=0;
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
}// 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.
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff