X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG3%2FvertexingHF%2FAliAnalysisTaskSEDplus.cxx;h=01b0ef6a45223f5303741a4d1499ce1b2809d7fa;hb=9877ffb324dfe769ef97e7bbefe80256989f5735;hp=7fabe861a6fd0b291a072a32c38dc72d2c12a09c;hpb=a96083b9325418d42cc7d5ffe2f653ef38f185e8;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx b/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx index 7fabe861a6f..01b0ef6a452 100644 --- a/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx +++ b/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx @@ -13,6 +13,8 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ + //************************************************************************* // Class AliAnalysisTaskSEDplus // AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and @@ -28,7 +30,6 @@ #include #include #include -#include #include #include "AliAnalysisManager.h" @@ -37,8 +38,6 @@ #include "AliAODEvent.h" #include "AliAODVertex.h" #include "AliAODTrack.h" -#include "AliAODMCHeader.h" -#include "AliAODMCParticle.h" #include "AliAODRecoDecayHF3Prong.h" #include "AliAnalysisVertexingHF.h" #include "AliAnalysisTaskSE.h" @@ -71,51 +70,101 @@ AliAnalysisTaskSE(), fFillNtuple(kFALSE), fReadMC(kFALSE), fUseStrangeness(kFALSE), - fDoLS(kFALSE) + fUseBit(kTRUE), + fCutsDistr(kFALSE), + fDoImpPar(kFALSE), + fNImpParBins(400), + fLowerImpPar(-2000.), + fHigherImpPar(2000.), + fDoLS(0) { // Default constructor } //________________________________________________________________________ AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,AliRDHFCutsDplustoKpipi *dpluscutsprod,Bool_t fillNtuple): -AliAnalysisTaskSE(name), -fOutput(0), -fHistNEvents(0), -fPtVsMass(0), -fPtVsMassTC(0), -fYVsPt(0), -fYVsPtTC(0), -fYVsPtSig(0), -fYVsPtSigTC(0), -fNtupleDplus(0), -fUpmasslimit(1.965), -fLowmasslimit(1.765), -fNPtBins(0), -fBinWidth(0.002), -fListCuts(0), -fRDCutsProduction(dpluscutsprod), -fRDCutsAnalysis(dpluscutsana), -fCounter(0), -fFillNtuple(fillNtuple), -fReadMC(kFALSE), -fUseStrangeness(kFALSE), -fDoLS(kFALSE) + AliAnalysisTaskSE(name), + fOutput(0), + fHistNEvents(0), + fPtVsMass(0), + fPtVsMassTC(0), + fYVsPt(0), + fYVsPtTC(0), + fYVsPtSig(0), + fYVsPtSigTC(0), + fNtupleDplus(0), + fUpmasslimit(1.965), + fLowmasslimit(1.765), + fNPtBins(0), + fBinWidth(0.002), + fListCuts(0), + fRDCutsProduction(dpluscutsprod), + fRDCutsAnalysis(dpluscutsana), + fCounter(0), + fFillNtuple(fillNtuple), + fReadMC(kFALSE), + fUseStrangeness(kFALSE), + fUseBit(kTRUE), + fCutsDistr(kFALSE), + fDoImpPar(kFALSE), + fNImpParBins(400), + fLowerImpPar(-2000.), + fHigherImpPar(2000.), + fDoLS(0) { // // Standrd constructor // - //Double_t ptlim[5]={0.,2.,3.,5,9999999.}; - //SetPtBinLimit(5, ptlim); - SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits()); + fNPtBins=fRDCutsAnalysis->GetNPtBins(); + + for(Int_t i=0;i<3;i++){ + if(fHistCentrality[i])fHistCentrality[i]=0; + if(fCorreld0Kd0pi[i])fCorreld0Kd0pi[i]=0; + } + + for(Int_t i=0; i<5; i++)fHistMassPtImpParTC[i]=0; + + + for(Int_t i=0;i<3*fNPtBins;i++){ + if(fMassHist[i])fMassHist[i]=0; + if(fCosPHist[i])fCosPHist[i]=0; + if(fDLenHist[i])fDLenHist[i]=0; + if(fSumd02Hist[i])fSumd02Hist[i]=0; + if(fSigVertHist[i])fSigVertHist[i]=0; + if(fPtMaxHist[i])fPtMaxHist[i]=0; + if(fPtKHist[i])fPtKHist[i]=0; + if(fPtpi1Hist[i])fPtpi1Hist[i]=0; + if(fPtpi2Hist[i])fPtpi2Hist[i]=0; + if(fDCAHist[i])fDCAHist[i]=0; + if(fMassHistTC[i])fMassHistTC[i]=0; + if(fMassHistTCPlus[i])fMassHistTCPlus[i]=0; + if(fMassHistTCMinus[i])fMassHistTCMinus[i]=0; + + if(fDLxy[i])fDLxy[i]=0; + if(fDLxyTC[i])fDLxyTC[i]=0; + if(fCosxy[i])fCosxy[i]=0; + if(fCosxyTC[i])fCosxyTC[i]=0; + if(fMassHistLS[i])fMassHistLS[i]=0; + if(fCosPHistLS[i])fCosPHistLS[i]=0; + if(fDLenHistLS[i])fDLenHistLS[i]=0; + if(fSumd02HistLS[i])fSumd02HistLS[i]=0; + if(fSigVertHistLS[i])fSigVertHistLS[i]=0; + if(fPtMaxHistLS[i])fPtMaxHistLS[i]=0; + if(fDCAHistLS[i])fDCAHistLS[i]=0; + if(fMassHistLSTC[i])fMassHistLSTC[i]=0; + } + + + // Default constructor - // Output slot #1 writes into a TList container + // Output slot #1 writes into a TList container DefineOutput(1,TList::Class()); //My private output // Output slot #2 writes cut to private output // DefineOutput(2,AliRDHFCutsDplustoKpipi::Class()); DefineOutput(2,TList::Class()); -// Output slot #3 writes cut to private output + // Output slot #3 writes cut to private output DefineOutput(3,AliNormalizationCounter::Class()); - + if(fFillNtuple){ // Output slot #4 writes into a TNtuple container DefineOutput(4,TNtuple::Class()); //My private output @@ -136,6 +185,10 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus() delete fHistNEvents; fHistNEvents=0; } + for(Int_t i=0;i<3;i++){ + if(fHistCentrality[i]){delete fHistCentrality[i]; fHistCentrality[i]=0;} + } + for(Int_t i=0;i<3*fNPtBins;i++){ if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;} if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;} @@ -143,11 +196,18 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus() if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;} if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;} if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;} + if(fPtKHist[i]){ delete fPtKHist[i]; fPtKHist[i]=0;} + if(fPtpi1Hist[i]){ delete fPtpi1Hist[i]; fPtpi1Hist[i]=0;} + if(fPtpi2Hist[i]){ delete fPtpi2Hist[i]; fPtpi2Hist[i]=0;} if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;} if(fMassHistTC[i]){ delete fMassHistTC[i]; fMassHistTC[i]=0;} if(fMassHistTCPlus[i]){ delete fMassHistTCPlus[i]; fMassHistTCPlus[i]=0;} if(fMassHistTCMinus[i]){ delete fMassHistTCMinus[i]; fMassHistTCMinus[i]=0;} + if(fDLxy[i]){delete fDLxy[i]; fDLxy[i]=0;} + if(fDLxyTC[i]){delete fDLxyTC[i]; fDLxyTC[i]=0;} + if(fCosxy[i]){delete fCosxy[i]; fCosxy[i]=0;} + if(fCosxyTC[i]){delete fCosxyTC[i]; fCosxyTC[i]=0;} if(fMassHistLS[i]){ delete fMassHistLS[i]; fMassHistLS[i]=0;} if(fCosPHistLS[i]){ delete fCosPHistLS[i]; fCosPHistLS[i]=0;} if(fDLenHistLS[i]){ delete fDLenHistLS[i]; fDLenHistLS[i]=0;} @@ -157,6 +217,11 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus() if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;} if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;} } + + for(Int_t i=0;i<3;i++){ + if(fCorreld0Kd0pi[i]){ delete fCorreld0Kd0pi[i]; fCorreld0Kd0pi[i]=0;} + } + if(fPtVsMass){ delete fPtVsMass; fPtVsMass=0; @@ -197,7 +262,9 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus() delete fRDCutsAnalysis; fRDCutsAnalysis = 0; } - + for(Int_t i=0; i<5; i++){ + delete fHistMassPtImpParTC[i]; + } if(fCounter){ delete fCounter; fCounter = 0; @@ -224,34 +291,6 @@ void AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){ SetBinWidth(bw); } } -//________________________________________________________________________ -void AliAnalysisTaskSEDplus::SetPtBinLimit(Int_t n, Float_t* lim){ - // define pt bins for analysis - if(n>kMaxPtBins){ - printf("Max. number of Pt bins = %d\n",kMaxPtBins); - fNPtBins=kMaxPtBins; - fArrayBinLimits[0]=0.; - fArrayBinLimits[1]=2.; - fArrayBinLimits[2]=3.; - fArrayBinLimits[3]=5.; - for(Int_t i=4; ifArrayBinLimits[i-1]){ - fArrayBinLimits[i]=lim[i]; - } - else { - fArrayBinLimits[i]=fArrayBinLimits[i-1]; - } - for(Int_t i=fNPtBins; i 1){ - printf("Number of Pt bins = %d\n",fNPtBins); - for(Int_t i=0; ifNPtBins)return -1; - return fArrayBinLimits[ibin]; -} -//_________________________________________________________________ Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){ return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5); } @@ -283,136 +316,103 @@ void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesA // // Fill the Like Sign histograms // - - //count pos/neg tracks + if(fDebug>1)printf("started LS\n"); + + //histograms for like sign + Int_t nbins=GetNBinsHistos();; + TH1F *histLSPlus = new TH1F("LSPlus","LSPlus",nbins,fLowmasslimit,fUpmasslimit); + TH1F *histLSMinus = new TH1F("LSMinus","LSMinus",nbins,fLowmasslimit,fUpmasslimit); + Int_t nPosTrks=0,nNegTrks=0; - //counter for particles passing single particle cuts - Int_t nspcplus=0; - Int_t nspcminus=0; - - for(Int_t it=0;itGetNumberOfTracks();it++) { - AliAODTrack *track = aod->GetTrack(it); - if(track->Charge()>0){ - nPosTrks++; - if(track->Pt()>=0.4){ - nspcplus++; - } - } - if(track->Charge()<0) - { - nNegTrks++; - if(track->Pt()>=0.4){ - nspcminus++; - } - } - } - Int_t nOStriplets = arrayOppositeSign->GetEntriesFast(); - Int_t nDplusLS=0; + Int_t nDminusLS=0; Int_t nLikeSign = arrayLikeSign->GetEntriesFast(); - Int_t index; - + Int_t index=0; + + // loop over like sign candidates for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) { AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign); + if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts))continue; Bool_t unsetvtx=kFALSE; if(!d->GetOwnPrimaryVtx()) { d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables unsetvtx=kTRUE; } - if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate))nDplusLS++; - if(unsetvtx) d->UnsetOwnPrimaryVtx(); - } - - Float_t wei2=0; - if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign; - Float_t wei3=0; - if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)nDplusLS; - - // loop over like sign candidates - for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) { - AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign); - Bool_t unsetvtx=kFALSE; - if(!d->GetOwnPrimaryVtx()) { - d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - unsetvtx=kTRUE; + Double_t ptCand = d->Pt(); + Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand); + if(iPtBin<0)continue; + Int_t sign= d->GetCharge(); + if(sign>0){ + nPosTrks++; + }else{ + nNegTrks++; } - - if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate)){ - - //set tight cuts values - Int_t iPtBin=-1; - Double_t ptCand = d->Pt(); - for(Int_t ibin=0;ibinfArrayBinLimits[0]&&ptCandIsSelected(d,AliRDHFCuts::kCandidate); - - Int_t sign= d->GetCharge(); - Float_t wei=1; - Float_t wei4=1; - if(sign>0&&nPosTrks>2&&nspcplus>2) { //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event - - wei=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.); - wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.); - } - - if(sign<0&&nNegTrks>2&&nspcminus>2){ - wei=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.); - wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.); - - } - + // if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)){ + fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod); + Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts(); + if(passTightCuts>0){ Float_t invMass = d->InvMassDplus(); - Double_t dlen=d->DecayLength(); - Double_t cosp=d->CosPointingAngle(); - Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2); - Double_t dca=d->GetDCA(); - Double_t sigvert=d->GetSigmaVert(); - Double_t ptmax=0; - for(Int_t i=0;i<3;i++){ - if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i); - } - index=GetLSHistoIndex(iPtBin); - fMassHistLS[index]->Fill(invMass,wei); fMassHistLS[index+1]->Fill(invMass); - fMassHistLS[index+2]->Fill(invMass,wei2); - fMassHistLS[index+3]->Fill(invMass,wei3); - fMassHistLS[index+4]->Fill(invMass,wei4); - - Int_t indexcut=GetHistoIndex(iPtBin); - fCosPHistLS[indexcut]->Fill(cosp); - fDLenHistLS[indexcut]->Fill(dlen); - fSumd02HistLS[indexcut]->Fill(sumD02); - fSigVertHistLS[indexcut]->Fill(sigvert); - fPtMaxHistLS[indexcut]->Fill(ptmax); - fDCAHistLS[indexcut]->Fill(dca); - - if(passTightCuts){ - fMassHistLSTC[index]->Fill(invMass,wei); - fMassHistLSTC[index+1]->Fill(invMass); - fMassHistLSTC[index+2]->Fill(invMass,wei2); - fMassHistLSTC[index+3]->Fill(invMass,wei3); - fMassHistLSTC[index+4]->Fill(invMass,wei4); + if(sign>0){ + histLSPlus->Fill(invMass); + nDplusLS++; + }else{ + histLSMinus->Fill(invMass); + nDminusLS++; + } + if(fCutsDistr){ + Double_t dlen=d->DecayLength(); + Double_t cosp=d->CosPointingAngle(); + Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2); + Double_t dca=d->GetDCA(); + Double_t sigvert=d->GetSigmaVert(); + Double_t ptmax=0; + for(Int_t i=0;i<3;i++){ + if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i); + } + fCosPHistLS[iPtBin]->Fill(cosp); + fDLenHistLS[iPtBin]->Fill(dlen); + fSumd02HistLS[iPtBin]->Fill(sumD02); + fSigVertHistLS[iPtBin]->Fill(sigvert); + fPtMaxHistLS[iPtBin]->Fill(ptmax); + fDCAHistLS[iPtBin]->Fill(dca); } } if(unsetvtx) d->UnsetOwnPrimaryVtx(); } + //wei:normalized to the number of combinations (production) + //wei2:normalized to the number of LS/OS (production) + //wei3:normalized to the number of LS/OS (analysis) + //wei4:normalized to the number of combinations (analysis) + Float_t wei2=0; + if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign; + Float_t wei3=0; + if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)(nDplusLS+nDminusLS); + Float_t weiplus=1.,weiminus=1.; + Float_t wei4plus=1.,wei4minus=1.; + //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event + if(nPosTrks>2)weiplus=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.); + if(nDplusLS>2)wei4plus=3.*(Float_t)nDminusLS/((Float_t)nDplusLS-2.); + if(nNegTrks>2)weiminus=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.); + if(nDminusLS>2)wei4minus=3.*(Float_t)nDplusLS/((Float_t)nDminusLS-2.); + + fMassHistLS[index]->Add(histLSPlus,weiplus); + fMassHistLS[index]->Add(histLSMinus,weiminus); + fMassHistLS[index+2]->Add(histLSPlus,wei2); + fMassHistLS[index+2]->Add(histLSMinus,wei2); + fMassHistLS[index+3]->Add(histLSPlus,wei3); + fMassHistLS[index+3]->Add(histLSMinus,wei3); + fMassHistLS[index+4]->Add(histLSPlus,wei4plus); + fMassHistLS[index+4]->Add(histLSMinus,wei4minus); + + delete histLSPlus;histLSPlus=0; + delete histLSMinus;histLSMinus=0; - //printf("------------ N. of positive tracks in Event ----- %d \n", nPosTrks); - //printf("------------ N. of negative tracks in Event ----- %d \n", nNegTrks); - - // printf("LS analysis...done\n"); - + if(fDebug>1) printf("LS analysis terminated\n"); } - //__________________________________________ void AliAnalysisTaskSEDplus::Init(){ // @@ -422,10 +422,10 @@ void AliAnalysisTaskSEDplus::Init(){ //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1 fListCuts=new TList(); - AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi(); - production=fRDCutsProduction; - AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(); - analysis=fRDCutsAnalysis; + AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi(*fRDCutsProduction); + production->SetName("ProductionCuts"); + AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis); + analysis->SetName("AnalysisCuts"); fListCuts->Add(production); fListCuts->Add(analysis); @@ -448,12 +448,17 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects() TString hisname; Int_t index=0; - Int_t indexLS=0; Int_t nbins=GetNBinsHistos(); + fHistCentrality[0]=new TH1F("centrality","centrality",100,0.5,30000.5); + fHistCentrality[1]=new TH1F("centrality(selectedCent)","centrality(selectedCent)",100,0.5,30000.5); + fHistCentrality[2]=new TH1F("centrality(OutofCent)","centrality(OutofCent)",100,0.5,30000.5); + for(Int_t i=0;i<3;i++){ + fHistCentrality[i]->Sumw2(); + fOutput->Add(fHistCentrality[i]); + } for(Int_t i=0;iSumw2(); - + hisname.Form("hPtKPt%d",i); + fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtKHist[index]->Sumw2(); + hisname.Form("hPtpi1Pt%d",i); + fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtpi1Hist[index]->Sumw2(); + hisname.Form("hPtpi2Pt%d",i); + fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtpi2Hist[index]->Sumw2(); hisname.Form("hDCAAllPt%d",i); fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); fDCAHist[index]->Sumw2(); - - + hisname.Form("hDLxyPt%d",i); + fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.); + fDLxy[index]->Sumw2(); + hisname.Form("hCosxyPt%d",i); + fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.); + fCosxy[index]->Sumw2(); + hisname.Form("hDLxyPt%dTC",i); + fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.); + fDLxyTC[index]->Sumw2(); + hisname.Form("hCosxyPt%dTC",i); + fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.); + fCosxyTC[index]->Sumw2(); + hisname.Form("hMassPt%dTC",i); fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); fMassHistTC[index]->Sumw2(); @@ -491,40 +515,8 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects() fMassHistTCMinus[index]->Sumw2(); - - - hisname.Form("hCosPAllPt%dLS",i); - fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.); - fCosPHistLS[index]->Sumw2(); - hisname.Form("hDLenAllPt%dLS",i); - fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5); - fDLenHistLS[index]->Sumw2(); - hisname.Form("hSumd02AllPt%dLS",i); - fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.); - fSumd02HistLS[index]->Sumw2(); - hisname.Form("hSigVertAllPt%dLS",i); - fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); - fSigVertHistLS[index]->Sumw2(); - hisname.Form("hPtMaxAllPt%dLS",i); - fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); - fPtMaxHistLS[index]->Sumw2(); - - hisname.Form("hDCAAllPt%dLS",i); - fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); - fDCAHistLS[index]->Sumw2(); - - hisname.Form("hLSPt%dLC",i); - fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLS[indexLS]->Sumw2(); - - hisname.Form("hLSPt%dTC",i); - fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLSTC[indexLS]->Sumw2(); - - - + index=GetSignalHistoIndex(i); - indexLS++; hisname.Form("hSigPt%d",i); fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); fMassHist[index]->Sumw2(); @@ -542,13 +534,33 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects() fSigVertHist[index]->Sumw2(); hisname.Form("hPtMaxSigPt%d",i); fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); - fPtMaxHist[index]->Sumw2(); + fPtMaxHist[index]->Sumw2(); + hisname.Form("hPtKSigPt%d",i); + fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtKHist[index]->Sumw2(); + hisname.Form("hPtpi1SigPt%d",i); + fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtpi1Hist[index]->Sumw2(); + hisname.Form("hPtpi2SigPt%d",i); + fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtpi2Hist[index]->Sumw2(); hisname.Form("hDCASigPt%d",i); fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); fDCAHist[index]->Sumw2(); - + hisname.Form("hDLxySigPt%d",i); + fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.); + fDLxy[index]->Sumw2(); + hisname.Form("hCosxySigPt%d",i); + fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.); + fCosxy[index]->Sumw2(); + hisname.Form("hDLxySigPt%dTC",i); + fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.); + fDLxyTC[index]->Sumw2(); + hisname.Form("hCosxySigPt%dTC",i); + fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.); + fCosxyTC[index]->Sumw2(); hisname.Form("hSigPt%dTC",i); fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); fMassHistTC[index]->Sumw2(); @@ -559,39 +571,8 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects() fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); fMassHistTCMinus[index]->Sumw2(); - hisname.Form("hLSPt%dLCnw",i); - fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLS[indexLS]->Sumw2(); - hisname.Form("hLSPt%dTCnw",i); - fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLSTC[indexLS]->Sumw2(); - - - - hisname.Form("hCosPSigPt%dLS",i); - fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.); - fCosPHistLS[index]->Sumw2(); - hisname.Form("hDLenSigPt%dLS",i); - fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5); - fDLenHistLS[index]->Sumw2(); - hisname.Form("hSumd02SigPt%dLS",i); - fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.); - fSumd02HistLS[index]->Sumw2(); - hisname.Form("hSigVertSigPt%dLS",i); - fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); - fSigVertHistLS[index]->Sumw2(); - hisname.Form("hPtMaxSigPt%dLS",i); - fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); - fPtMaxHistLS[index]->Sumw2(); - - hisname.Form("hDCASigPt%dLS",i); - fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); - fDCAHistLS[index]->Sumw2(); - - index=GetBackgroundHistoIndex(i); - indexLS++; hisname.Form("hBkgPt%d",i); fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); fMassHist[index]->Sumw2(); @@ -610,11 +591,32 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects() hisname.Form("hPtMaxBkgPt%d",i); fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); fPtMaxHist[index]->Sumw2(); - + hisname.Form("hPtKBkgPt%d",i); + fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtKHist[index]->Sumw2(); + hisname.Form("hPtpi1BkgPt%d",i); + fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtpi1Hist[index]->Sumw2(); + hisname.Form("hPtpi2BkgPt%d",i); + fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtpi2Hist[index]->Sumw2(); hisname.Form("hDCABkgPt%d",i); fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); fDCAHist[index]->Sumw2(); + hisname.Form("hDLxyBkgPt%d",i); + fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.); + fDLxy[index]->Sumw2(); + hisname.Form("hCosxyBkgPt%d",i); + fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.); + fCosxy[index]->Sumw2(); + hisname.Form("hDLxyBkgPt%dTC",i); + fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.); + fDLxyTC[index]->Sumw2(); + hisname.Form("hCosxyBkgPt%dTC",i); + fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.); + fCosxyTC[index]->Sumw2(); + hisname.Form("hBkgPt%dTC",i); fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); @@ -625,85 +627,59 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects() hisname.Form("hBkgPt%dTCMinus",i); fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); fMassHistTCMinus[index]->Sumw2(); - - hisname.Form("hLSPt%dLCntrip",i); - fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLS[indexLS]->Sumw2(); - hisname.Form("hLSPt%dTCntrip",i); - fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLSTC[indexLS]->Sumw2(); - - - hisname.Form("hCosPBkgPt%dLS",i); - fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.); - fCosPHistLS[index]->Sumw2(); - hisname.Form("hDLenBkgPt%dLS",i); - fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5); - fDLenHistLS[index]->Sumw2(); - hisname.Form("hSumd02BkgPt%dLS",i); - fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.); - fSumd02HistLS[index]->Sumw2(); - hisname.Form("hSigVertBkgPt%dLS",i); - fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); - fSigVertHistLS[index]->Sumw2(); - hisname.Form("hPtMaxBkgPt%dLS",i); - fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); - fPtMaxHistLS[index]->Sumw2(); - hisname.Form("hDCABkgPt%dLS",i); - fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); - fDCAHistLS[index]->Sumw2(); + } - indexLS++; - hisname.Form("hLSPt%dLCntripsinglecut",i); - fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLS[indexLS]->Sumw2(); - hisname.Form("hLSPt%dTCntripsinglecut",i); - fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLSTC[indexLS]->Sumw2(); - - indexLS++; - hisname.Form("hLSPt%dLCspc",i); - fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLS[indexLS]->Sumw2(); - hisname.Form("hLSPt%dTCspc",i); - fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); - fMassHistLSTC[indexLS]->Sumw2(); - } - for(Int_t i=0; i<3*fNPtBins; i++){ fOutput->Add(fMassHist[i]); - fOutput->Add(fCosPHist[i]); - fOutput->Add(fDLenHist[i]); - fOutput->Add(fSumd02Hist[i]); - fOutput->Add(fSigVertHist[i]); - fOutput->Add(fPtMaxHist[i]); - fOutput->Add(fDCAHist[i]); + if(fCutsDistr){ + fOutput->Add(fCosPHist[i]); + fOutput->Add(fDLenHist[i]); + fOutput->Add(fSumd02Hist[i]); + fOutput->Add(fSigVertHist[i]); + fOutput->Add(fPtMaxHist[i]); + fOutput->Add(fPtKHist[i]); + fOutput->Add(fPtpi1Hist[i]); + fOutput->Add(fPtpi2Hist[i]); + fOutput->Add(fDCAHist[i]); + fOutput->Add(fDLxy[i]); + fOutput->Add(fDLxyTC[i]); + fOutput->Add(fCosxy[i]); + fOutput->Add(fCosxyTC[i]); + } fOutput->Add(fMassHistTC[i]); fOutput->Add(fMassHistTCPlus[i]); fOutput->Add(fMassHistTCMinus[i]); + } - for(Int_t i=0; i<3*fNPtBins&&fDoLS; i++){ - fOutput->Add(fCosPHistLS[i]); - fOutput->Add(fDLenHistLS[i]); - fOutput->Add(fSumd02HistLS[i]); - fOutput->Add(fSigVertHistLS[i]); - fOutput->Add(fPtMaxHistLS[i]); - fOutput->Add(fDCAHistLS[i]); - } - for(Int_t i=0; i<5*fNPtBins&&fDoLS; i++){ - fOutput->Add(fMassHistLS[i]); - fOutput->Add(fMassHistLSTC[i]); + + if(fCutsDistr){ + fCorreld0Kd0pi[0]=new TH2F("hCorreld0Kd0piAll","",100,-0.02,0.02,100,-0.02,0.02); + fCorreld0Kd0pi[1]=new TH2F("hCorreld0Kd0piSig","",100,-0.02,0.02,100,-0.02,0.02); + fCorreld0Kd0pi[2]=new TH2F("hCorreld0Kd0piBkg","",100,-0.02,0.02,100,-0.02,0.02); + for(Int_t i=0; i<3; i++){ + fCorreld0Kd0pi[i]->Sumw2(); + fOutput->Add(fCorreld0Kd0pi[i]); + } } - - fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5); + fHistNEvents = new TH1F("fHistNEvents", "number of events ",9,-0.5,8.5); + fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal"); + fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with good vertex"); + fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with PbPb HM trigger"); + fHistNEvents->GetXaxis()->SetBinLabel(4,"no. of Rejected pileup events"); + fHistNEvents->GetXaxis()->SetBinLabel(5,"no. of candidate"); + fHistNEvents->GetXaxis()->SetBinLabel(6,"no. of D+ after loose cuts"); + fHistNEvents->GetXaxis()->SetBinLabel(7,"no. of D+ after tight cuts"); + fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events"); + fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask"); + fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE); fHistNEvents->Sumw2(); fHistNEvents->SetMinimum(0); fOutput->Add(fHistNEvents); - fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,40,0.,20.); - fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,40,0.,20.); + fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.); + fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.); fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.); fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.); fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.); @@ -718,10 +694,17 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects() //Counter for Normalization - fCounter = new AliNormalizationCounter("NormalizationCounter");//new line + TString normName="NormalizationCounter"; + AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer(); + if(cont)normName=(TString)cont->GetName(); + fCounter = new AliNormalizationCounter(normName.Data()); + fCounter->Init(); + + if(fDoLS) CreateLikeSignHistos(); + if(fDoImpPar) CreateImpactParameterHistos(); if(fFillNtuple){ - OpenFile(2); // 2 is the slot number of the ntuple + OpenFile(4); // 4 is the slot number of the ntuple fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square"); @@ -756,16 +739,16 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/) array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong"); arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong"); } - } else { + } else if(aod) { array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong"); arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong"); } - if(!array3Prong) { + if(!aod || !array3Prong) { printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n"); return; } - if(!arrayLikeSign) { + if(!arrayLikeSign && fDoLS) { printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n"); return; } @@ -774,25 +757,44 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/) // fix for temporary bug in ESDfilter // the AODs with null vertex pointer didn't pass the PhysSel if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return; - fCounter->StoreEvent(aod,fReadMC); + fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC); fHistNEvents->Fill(0); // count event - // Post the data already here + + Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod); + Bool_t isEvSelP=kTRUE; + isEvSelP=fRDCutsProduction->IsEventSelected(aod); // to have proper PID object settings + + Float_t centrality=aod->GetNTracks();//fRDCutsAnalysis->GetCentrality(aod); + fHistCentrality[0]->Fill(centrality); + // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD + TString trigclass=aod->GetFiredTriggerClasses(); + if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(2); + if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(3); + if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(7);fHistCentrality[2]->Fill(centrality);} + + // Post the data already here PostData(1,fOutput); - + if(!isEvSel)return; + + fHistCentrality[1]->Fill(centrality); + fHistNEvents->Fill(1); + TClonesArray *arrayMC=0; AliAODMCHeader *mcHeader=0; // AOD primary vertex AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); // vtx1->Print(); - + TString primTitle = vtx1->GetTitle(); + //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2); + // load MC particles if(fReadMC){ arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()); if(!arrayMC) { printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n"); - // return; + return; } // load MC header @@ -804,40 +806,53 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/) } Int_t n3Prong = array3Prong->GetEntriesFast(); - if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong); + // printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks()); Int_t nOS=0; Int_t index; Int_t pdgDgDplustoKpipi[3]={321,211,211}; + + if(fDoLS>1){ + for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) { + AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong); + if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){ + if(fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod))nOS++; + } + } + }else{ // Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};//TO REMOVE //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts(); Int_t nSelectedloose=0,nSelectedtight=0; for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) { AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong); - - + fHistNEvents->Fill(4); + if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){ + fHistNEvents->Fill(8); + continue; + } Bool_t unsetvtx=kFALSE; if(!d->GetOwnPrimaryVtx()){ d->SetOwnPrimaryVtx(vtx1); unsetvtx=kTRUE; } - if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate)) { + if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) { - - - Int_t iPtBin = -1; Double_t ptCand = d->Pt(); - - for(Int_t ibin=0;ibinfArrayBinLimits[0]&&ptCandPtBin(ptCand); + + Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod); + Bool_t recVtx=kFALSE; + AliAODVertex *origownvtx=0x0; + if(fRDCutsProduction->GetIsPrimaryWithoutDaughters()){ + if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx()); + if(fRDCutsProduction->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE; + else fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx); } - Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate); - - Int_t labDp=-1; + Bool_t isPrimary=kTRUE; Float_t deltaPx=0.; Float_t deltaPy=0.; Float_t deltaPz=0.; @@ -846,10 +861,12 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/) Float_t yDecay=0.; Float_t zDecay=0.; Float_t pdgCode=-2; + Float_t trueImpParXY=0.; if(fReadMC){ labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi); if(labDp>=0){ AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp); + if(CheckOrigin(arrayMC,partDp)==5) isPrimary=kFALSE; AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0)); deltaPx=partDp->Px()-d->Px(); deltaPy=partDp->Py()-d->Py(); @@ -859,21 +876,37 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/) yDecay=dg0->Yv(); zDecay=dg0->Zv(); pdgCode=TMath::Abs(partDp->GetPdgCode()); + if(!isPrimary){ + trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.; + } }else{ pdgCode=-1; } } + Double_t invMass=d->InvMassDplus(); Double_t rapid=d->YDplus(); fYVsPt->Fill(ptCand,rapid); - if(passTightCuts) fYVsPtTC->Fill(ptCand,rapid); + if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid);nOS++;} Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid); if(isFidAcc){ fPtVsMass->Fill(invMass,ptCand); if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand); } Float_t tmp[24]; - if(fFillNtuple){ + Double_t dlen=d->DecayLength(); + Double_t cosp=d->CosPointingAngle(); + Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2); + Double_t dca=d->GetDCA(); + Double_t sigvert=d->GetSigmaVert(); + Double_t ptmax=0; + for(Int_t i=0;i<3;i++){ + if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i); + } + Double_t impparXY=d->ImpParXY()*10000.; + Double_t arrayForSparse[3]={invMass,ptCand,impparXY}; + Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY}; + if(fFillNtuple){ tmp[0]=pdgCode; tmp[1]=deltaPx; tmp[2]=deltaPy; @@ -886,52 +919,61 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/) tmp[9]=d->PtProng(1); tmp[10]=d->PtProng(2); tmp[11]=d->Pt(); - tmp[12]=d->CosPointingAngle(); - tmp[13]=d->DecayLength(); + tmp[12]=cosp; + tmp[13]=dlen; tmp[14]=d->Xv(); tmp[15]=d->Yv(); tmp[16]=d->Zv(); tmp[17]=d->InvMassDplus(); - tmp[18]=d->GetSigmaVert(); + tmp[18]=sigvert; tmp[19]=d->Getd0Prong(0); tmp[20]=d->Getd0Prong(1); tmp[21]=d->Getd0Prong(2); - tmp[22]=d->GetDCA(); + tmp[22]=dca; tmp[23]=d->Prodd0d0(); fNtupleDplus->Fill(tmp); - PostData(3,fNtupleDplus); - } - Double_t dlen=d->DecayLength(); - Double_t cosp=d->CosPointingAngle(); - Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2); - Double_t dca=d->GetDCA(); - Double_t sigvert=d->GetSigmaVert(); - Double_t ptmax=0; - for(Int_t i=0;i<3;i++){ - if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i); + PostData(4,fNtupleDplus); } if(iPtBin>=0){ - + Float_t dlxy=d->NormalizedDecayLengthXY(); + Float_t cxy=d->CosPointingAngleXY(); index=GetHistoIndex(iPtBin); if(isFidAcc){ + fHistNEvents->Fill(5); nSelectedloose++; fMassHist[index]->Fill(invMass); - fCosPHist[index]->Fill(cosp); - fDLenHist[index]->Fill(dlen); - fSumd02Hist[index]->Fill(sumD02); - fSigVertHist[index]->Fill(sigvert); - fPtMaxHist[index]->Fill(ptmax); - fDCAHist[index]->Fill(dca); - - if(passTightCuts){ + if(fCutsDistr){ + fCosPHist[index]->Fill(cosp); + fDLenHist[index]->Fill(dlen); + fSumd02Hist[index]->Fill(sumD02); + fSigVertHist[index]->Fill(sigvert); + fPtMaxHist[index]->Fill(ptmax); + fPtKHist[index]->Fill(d->PtProng(1)); + fPtpi1Hist[index]->Fill(d->PtProng(0)); + fPtpi2Hist[index]->Fill(d->PtProng(2)); + fDCAHist[index]->Fill(dca); + fDLxy[index]->Fill(dlxy); + fCosxy[index]->Fill(cxy); + fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1), + d->Getd0Prong(2)*d->Getd0Prong(1)); + } + if(passTightCuts){ fHistNEvents->Fill(6); nSelectedtight++; fMassHistTC[index]->Fill(invMass); + if(fCutsDistr){ + fDLxyTC[index]->Fill(dlxy); + fCosxyTC[index]->Fill(cxy); + } if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass); else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass); - } + if(fDoImpPar){ + fHistMassPtImpParTC[0]->Fill(arrayForSparse); + } + } } - + if(fReadMC){ + // if(fCutsDistr){ if(labDp>=0) { index=GetSignalHistoIndex(iPtBin); if(isFidAcc){ @@ -964,17 +1006,38 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/) } Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k]; fMassHist[index]->Fill(invMass); - fCosPHist[index]->Fill(cosp,fact); - fDLenHist[index]->Fill(dlen,fact); - Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2]; - fSumd02Hist[index]->Fill(sumd02s); - fSigVertHist[index]->Fill(sigvert,fact); - fPtMaxHist[index]->Fill(ptmax,fact); - fDCAHist[index]->Fill(dca,fact); + if(fCutsDistr){ + fCosPHist[index]->Fill(cosp,fact); + fDLenHist[index]->Fill(dlen,fact); + fDLxy[index]->Fill(dlxy); + fCosxy[index]->Fill(cxy); + + Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2]; + fSumd02Hist[index]->Fill(sumd02s); + fSigVertHist[index]->Fill(sigvert,fact); + fPtMaxHist[index]->Fill(ptmax,fact); + fPtKHist[index]->Fill(d->PtProng(1),fact); + fPtpi1Hist[index]->Fill(d->PtProng(0),fact); + fPtpi2Hist[index]->Fill(d->PtProng(2),fact); + fDCAHist[index]->Fill(dca,fact); + fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1), + d->Getd0Prong(2)*d->Getd0Prong(1)); + } if(passTightCuts){ fMassHistTC[index]->Fill(invMass); + if(fCutsDistr){ + fDLxyTC[index]->Fill(dlxy); + fCosxyTC[index]->Fill(cxy); + } if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass); else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass); + if(fDoImpPar){ + if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse); + else{ + fHistMassPtImpParTC[2]->Fill(arrayForSparse); + fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue); + } + } } } fYVsPtSig->Fill(ptCand,rapid); @@ -1009,39 +1072,217 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/) }//if labd>=0 }//prong loop } + Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k]; fMassHist[index]->Fill(invMass); - fCosPHist[index]->Fill(cosp,fact); - fDLenHist[index]->Fill(dlen,fact); - Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2]; - fSumd02Hist[index]->Fill(sumd02s); - fSigVertHist[index]->Fill(sigvert,fact); - fPtMaxHist[index]->Fill(ptmax,fact); - fDCAHist[index]->Fill(dca,fact); + if(fCutsDistr){ + fCosPHist[index]->Fill(cosp,fact); + fDLenHist[index]->Fill(dlen,fact); + fDLxy[index]->Fill(dlxy); + fCosxy[index]->Fill(cxy); + + Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2]; + fSumd02Hist[index]->Fill(sumd02s); + fSigVertHist[index]->Fill(sigvert,fact); + fPtMaxHist[index]->Fill(ptmax,fact); + fPtKHist[index]->Fill(d->PtProng(1),fact); + fPtpi1Hist[index]->Fill(d->PtProng(0),fact); + fPtpi2Hist[index]->Fill(d->PtProng(2),fact); + fDCAHist[index]->Fill(dca,fact); + fCorreld0Kd0pi[2]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1), + d->Getd0Prong(2)*d->Getd0Prong(1)); + } if(passTightCuts){ fMassHistTC[index]->Fill(invMass); + if(fCutsDistr){ + fDLxyTC[index]->Fill(dlxy); + fCosxyTC[index]->Fill(cxy); + } if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass); else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass); + if(fDoImpPar){ + fHistMassPtImpParTC[4]->Fill(arrayForSparse); + } } - } + } } + } - } + } + + if(recVtx)fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx); } if(unsetvtx) d->UnsetOwnPrimaryVtx(); } fCounter->StoreCandidates(aod,nSelectedloose,kTRUE); fCounter->StoreCandidates(aod,nSelectedtight,kFALSE); - + } //start LS analysis if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS); - PostData(1,fOutput); + PostData(1,fOutput); + PostData(2,fListCuts); PostData(3,fCounter); return; } +//________________________________________________________________________ +void AliAnalysisTaskSEDplus::CreateLikeSignHistos(){ + // Histos for Like Sign bckground + + TString hisname; + Int_t indexLS=0; + Int_t index=0; + Int_t nbins=GetNBinsHistos(); + for(Int_t i=0;iSumw2(); + hisname.Form("hLSPt%dTC",i); + fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistLSTC[indexLS]->Sumw2(); + + hisname.Form("hCosPAllPt%dLS",i); + fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.); + fCosPHistLS[index]->Sumw2(); + hisname.Form("hDLenAllPt%dLS",i); + fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5); + fDLenHistLS[index]->Sumw2(); + hisname.Form("hSumd02AllPt%dLS",i); + fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.); + fSumd02HistLS[index]->Sumw2(); + hisname.Form("hSigVertAllPt%dLS",i); + fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fSigVertHistLS[index]->Sumw2(); + hisname.Form("hPtMaxAllPt%dLS",i); + fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtMaxHistLS[index]->Sumw2(); + hisname.Form("hDCAAllPt%dLS",i); + fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fDCAHistLS[index]->Sumw2(); + + index=GetSignalHistoIndex(i); + indexLS++; + + hisname.Form("hLSPt%dLCnw",i); + fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistLS[indexLS]->Sumw2(); + hisname.Form("hLSPt%dTCnw",i); + fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistLSTC[indexLS]->Sumw2(); + + hisname.Form("hCosPSigPt%dLS",i); + fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.); + fCosPHistLS[index]->Sumw2(); + hisname.Form("hDLenSigPt%dLS",i); + fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5); + fDLenHistLS[index]->Sumw2(); + hisname.Form("hSumd02SigPt%dLS",i); + fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.); + fSumd02HistLS[index]->Sumw2(); + hisname.Form("hSigVertSigPt%dLS",i); + fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fSigVertHistLS[index]->Sumw2(); + hisname.Form("hPtMaxSigPt%dLS",i); + fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtMaxHistLS[index]->Sumw2(); + hisname.Form("hDCASigPt%dLS",i); + fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fDCAHistLS[index]->Sumw2(); + + index=GetBackgroundHistoIndex(i); + indexLS++; + + hisname.Form("hLSPt%dLCntrip",i); + fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistLS[indexLS]->Sumw2(); + hisname.Form("hLSPt%dTCntrip",i); + fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistLSTC[indexLS]->Sumw2(); + + hisname.Form("hCosPBkgPt%dLS",i); + fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.); + fCosPHistLS[index]->Sumw2(); + hisname.Form("hDLenBkgPt%dLS",i); + fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5); + fDLenHistLS[index]->Sumw2(); + hisname.Form("hSumd02BkgPt%dLS",i); + fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.); + fSumd02HistLS[index]->Sumw2(); + hisname.Form("hSigVertBkgPt%dLS",i); + fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fSigVertHistLS[index]->Sumw2(); + hisname.Form("hPtMaxBkgPt%dLS",i); + fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtMaxHistLS[index]->Sumw2(); + hisname.Form("hDCABkgPt%dLS",i); + fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fDCAHistLS[index]->Sumw2(); + + indexLS++; + hisname.Form("hLSPt%dLCntripsinglecut",i); + fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistLS[indexLS]->Sumw2(); + hisname.Form("hLSPt%dTCntripsinglecut",i); + fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistLSTC[indexLS]->Sumw2(); + + indexLS++; + hisname.Form("hLSPt%dLCspc",i); + fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistLS[indexLS]->Sumw2(); + hisname.Form("hLSPt%dTCspc",i); + fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistLSTC[indexLS]->Sumw2(); + } + + for(Int_t i=0; i<3*fNPtBins; i++){ + fOutput->Add(fCosPHistLS[i]); + fOutput->Add(fDLenHistLS[i]); + fOutput->Add(fSumd02HistLS[i]); + fOutput->Add(fSigVertHistLS[i]); + fOutput->Add(fPtMaxHistLS[i]); + fOutput->Add(fDCAHistLS[i]); + } + for(Int_t i=0; i<5*fNPtBins; i++){ + fOutput->Add(fMassHistLS[i]); + fOutput->Add(fMassHistLSTC[i]); + } +} +//________________________________________________________________________ +void AliAnalysisTaskSEDplus::CreateImpactParameterHistos(){ + // Histos for impact paramter study + + Int_t nmassbins=GetNBinsHistos(); + Int_t nbins[3]={nmassbins,200,fNImpParBins}; + Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar}; + Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar}; + + fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll", + "Mass vs. pt vs.imppar - All", + 3,nbins,xmin,xmax); + fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt", + "Mass vs. pt vs.imppar - promptD", + 3,nbins,xmin,xmax); + fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed", + "Mass vs. pt vs.imppar - DfromB", + 3,nbins,xmin,xmax); + fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed", + "Mass vs. pt vs.true imppar -DfromB", + 3,nbins,xmin,xmax); + fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg", + "Mass vs. pt vs.imppar - backgr.", + 3,nbins,xmin,xmax); + + for(Int_t i=0; i<5;i++){ + fOutput->Add(fHistMassPtImpParTC[i]); + } +} //________________________________________________________________________ void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/) @@ -1056,172 +1297,17 @@ void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/) return; } fHistNEvents = dynamic_cast(fOutput->FindObject("fHistNEvents")); - fYVsPt = dynamic_cast(fOutput->FindObject("hYVsPt")); - fYVsPtTC = dynamic_cast(fOutput->FindObject("hYVsPtTC")); - fYVsPtSig = dynamic_cast(fOutput->FindObject("hYVsPtSig")); - fYVsPtSigTC = dynamic_cast(fOutput->FindObject("hYVsPtSigTC")); - fPtVsMass = dynamic_cast(fOutput->FindObject("hPtVsMass")); - fPtVsMassTC = dynamic_cast(fOutput->FindObject("hPtVsMassTC")); TString hisname; Int_t index=0; - - Int_t indexLS=0; for(Int_t i=0;i(fOutput->FindObject(hisname.Data())); - hisname.Form("hCosPAllPt%d",i); - fCosPHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDLenAllPt%d",i); - fDLenHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSumd02AllPt%d",i); - fSumd02Hist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSigVertAllPt%d",i); - fSigVertHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hPtMaxAllPt%d",i); - fPtMaxHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDCAAllPt%d",i); - fDCAHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hMassPt%dTC",i); - fMassHistTC[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hMassPt%dTCPlus",i); - fMassHistTCPlus[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hMassPt%dTCMinus",i); - fMassHistTCMinus[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - if(fDoLS){ - hisname.Form("hLSPt%dLC",i); - fMassHistLS[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hCosPAllPt%dLS",i); - fCosPHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDLenAllPt%dLS",i); - fDLenHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSumd02AllPt%dLS",i); - fSumd02HistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSigVertAllPt%dLS",i); - fSigVertHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hPtMaxAllPt%dLS",i); - fPtMaxHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDCAAllPt%dLS",i); - fDCAHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - - hisname.Form("hLSPt%dTC",i); - fMassHistLSTC[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - - } - - index=GetSignalHistoIndex(i); - if(fDoLS)indexLS++; - hisname.Form("hSigPt%d",i); - fMassHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hCosPSigPt%d",i); - fCosPHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDLenSigPt%d",i); - fDLenHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSumd02SigPt%d",i); - fSumd02Hist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSigVertSigPt%d",i); - fSigVertHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hPtMaxSigPt%d",i); - fPtMaxHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDCASigPt%d",i); - fDCAHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSigPt%dTC",i); + hisname.Form("hMassPt%dTC",i); fMassHistTC[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSigPt%dTCPlus",i); - fMassHistTCPlus[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSigPt%dTCMinus",i); - fMassHistTCMinus[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - if(fDoLS){ - hisname.Form("hLSPt%dLCnw",i); - fMassHistLS[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hCosPSigPt%dLS",i); - fCosPHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDLenSigPt%dLS",i); - fDLenHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSumd02SigPt%dLS",i); - fSumd02HistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSigVertSigPt%dLS",i); - fSigVertHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hPtMaxSigPt%dLS",i); - fPtMaxHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDCASigPt%dLS",i); - fDCAHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - - hisname.Form("hLSPt%dTCnw",i); - fMassHistLSTC[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - } + } - index=GetBackgroundHistoIndex(i); - if(fDoLS)indexLS++; - hisname.Form("hBkgPt%d",i); - fMassHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hCosPBkgPt%d",i); - fCosPHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDLenBkgPt%d",i); - fDLenHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSumd02BkgPt%d",i); - fSumd02Hist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSigVertBkgPt%d",i); - fSigVertHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hPtMaxBkgPt%d",i); - fPtMaxHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDCABkgPt%d",i); - fDCAHist[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hBkgPt%dTC",i); - fMassHistTC[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hBkgPt%dTCPlus",i); - fMassHistTCPlus[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hBkgPt%dTCMinus",i); - fMassHistTCMinus[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - if(fDoLS){ - hisname.Form("hLSPt%dLCntrip",i); - fMassHistLS[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - - - hisname.Form("hCosPBkgPt%dLS",i); - fCosPHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDLenBkgPt%dLS",i); - fDLenHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSumd02BkgPt%dLS",i); - fSumd02HistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hSigVertBkgPt%dLS",i); - fSigVertHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hPtMaxBkgPt%dLS",i); - fPtMaxHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - hisname.Form("hDCABkgPt%dLS",i); - fDCAHistLS[index]=dynamic_cast(fOutput->FindObject(hisname.Data())); - - hisname.Form("hLSPt%dTCntrip",i); - fMassHistLSTC[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - - - indexLS++; - hisname.Form("hLSPt%dLCntripsinglecut",i); - fMassHistLS[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - - hisname.Form("hLSPt%dTCntripsinglecut",i); - fMassHistLSTC[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - - - indexLS++; - hisname.Form("hLSPt%dLCspc",i); - fMassHistLS[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - - hisname.Form("hLSPt%dTCspc",i); - fMassHistLSTC[indexLS]=dynamic_cast(fOutput->FindObject(hisname.Data())); - } - - } - fCounter = dynamic_cast(GetOutputData(3)); - - if(fFillNtuple){ - fNtupleDplus = dynamic_cast(GetOutputData(4)); - } - TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500); c1->cd(); TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC"); @@ -1231,3 +1317,115 @@ void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/) return; } +//_________________________________________________________________________________________________ +Int_t AliAnalysisTaskSEDplus::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const { + // + // checking whether the mother of the particles come from a charm or a bottom quark + // + + Int_t pdgGranma = 0; + Int_t mother = 0; + mother = mcPartCandidate->GetMother(); + Int_t istep = 0; + Int_t abspdgGranma =0; + Bool_t isFromB=kFALSE; + Bool_t isQuarkFound=kFALSE; + while (mother >0 ){ + istep++; + AliAODMCParticle* mcGranma = dynamic_cast(arrayMC->At(mother)); + if (mcGranma){ + pdgGranma = mcGranma->GetPdgCode(); + abspdgGranma = TMath::Abs(pdgGranma); + if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){ + isFromB=kTRUE; + } + if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE; + mother = mcGranma->GetMother(); + }else{ + AliError("Failed casting the mother particle!"); + break; + } + } + + if(isFromB) return 5; + else return 4; +} +//_________________________________________________________________________________________________ +Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partDp) const { + // true impact parameter calculation + + Double_t vtxTrue[3]; + mcHeader->GetVertex(vtxTrue); + Double_t origD[3]; + partDp->XvYvZv(origD); + Short_t charge=partDp->Charge(); + Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3]; + for(Int_t iDau=0; iDau<3; iDau++){ + pXdauTrue[iDau]=0.; + pYdauTrue[iDau]=0.; + pZdauTrue[iDau]=0.; + } + + Int_t nDau=partDp->GetNDaughters(); + Int_t labelFirstDau = partDp->GetDaughter(0); + if(nDau==3){ + for(Int_t iDau=0; iDau<3; iDau++){ + Int_t ind = labelFirstDau+iDau; + AliAODMCParticle* part = dynamic_cast(arrayMC->At(ind)); + if(!part){ + AliError("Daughter particle not found in MC array"); + return 99999.; + } + pXdauTrue[iDau]=part->Px(); + pYdauTrue[iDau]=part->Py(); + pZdauTrue[iDau]=part->Pz(); + } + }else if(nDau==2){ + Int_t theDau=0; + for(Int_t iDau=0; iDau<2; iDau++){ + Int_t ind = labelFirstDau+iDau; + AliAODMCParticle* part = dynamic_cast(arrayMC->At(ind)); + if(!part){ + AliError("Daughter particle not found in MC array"); + return 99999.; + } + Int_t pdgCode=TMath::Abs(part->GetPdgCode()); + if(pdgCode==211 || pdgCode==321){ + pXdauTrue[theDau]=part->Px(); + pYdauTrue[theDau]=part->Py(); + pZdauTrue[theDau]=part->Pz(); + ++theDau; + }else{ + Int_t nDauRes=part->GetNDaughters(); + if(nDauRes==2){ + Int_t labelFirstDauRes = part->GetDaughter(0); + for(Int_t iDauRes=0; iDauRes<2; iDauRes++){ + Int_t indDR = labelFirstDauRes+iDauRes; + AliAODMCParticle* partDR = dynamic_cast(arrayMC->At(indDR)); + if(!partDR){ + AliError("Daughter particle not found in MC array"); + return 99999.; + } + + Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode()); + if(pdgCodeDR==211 || pdgCodeDR==321){ + pXdauTrue[theDau]=partDR->Px(); + pYdauTrue[theDau]=partDR->Py(); + pZdauTrue[theDau]=partDR->Pz(); + ++theDau; + } + } + } + } + } + if(theDau!=3){ + AliError("Wrong number of decay prongs"); + return 99999.; + } + } + + Double_t d0dummy[3]={0.,0.,0.}; + AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy); + return aodDplusMC.ImpParXY(); + +}