X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG3%2FvertexingHF%2FAliAnalysisTaskSEDplus.cxx;h=36dad4b714e0bffc9072881b2c39bfe04d8c2938;hb=9daf22f89a6ff15396be26c8f46d2ad43910561e;hp=bde05937b4655ca879014ad5740d1f387f0b70c0;hpb=fc8d975b37e0b4deb892193b79c96c6a61a2b9fb;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx b/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx index bde05937b46..36dad4b714e 100644 --- a/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx +++ b/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx @@ -13,17 +13,27 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -///////////////////////////////////////////////////////////// -// -// AliAnalysisTaskSE for the extraction of signal(e.g D+) of heavy flavor -// decay candidates with the MC truth. - +//************************************************************************* +// Class AliAnalysisTaskSEDplus +// AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and +//comparison of heavy-flavour decay candidates +// to MC truth (kinematics stored in the AOD) +// Authors: Renu Bala, bala@to.infn.it +// F. Prino, prino@to.infn.it +// G. Ortona, ortona@to.infn.it ///////////////////////////////////////////////////////////// #include #include +#include #include +#include #include +#include + +#include "AliAnalysisManager.h" +#include "AliRDHFCutsDplustoKpipi.h" +#include "AliAODHandler.h" #include "AliAODEvent.h" #include "AliAODVertex.h" #include "AliAODTrack.h" @@ -33,6 +43,7 @@ #include "AliAnalysisVertexingHF.h" #include "AliAnalysisTaskSE.h" #include "AliAnalysisTaskSEDplus.h" +#include "AliNormalizationCounter.h" ClassImp(AliAnalysisTaskSEDplus) @@ -40,65 +51,386 @@ ClassImp(AliAnalysisTaskSEDplus) //________________________________________________________________________ AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(): AliAnalysisTaskSE(), -fOutput(0), -fNtupleDplus(0), -fNtupleDplusbackg(0), -fHistMass(0), -fHistSignal(0), -fHistBackground(0), -fVHF(0) + 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(0), + fRDCutsAnalysis(0), + fCounter(0), + fFillNtuple(kFALSE), + fReadMC(kFALSE), + fUseStrangeness(kFALSE), + fDoLS(kFALSE) { - // Default constructor - + // Default constructor } //________________________________________________________________________ -AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name): +AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,AliRDHFCutsDplustoKpipi *dpluscutsprod,Bool_t fillNtuple): AliAnalysisTaskSE(name), -fOutput(0), +fOutput(0), +fHistNEvents(0), +fPtVsMass(0), +fPtVsMassTC(0), +fYVsPt(0), +fYVsPtTC(0), +fYVsPtSig(0), +fYVsPtSigTC(0), fNtupleDplus(0), -fNtupleDplusbackg(0), -fHistMass(0), -fHistSignal(0), -fHistBackground(0), -fVHF(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) { + // + // Standrd constructor + // + //Double_t ptlim[5]={0.,2.,3.,5,9999999.}; + //SetPtBinLimit(5, ptlim); + SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits()); // 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 into a TNtuple container - DefineOutput(2,TNtuple::Class()); //My private output //AD - // Output slot #3 writes into a TNtuple container - DefineOutput(3,TNtuple::Class()); //My private output //AD + // Output slot #2 writes cut to private output + // DefineOutput(2,AliRDHFCutsDplustoKpipi::Class()); + DefineOutput(2,TList::Class()); +// 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 + } } //________________________________________________________________________ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus() { + // // Destructor + // if (fOutput) { delete fOutput; fOutput = 0; } - if (fVHF) { - delete fVHF; - fVHF = 0; + if(fHistNEvents){ + delete fHistNEvents; + fHistNEvents=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;} + if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;} + 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(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(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;} + if(fSumd02HistLS[i]){ delete fSumd02HistLS[i]; fSumd02HistLS[i]=0;} + if(fSigVertHistLS[i]){ delete fSigVertHistLS[i]; fSigVertHistLS[i]=0;} + if(fPtMaxHistLS[i]){ delete fPtMaxHistLS[i]; fPtMaxHistLS[i]=0;} + if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;} + if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;} } -} + if(fPtVsMass){ + delete fPtVsMass; + fPtVsMass=0; + } + if(fPtVsMassTC){ + delete fPtVsMassTC; + fPtVsMassTC=0; + } + if(fYVsPt){ + delete fYVsPt; + fYVsPt=0; + } + if(fYVsPtTC){ + delete fYVsPtTC; + fYVsPtTC=0; + } + if(fYVsPtSig){ + delete fYVsPtSig; + fYVsPtSig=0; + } + if(fYVsPtSigTC){ + delete fYVsPtSigTC; + fYVsPtSigTC=0; + } + if(fNtupleDplus){ + delete fNtupleDplus; + fNtupleDplus=0; + } + if (fListCuts) { + delete fListCuts; + fListCuts = 0; + } + if(fRDCutsProduction){ + delete fRDCutsProduction; + fRDCutsProduction = 0; + } + if(fRDCutsAnalysis){ + delete fRDCutsAnalysis; + fRDCutsAnalysis = 0; + } + + if(fCounter){ + delete fCounter; + fCounter = 0; + } + +} +//_________________________________________________________________ +void AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){ + // set invariant mass limits + Float_t bw=GetBinWidth(); + fUpmasslimit = 1.865+range; + fLowmasslimit = 1.865-range; + SetBinWidth(bw); +} +//_________________________________________________________________ +void AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){ + // set invariant mass limits + if(uplimit>lowlimit) + { + Float_t bw=GetBinWidth(); + fUpmasslimit = lowlimit; + fLowmasslimit = uplimit; + SetBinWidth(bw); + } +} //________________________________________________________________________ -void AliAnalysisTaskSEDplus::Init() -{ - // Initialization +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; i1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width); + } + fBinWidth=width; +} +//_________________________________________________________________ +Double_t AliAnalysisTaskSEDplus::GetPtBinLimit(Int_t ibin){ + // get pt bin limit + if(ibin>fNPtBins)return -1; + return fArrayBinLimits[ibin]; +} +//_________________________________________________________________ +Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){ + return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5); +} +//_________________________________________________________________ +void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){ + // + // + // Fill the Like Sign histograms + // - if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n"); + //count pos/neg tracks + 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 nLikeSign = arrayLikeSign->GetEntriesFast(); + Int_t index; + + 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; + } + if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod))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; + } + + if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)){ + + //set tight cuts values + Int_t iPtBin=-1; + Double_t ptCand = d->Pt(); + for(Int_t ibin=0;ibinfArrayBinLimits[0]&&ptCandIsSelected(d,AliRDHFCuts::kCandidate,aod); - gROOT->LoadMacro("ConfigVertexingHF.C"); + 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.); + + } - fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()"); - fVHF->PrintStatus(); + 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(unsetvtx) d->UnsetOwnPrimaryVtx(); + } + + //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"); + +} + +//__________________________________________ +void AliAnalysisTaskSEDplus::Init(){ + // + // Initialization + // + if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n"); + + //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; + + fListCuts->Add(production); + fListCuts->Add(analysis); + PostData(2,fListCuts); + return; } @@ -112,33 +444,303 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects() // Several histograms are more conveniently managed in a TList fOutput = new TList(); fOutput->SetOwner(); + fOutput->SetName("OutputHistos"); + + TString hisname; + Int_t index=0; + Int_t indexLS=0; + Int_t nbins=GetNBinsHistos(); + for(Int_t i=0;iSumw2(); + hisname.Form("hCosPAllPt%d",i); + fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.); + fCosPHist[index]->Sumw2(); + hisname.Form("hDLenAllPt%d",i); + fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5); + fDLenHist[index]->Sumw2(); + hisname.Form("hSumd02AllPt%d",i); + fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.); + fSumd02Hist[index]->Sumw2(); + hisname.Form("hSigVertAllPt%d",i); + fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fSigVertHist[index]->Sumw2(); + hisname.Form("hPtMaxAllPt%d",i); + fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtMaxHist[index]->Sumw2(); + hisname.Form("hDCAAllPt%d",i); + fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fDCAHist[index]->Sumw2(); - fHistMass->Sumw2(); - fHistSignal->Sumw2(); - fHistBackground->Sumw2(); + hisname.Form("hMassPt%dTC",i); + fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistTC[index]->Sumw2(); + hisname.Form("hMassPt%dTCPlus",i); + fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistTCPlus[index]->Sumw2(); + hisname.Form("hMassPt%dTCMinus",i); + fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistTCMinus[index]->Sumw2(); - //fHistMass->SetMinimum(0); - //fHistSignal->SetMinimum(0); - //fHistBackground->SetMinimum(0); - fOutput->Add(fHistMass); - fOutput->Add(fHistSignal); - fOutput->Add(fHistBackground); + + + 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(); + hisname.Form("hCosPSigPt%d",i); + fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.); + fCosPHist[index]->Sumw2(); + hisname.Form("hDLenSigPt%d",i); + fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5); + fDLenHist[index]->Sumw2(); + hisname.Form("hSumd02SigPt%d",i); + fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.); + fSumd02Hist[index]->Sumw2(); + hisname.Form("hSigVertSigPt%d",i); + fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fSigVertHist[index]->Sumw2(); + hisname.Form("hPtMaxSigPt%d",i); + fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtMaxHist[index]->Sumw2(); + + hisname.Form("hDCASigPt%d",i); + fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fDCAHist[index]->Sumw2(); + + + hisname.Form("hSigPt%dTC",i); + fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistTC[index]->Sumw2(); + hisname.Form("hSigPt%dTCPlus",i); + fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistTCPlus[index]->Sumw2(); + hisname.Form("hSigPt%dTCMinus",i); + 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(); + hisname.Form("hCosPBkgPt%d",i); + fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.); + fCosPHist[index]->Sumw2(); + hisname.Form("hDLenBkgPt%d",i); + fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5); + fDLenHist[index]->Sumw2(); + hisname.Form("hSumd02BkgPt%d",i); + fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.); + fSumd02Hist[index]->Sumw2(); + hisname.Form("hSigVertBkgPt%d",i); + fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fSigVertHist[index]->Sumw2(); + hisname.Form("hPtMaxBkgPt%d",i); + fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.); + fPtMaxHist[index]->Sumw2(); + + hisname.Form("hDCABkgPt%d",i); + fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1); + fDCAHist[index]->Sumw2(); + + + hisname.Form("hBkgPt%dTC",i); + fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistTC[index]->Sumw2(); + hisname.Form("hBkgPt%dTCPlus",i); + fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit); + fMassHistTCPlus[index]->Sumw2(); + 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]); + 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]); + } + + + fHistNEvents = new TH1F("fHistNEvents", "number of events ",7,-0.5,6.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()->SetNdivisions(1,kFALSE); + + fHistNEvents->Sumw2(); + fHistNEvents->SetMinimum(0); + fOutput->Add(fHistNEvents); + + 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.); + fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.); + + fOutput->Add(fPtVsMass); + fOutput->Add(fPtVsMassTC); + fOutput->Add(fYVsPt); + fOutput->Add(fYVsPtTC); + fOutput->Add(fYVsPtSig); + fOutput->Add(fYVsPtSigTC); - fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Ptpi:Ptpi2:PtK:PtRec:PtTrue:PointingAngle:DecLeng:VxTrue:VxRec:InvMass:sigvert"); - fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptpi2bkg:PtKbkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:InvMassbkg:sigvertbkg"); - //fOutput->Add(fNtupleDplus); // AD - //fOutput->Add(fNtupleDplusbackg); //AD + //Counter for Normalization + TString normName="NormalizationCounter"; + AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer(); + if(cont)normName=(TString)cont->GetName(); + fCounter = new AliNormalizationCounter(normName.Data()); + fCounter->SetRejectPileUp(fRDCutsProduction->GetOptPileUp()); + if(fFillNtuple){ + 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"); + + } + return; } @@ -148,93 +750,322 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/) // Execute analysis for current event: // heavy flavor candidates association to MC truth - - AliAODEvent *aod = dynamic_cast (InputEvent()); + + + + TClonesArray *array3Prong = 0; + TClonesArray *arrayLikeSign =0; + if(!aod && AODEvent() && IsStandardAOD()) { + // In case there is an AOD handler writing a standard AOD, use the AOD + // event in memory rather than the input (ESD) event. + aod = dynamic_cast (AODEvent()); + // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) + // have to taken from the AOD event hold by the AliAODExtension + AliAODHandler* aodHandler = (AliAODHandler*) + ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); + if(aodHandler->GetExtensions()) { + AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); + AliAODEvent *aodFromExt = ext->GetAOD(); + array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong"); + arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong"); + } + } else { + array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong"); + arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong"); + } - // load Dplus->Kpipi candidates - TClonesArray *array3Prong = - (TClonesArray*)aod->GetList()->FindObject("Charm3Prong"); if(!array3Prong) { printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n"); return; } - - // AOD primary vertex - AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); - // vtx1->Print(); - - // load MC particles - TClonesArray *arrayMC = - (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()); - if(!arrayMC) { - printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n"); + if(!arrayLikeSign && fDoLS) { + printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n"); return; } + + // 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); + fHistNEvents->Fill(0); // count event + // Post the data already here + if(fRDCutsAnalysis->IsEventSelected(aod))fHistNEvents->Fill(1); + // 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); + + PostData(1,fOutput); + + 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; + } + // load MC header - AliAODMCHeader *mcHeader = - (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()); - if(!mcHeader) { + mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()); + if(!mcHeader) { printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n"); return; + } } - Int_t n3Prong = array3Prong->GetEntriesFast(); - printf("Number of D+->Kpipi: %d\n",n3Prong); + + Int_t n3Prong = array3Prong->GetEntriesFast(); + if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong); + + + Int_t nOS=0; + Int_t index; + Int_t pdgDgDplustoKpipi[3]={321,211,211}; + // 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); + + Bool_t unsetvtx=kFALSE; + if(!d->GetOwnPrimaryVtx()){ + d->SetOwnPrimaryVtx(vtx1); + unsetvtx=kTRUE; + } + if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) { - for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) { - AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong); + + Int_t iPtBin = -1; + Double_t ptCand = d->Pt(); - Bool_t unsetvtx=kFALSE; - if(!d->GetOwnPrimaryVtx()){ - d->SetOwnPrimaryVtx(vtx1); - unsetvtx=kTRUE; + for(Int_t ibin=0;ibinfArrayBinLimits[0]&&ptCandSelectDplus(fVHF->GetDplusCuts())) - { - - - Int_t labDp = d->MatchToMC(411,arrayMC); - // if(labDp>=0) { - AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp); - if(labDp>=0) { - Int_t pdgDp = TMath::Abs(partDp->GetPdgCode()); - if(pdgDp==411){ - - fHistSignal->Fill(d->InvMassDplus()); - fHistMass->Fill(d->InvMassDplus()); - - // Post the data already here - PostData(1,fOutput); - - fNtupleDplus->Fill(pdgDp,partDp->Px()-d->Px(),partDp->Py()-d->Py(),partDp->Pz()-d->Pz(),d->PtProng(0),d->PtProng(2),d->PtProng(1),d->Pt(),partDp->Pt(),d->CosPointingAngle(),d->DecayLength(),partDp->Xv(),d->Xv(),d->InvMassDplus(),d->GetSigmaVert()); - PostData(2,fNtupleDplus); //AD - } - } -else { - fHistBackground->Fill(d->InvMassDplus()); - fNtupleDplusbackg->Fill(d->PtProng(0),d->PtProng(2),d->PtProng(1),d->Pt(),d->CosPointingAngle(),d->DecayLength(),d->Xv(),d->InvMassDplus(),d->GetSigmaVert()); - PostData(3,fNtupleDplusbackg); //AD - - fHistMass->Fill(d->InvMassDplus()); - } + + Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod); + + + Int_t labDp=-1; + Float_t deltaPx=0.; + Float_t deltaPy=0.; + Float_t deltaPz=0.; + Float_t truePt=0.; + Float_t xDecay=0.; + Float_t yDecay=0.; + Float_t zDecay=0.; + Float_t pdgCode=-2; + if(fReadMC){ + labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi); + if(labDp>=0){ + AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp); + AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0)); + deltaPx=partDp->Px()-d->Px(); + deltaPy=partDp->Py()-d->Py(); + deltaPz=partDp->Pz()-d->Pz(); + truePt=partDp->Pt(); + xDecay=dg0->Xv(); + yDecay=dg0->Yv(); + zDecay=dg0->Zv(); + pdgCode=TMath::Abs(partDp->GetPdgCode()); + }else{ + pdgCode=-1; + } + } + Double_t invMass=d->InvMassDplus(); + Double_t rapid=d->YDplus(); + fYVsPt->Fill(ptCand,rapid); + if(passTightCuts) fYVsPtTC->Fill(ptCand,rapid); + 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){ + tmp[0]=pdgCode; + tmp[1]=deltaPx; + tmp[2]=deltaPy; + tmp[3]=deltaPz; + tmp[4]=truePt; + tmp[5]=xDecay; + tmp[6]=yDecay; + tmp[7]=zDecay; + tmp[8]=d->PtProng(0); + tmp[9]=d->PtProng(1); + tmp[10]=d->PtProng(2); + tmp[11]=d->Pt(); + tmp[12]=d->CosPointingAngle(); + tmp[13]=d->DecayLength(); + tmp[14]=d->Xv(); + tmp[15]=d->Yv(); + tmp[16]=d->Zv(); + tmp[17]=d->InvMassDplus(); + tmp[18]=d->GetSigmaVert(); + tmp[19]=d->Getd0Prong(0); + tmp[20]=d->Getd0Prong(1); + tmp[21]=d->Getd0Prong(2); + tmp[22]=d->GetDCA(); + tmp[23]=d->Prodd0d0(); + fNtupleDplus->Fill(tmp); + PostData(4,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); + } + if(iPtBin>=0){ + + 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){ fHistNEvents->Fill(6); + nSelectedtight++; + fMassHistTC[index]->Fill(invMass); + if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass); + else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass); + } } - - - if(unsetvtx) d->UnsetOwnPrimaryVtx(); - + if(fReadMC){ + if(labDp>=0) { + index=GetSignalHistoIndex(iPtBin); + if(isFidAcc){ + Float_t factor[3]={1.,1.,1.}; + if(fUseStrangeness){ + for(Int_t iprong=0;iprong<3;iprong++){ + AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong); + Int_t labd= trad->GetLabel(); + if(labd>=0){ + AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd); + if(dau){ + Int_t labm = dau->GetMother(); + if(labm>=0){ + AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm); + if(mot){ + if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+- + if(d->PtProng(iprong)<=1)factor[iprong]=1./.7; + else factor[iprong]=1./.6; + // fNentries->Fill(12); + } + if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda + factor[iprong]=1./0.25; + // fNentries->Fill(13); + }//if 3122 + }//if(mot) + }//if labm>0 + }//if(dau) + }//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(passTightCuts){ + fMassHistTC[index]->Fill(invMass); + if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass); + else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass); + } + } + fYVsPtSig->Fill(ptCand,rapid); + if(passTightCuts) fYVsPtSigTC->Fill(ptCand,rapid); + }else{ + index=GetBackgroundHistoIndex(iPtBin); + if(isFidAcc){ + Float_t factor[3]={1.,1.,1.}; + if(fUseStrangeness){ + for(Int_t iprong=0;iprong<3;iprong++){ + AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong); + Int_t labd= trad->GetLabel(); + if(labd>=0){ + AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd); + if(dau){ + Int_t labm = dau->GetMother(); + if(labm>=0){ + AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm); + if(mot){ + if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+- + if(d->PtProng(iprong)<=1)factor[iprong]=1./.7; + else factor[iprong]=1./.6; + // fNentries->Fill(12); + } + if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda + factor[iprong]=1./0.25; + // fNentries->Fill(13); + }//if 3122 + }//if(mot) + }//if labm>0 + }//if(dau) + }//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(passTightCuts){ + fMassHistTC[index]->Fill(invMass); + if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass); + else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass); + } + } + } + } + } } - // end loop on D+->Kpipi - - + 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(3,fCounter); return; } + + //________________________________________________________________________ void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/) { @@ -247,17 +1078,179 @@ void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/) printf("ERROR: fOutput not available\n"); 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")); - //fNtupleDplus = dynamic_cast(fOutput->FindObject("fNtupleDplus"));//AD - fHistMass = dynamic_cast(fOutput->FindObject("fHistMass")); - fHistSignal = dynamic_cast(fOutput->FindObject("fHistSignal")); - fHistBackground = dynamic_cast(fOutput->FindObject("fHistBackground")); - //fNtupleDplusbackg = dynamic_cast(fOutput->FindObject("fNtupleDplusbackg")); //AD + TString hisname; + Int_t index=0; + - fNtupleDplus = dynamic_cast (GetOutputData(2)); //AD - fNtupleDplusbackg = dynamic_cast (GetOutputData(3)); //AD + 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); + 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())); + - return; -} + 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"); + hMassPt->SetLineColor(kBlue); + hMassPt->SetXTitle("M[GeV/c^{2}]"); + hMassPt->Draw(); + + return; +}