]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVarMultiStep.cxx
Added LHC10h run list for flow analysis (Giacomo)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFHeavyFlavourTaskMultiVarMultiStep.cxx
index ee5dbe3ec7808325851733c81f8c245bb07730ee..df07265b1295ca6bc9694e0d6a05c3472f7c81a4 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 //-----------------------------------------------------------------------
 // Class for HF corrections as a function of many variables
 // 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl, 
 #include "AliAODMCParticle.h"
 #include "AliAODMCHeader.h"
 #include "AliESDtrack.h"
+#include "AliRDHFCutsD0toKpi.h"
 #include "TChain.h"
 #include "THnSparse.h"
 #include "TH2D.h"
+#include "AliAnalysisDataSlot.h"
+#include "AliAnalysisDataContainer.h"
 
 //__________________________________________________________________________
 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep() :
@@ -69,18 +74,24 @@ AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep()
        fCountRecoAcc(0),
        fCountRecoITSClusters(0),
        fCountRecoPPR(0),
+       fCountRecoPID(0),
        fEvents(0),
        fFillFromGenerated(kFALSE),
        fMinITSClusters(5),
         fAcceptanceUnf(kTRUE),
-       fKeepD0fromB(kFALSE)
+       fKeepD0fromB(kFALSE),
+       fKeepD0fromBOnly(kFALSE),
+       fCuts(0),
+       fUseWeight(kFALSE),
+       fWeight(1.),
+       fSign(2)
 {
        //
        //Default ctor
        //
 }
 //___________________________________________________________________________
-AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name) :
+AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name, AliRDHFCutsD0toKpi* cuts) :
        AliAnalysisTaskSE(name),
        fPDG(0),
        fCFManager(0x0),
@@ -94,11 +105,17 @@ AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(c
        fCountRecoAcc(0),
        fCountRecoITSClusters(0),
        fCountRecoPPR(0),
+       fCountRecoPID(0),
        fEvents(0),
        fFillFromGenerated(kFALSE),
        fMinITSClusters(5),
         fAcceptanceUnf(kTRUE),
-       fKeepD0fromB(kFALSE)
+       fKeepD0fromB(kFALSE),
+        fKeepD0fromBOnly(kFALSE),
+       fCuts(cuts),
+       fUseWeight(kFALSE),
+       fWeight(1.),
+       fSign(2)
 {
        //
        // Constructor. Initialization of Inputs and Outputs
@@ -111,6 +128,9 @@ AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(c
        DefineOutput(1,TH1I::Class());
        DefineOutput(2,AliCFContainer::Class());
         DefineOutput(3,THnSparseD::Class());
+        DefineOutput(4,AliRDHFCutsD0toKpi::Class());
+
+       fCuts->PrintAll();
 }
 
 //___________________________________________________________________________
@@ -124,6 +144,7 @@ AliCFHeavyFlavourTaskMultiVarMultiStep& AliCFHeavyFlavourTaskMultiVarMultiStep::
                fPDG      = c.fPDG;
                fCFManager  = c.fCFManager;
                fHistEventsProcessed = c.fHistEventsProcessed;
+               fCuts = c.fCuts;
        }
        return *this;
 }
@@ -143,11 +164,17 @@ AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(c
        fCountRecoAcc(c.fCountRecoAcc),
        fCountRecoITSClusters(c.fCountRecoITSClusters),
        fCountRecoPPR(c.fCountRecoPPR),
+       fCountRecoPID(c.fCountRecoPID),
        fEvents(c.fEvents),
        fFillFromGenerated(c.fFillFromGenerated),
        fMinITSClusters(c.fMinITSClusters),
         fAcceptanceUnf(c.fAcceptanceUnf),
-       fKeepD0fromB(c.fKeepD0fromB)
+       fKeepD0fromB(c.fKeepD0fromB),
+        fKeepD0fromBOnly(c.fKeepD0fromBOnly),
+       fCuts(c.fCuts),
+       fUseWeight(c.fUseWeight),
+       fWeight(c.fWeight),
+       fSign(c.fSign)
 {
        //
        // Copy Constructor
@@ -162,6 +189,26 @@ AliCFHeavyFlavourTaskMultiVarMultiStep::~AliCFHeavyFlavourTaskMultiVarMultiStep(
        if (fCFManager)           delete fCFManager ;
        if (fHistEventsProcessed) delete fHistEventsProcessed ;
         if (fCorrelation) delete fCorrelation ;
+       //        if (fCuts) delete fCuts ;
+}
+
+//________________________________________________________________________
+void AliCFHeavyFlavourTaskMultiVarMultiStep::Init(){
+
+       //
+       // Initialization
+       //
+
+       if(fDebug > 1) printf("AliCFHeavyFlavourTaskMultiVarMultiStep::Init() \n");
+       
+       fMinITSClusters = fCuts->GetTrackCuts()->GetMinNClustersITS();
+       AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
+       const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
+       copyfCuts->SetName(nameoutput);
+       // Post the data
+       PostData(4,copyfCuts);
+       
+       return;
 }
 
 //_________________________________________________
@@ -171,6 +218,12 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
        // Main loop function
        //
        
+       PostData(1,fHistEventsProcessed) ;
+       PostData(2,fCFManager->GetParticleContainer()) ;
+        PostData(3,fCorrelation) ;
+
+       AliESDtrackCuts* trackCuts = fCuts->GetTrackCuts(); // track cuts
+
        if (fFillFromGenerated){
                AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
        }
@@ -180,8 +233,12 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                return;
        }
        
-       fEvents++;
-       if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
+       // check that the fKeepD0fromB flag is set to true when the fKeepD0fromBOnly flag is true
+       if(fKeepD0fromBOnly) { 
+         fKeepD0fromB=true;   
+         if(fEvents<2) AliInfo(Form("Both fKeepD0fromB and fKeepD0fromBOnly flags are true, looking _ONLY_ at D0 FROM B"));
+       }
+
        AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 
        TClonesArray *arrayD0toKpi=0;
@@ -209,6 +266,11 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
          return;
        }
 
+       // fix for temporary bug in ESDfilter 
+       // the AODs with null vertex pointer didn't pass the PhysSel
+       if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
+
+       fEvents++;
 
        fCFManager->SetRecEventInfo(aodEvent);
        fCFManager->SetMCEventInfo(aodEvent);
@@ -232,6 +294,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
        Int_t icountRecoAcc = 0;
        Int_t icountRecoITSClusters = 0;
        Int_t icountRecoPPR = 0;
+       Int_t icountRecoPID = 0;
        
        AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
        if (!mcHeader) {
@@ -243,6 +306,10 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                
        // AOD primary vertex
        AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
+       if(!vtx1) { 
+         AliError("There is no primary vertex !"); 
+         return; 
+       }
        Double_t zPrimVertex = vtx1->GetZ();
        Double_t zMCVertex = mcHeader->GetVtxZ();
        Bool_t vtxFlag = kTRUE;
@@ -251,12 +318,12 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
 
        for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { 
                AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
-               if (mcPart->GetPdgCode() == 4) cquarks++; 
-               if (mcPart->GetPdgCode() == -4) cquarks++; 
                if (!mcPart) {
                        AliWarning("Particle not found in tree, skipping"); 
                        continue;
                } 
+               if (mcPart->GetPdgCode() == 4) cquarks++; 
+               if (mcPart->GetPdgCode() == -4) cquarks++; 
                
                // check the MC-level cuts
 
@@ -267,11 +334,13 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                        AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
                        if (!fKeepD0fromB) continue;  // skipping particles that don't come from c quark
                }
+               else { if(fKeepD0fromBOnly) continue; } // skipping particles that don't come from b quark
                
                //              if (TMath::Abs(pdgGranma)!=4) {
 
                // fill the container for Gen-level selection
                Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
+
                if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
                        containerInputMC[0] = vectorMC[0];
                        containerInputMC[1] = vectorMC[1] ;
@@ -286,7 +355,13 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                        containerInputMC[10] = 1.01;    // dummy value, meaningless in MC
                        containerInputMC[11] = vectorMC[6];    // dummy value, meaningless in MC
                        containerInputMC[12] = zMCVertex;    // z of reconstructed of primary vertex
-                       fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated);
+                       if (fUseWeight) fWeight = GetWeight(vectorMC[0]); // setting the weight according to the function defined in AliCFHeavyFlavourTaskMultiVarMultiStep::GetWeight(Float_t pt)
+                       AliDebug(3,Form("weight = %f",fWeight));
+                       if (!fCuts->IsInFiducialAcceptance(vectorMC[0],vectorMC[1])) continue;
+                       if (TMath::Abs(vectorMC[1]) < 0.5) {
+                               fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc,fWeight);
+                       }
+                       fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated,fWeight);
                        icountMC++;
 
                        // check the MC-Acceptance level cuts
@@ -305,6 +380,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                        AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
                        if (!mcPartDaughter0 || !mcPartDaughter1) {
                                AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done..."); 
+                               continue;
                        }
                        Double_t eta0 = mcPartDaughter0->Eta();
                        Double_t eta1 = mcPartDaughter1->Eta();
@@ -325,29 +401,46 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
                                        AliDebug(2,"Inconsistency with CF cut for daughter 1!");
                                } 
-                               fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
+                               fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance,fWeight);
                                icountAcc++;
 
                                // check on the vertex
-                               if (vtxFlag){
-                                       printf("Vertex cut passed\n");
+                               if (fCuts->IsEventSelected(aodEvent)){
+                                       AliDebug(2,"Vertex cut passed\n");
                                        // filling the container if the vertex is ok
-                                       fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
-                                       printf("Vertex cut passed 2\n");
+                                       fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex,fWeight) ;
                                        icountVertex++;
                                        // check on the kTPCrefit and kITSrefit conditions of the daughters
                                        Bool_t refitFlag = kTRUE;
-                                       for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
-                                               AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
-                                               if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1)) {
-                                                 if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
-                                                   refitFlag = kFALSE;
-                                                 }
+                                       if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
+                                               Int_t foundDaughters = 0;
+                                               for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
+                                                       AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
+                                                       if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
+                                                               if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1)) {
+                                                               foundDaughters++;
+                                                               if (trackCuts->GetRequireTPCRefit()) {
+                                                                           if(!(track->GetStatus()&AliESDtrack::kTPCrefit)){
+                                                                                   refitFlag = kFALSE;
+                                                                                   break;
+                                                                           }
+                                                               }
+                                                               if (trackCuts->GetRequireITSRefit()) {
+                                                                           if(!(track->GetStatus()&AliESDtrack::kITSrefit)){
+                                                                                   refitFlag = kFALSE;
+                                                                                   break;
+                                                                           }
+                                                               }
+                                                       }
+                                                       if (foundDaughters == 2){  // both daughters have been checked
+                                                               break;
+                                                       }
                                                }
+                                               if (foundDaughters != 2) refitFlag = kFALSE;
                                        }
                                        if (refitFlag){
-                                               printf("Refit cut passed\n");
-                                               fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
+                                               AliDebug(3,"Refit cut passed\n");
+                                               fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit,fWeight);
                                                icountRefit++;
                                        }
                                        else{
@@ -382,10 +475,12 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
        AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
 
        Int_t pdgDgD0toKpi[2]={321,211};
+       Int_t isD0D0bar=1;// 1 for D0, 2 for D0bar, to be used for the PPR and PID selection steps
 
        for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
                
                AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
+               if(!d0tokpi) continue;
                Bool_t unsetvtx=kFALSE;
                if(!d0tokpi->GetOwnPrimaryVtx()) {
                  d0tokpi->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
@@ -397,6 +492,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                if (mcLabel == -1) 
                        {
                                AliDebug(2,"No MC particle found");
+                               continue;
                        }
                else {
                        AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
@@ -404,16 +500,42 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                AliWarning("Could not find associated MC in AOD MC tree");
                                continue;
                        }
+
+                       if (mcVtxHF->GetPdgCode() == 421){  // particle is D0
+                               if (fSign == 1){ // I ask for D0bar only
+                                       AliDebug(2,"particle is D0, I ask for D0bar only");
+                                       continue;
+                               }
+                               else{
+                                       isD0D0bar=1;
+                               }
+                       }
+                       else if (mcVtxHF->GetPdgCode()== -421){ // particle is D0bar
+                               if (fSign == 0){ // I ask for D0 only
+                                       AliDebug(2,"particle is D0bar, I ask for D0 only");
+                                       continue;
+                               }
+                               else{
+                                       isD0D0bar=2;
+                               }
+                       } 
+                       else continue;
+
                        // check whether the daughters have kTPCrefit and kITSrefit set
                        AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
                        AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
-                       if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
-                               // skipping if at least one daughter does not have kTPCrefit or kITSrefit
+                       if( !track0 || !track1 ) {
+                         AliWarning("Could not find associated MC daughter tracks in AOD MC tree");
+                         continue;
+                       }
+                       if ((trackCuts->GetRequireTPCRefit() && (!(track0->GetStatus()&AliESDtrack::kTPCrefit) || !(track1->GetStatus()&AliESDtrack::kTPCrefit))) || 
+                           (trackCuts->GetRequireITSRefit() && (!(track0->GetStatus()&AliESDtrack::kITSrefit) || !(track1->GetStatus()&AliESDtrack::kITSrefit)))){
+                               // skipping if at least one daughter does not have kTPCrefit or kITSrefit, if they were required
                                continue;
                        }
 
                        // check on the vertex -- could be moved outside the loop on the reconstructed D0... 
-                       if(!vtxFlag) {
+                       if(!fCuts->IsEventSelected(aodEvent)) {
                                // skipping cause vertex is not ok
                                continue;
                        }
@@ -436,6 +558,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                AliDebug(2,Form("At Reco level, from MC info: Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
                                if (!fKeepD0fromB) continue;  // skipping particles that don't come from c quark
                        }
+                       else { if(fKeepD0fromBOnly) continue; } // skipping particles that don't come from b quark
 
                        // fill the container...
                        Double_t pt = d0tokpi->Pt();
@@ -509,16 +632,23 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                        continue;
                                }
                        }
+
+                       if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue; // fiducial region
+
                        AliDebug(2, Form("Filling the container with pt = %f, rapidity = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, cT = %f", containerInput[0], containerInput[1], containerInput[2], containerInput[3], containerInput[4], containerInput[5]));
                        icountReco++;
-                       fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
+                       AliDebug(2,Form("%d: filling RECO step\n",iD0toKpi));
+                       fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed,fWeight) ;   
 
                        // cut in acceptance
-                       Bool_t acceptanceProng0 = (TMath::Abs(d0tokpi->EtaProng(0))< 0.9 && d0tokpi->PtProng(0) > 0.1);
-                       Bool_t acceptanceProng1 = (TMath::Abs(d0tokpi->EtaProng(1))< 0.9 && d0tokpi->PtProng(1) > 0.1);
+                       Float_t etaCutMin, etaCutMax, ptCutMin, ptCutMax;
+                       trackCuts->GetEtaRange(etaCutMin,etaCutMax);
+                       trackCuts->GetPtRange(ptCutMin,ptCutMax);
+                       Bool_t acceptanceProng0 = (d0tokpi->EtaProng(0)>etaCutMin && d0tokpi->EtaProng(0)<etaCutMax && d0tokpi->PtProng(0) > ptCutMin && d0tokpi->PtProng(0) < ptCutMax);
+                       Bool_t acceptanceProng1 = (d0tokpi->EtaProng(1)>etaCutMin && d0tokpi->EtaProng(1)<etaCutMax && d0tokpi->PtProng(1) > ptCutMin && d0tokpi->PtProng(1) < ptCutMax);
                        if (acceptanceProng0 && acceptanceProng1) {
                                AliDebug(2,"D0 reco daughters in acceptance");
-                               fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
+                               fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance,fWeight) ;
                                icountRecoAcc++; 
                                
                                if(fAcceptanceUnf){
@@ -540,133 +670,23 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                }  
                                
                                // cut on the min n. of clusters in ITS
-                               Int_t ncls0=0,ncls1=0;
-                               for(Int_t l=0;l<6;l++) {
-                                 if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
-                                 if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
-                               }
-                               AliDebug(2, Form("n clusters = %d", ncls0));
-                               if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters) {
+                               if (fCuts->IsSelected(d0tokpi,AliRDHFCuts::kTracks)){
                                        fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
                                        icountRecoITSClusters++;   
                                        AliDebug(2,Form("pT = %f, dca = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, d0pi = %f, d0K = %f, d0xd0 = %f, cosPointingAngle = %f", pt, dca, cosThetaStar,pTpi, pTK, d0pi*1E4, d0K*1E4, d0xd0*1E8, cosPointingAngle));
-                                       
-                                       // PPR cuts 
-                                       //added mass cut for D0 analysis
-                                       Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027}; 
-                                       //cuts of the D0 analysis (looser) see AliAnalysisTaskSED0Mass.cxx
-                                       if (pt <= 1){
-                                         
-                                         //coutputmassD02
-                                         cuts[0] = 400; //dca (um)
-                                         cuts[1] = 0.8; //cosTstar
-                                         cuts[2] = 0.5; //pt  (GeV/c)
-                                         cuts[3] = 500; //d0  (um)
-                                         cuts[4] = -25000; //d0xd0 (um^2)
-                                         cuts[5] = 0.7; //cosTpointing
-                                         /*
-                                         //PPR
-                                         cuts[0] = 400; //dca (um)
-                                         cuts[1] = 0.8; //cosTstar
-                                         cuts[2] = 0.5; //pt  (GeV/c)
-                                         cuts[3] = 500; //d0  (um)
-                                         cuts[4] = -20000; //d0xd0 (um^2)
-                                         cuts[5] = 0.5; //cosTpointing
-                                         */
-                                       }
-                                       /* //not same cuts for pt = 1 to 2 and 1 to 3 GeV in case must match them because of poor stat
-                                       else if (pt > 1 && pt <= 2){
-                                        
-                                         //coutputmassD02
-                                         cuts[0] = 300;
-                                         cuts[1] = 0.8;
-                                         cuts[2] = 0.6;
-                                         cuts[3] = 1000;
-                                         cuts[4] = -25000;
-                                         cuts[5] = 0.7;
-                                         
-                                         //PPR
-                                         cuts[0] = 300;
-                                         cuts[1] = 0.8;
-                                         cuts[2] = 0.6;
-                                         cuts[3] = 500;
-                                         cuts[4] = -20000;
-                                         cuts[5] = 0.6;
-                                         
-                                         }
-                                       */
-                                       else if (pt > 1 && pt <= 3){
-                                         
-                                         //coutputmassD02
-                                         cuts[0] = 200;
-                                         cuts[1] = 0.8;
-                                         cuts[2] = 0.7;
-                                         cuts[3] = 1000;
-                                         cuts[4] = -25000;
-                                         cuts[5] = 0.8;
-                                         /*
-                                         //PPR
-                                         cuts[0] = 300;
-                                         cuts[1] = 0.8;
-                                         cuts[2] = 0.6;
-                                         cuts[3] = 500;
-                                         cuts[4] = -20000;
-                                         cuts[5] = 0.6;
-                                         */
-                                       }
-                                       else if (pt > 3 && pt <= 5){
-                                         
-                                         //coutputmassD02
-                                         cuts[0] = 200;
-                                         cuts[1] = 0.8;
-                                         cuts[2] = 0.7;
-                                         cuts[3] = 500;
-                                         cuts[4] = -15000;
-                                         cuts[5] = 0.8;
-                                         /*
-                                         //PPR
-                                         cuts[0] = 200;
-                                         cuts[1] = 0.8;
-                                         cuts[2] = 0.7;
-                                         cuts[3] = 500;
-                                         cuts[4] = -10000;
-                                         cuts[5] = 0.8;
-                                         */
-                                       }
-                                       else if (pt > 5){
-                                         
-                                         //coutputmassD02
-                                         cuts[0] = 200;
-                                         cuts[1] = 0.8;
-                                         cuts[2] = 0.7;
-                                         cuts[3] = 500;
-                                         cuts[4] = -15000;
-                                         cuts[5] = 0.9;
-                                         /*
-                                         //PPR
-                                         cuts[0] = 200;
-                                         cuts[1] = 0.8;
-                                         cuts[2] = 0.7;
-                                         cuts[3] = 500;
-                                         cuts[4] = -5000;
-                                         cuts[5] = 0.8;
-                                         */
-                                       }
-                               
-                                       Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
-                                       if (TMath::Abs(invMass-mD0PDG) < cuts[6]
-                                           && dca*1E4 < cuts[0] 
-                                           && TMath::Abs(cosThetaStar) < cuts[1]  
-                                           && pTpi > cuts[2] 
-                                           && pTK > cuts[2]  
-                                           && TMath::Abs(d0pi*1E4) < cuts[3] 
-                                           && TMath::Abs(d0K*1E4) < cuts[3]  
-                                           && d0xd0*1E8 < cuts[4] 
-                                           && cosPointingAngle > cuts[5]
-                                           ){
-                                               
+               
+                                       // setting the use of the PID cut when applying the selection to FALSE - whatever it was. Keeping track of the original value
+                                       Bool_t iscutusingpid=fCuts->GetIsUsePID();
+                                       Int_t isselcuts=-1,isselpid=-1;
+                                       fCuts->SetUsePID(kFALSE);       
+                                       //Bool_t origFlag = fCuts->GetIsPrimaryWithoutDaughters();
+                                       //fCuts->SetRemoveDaughtersFromPrim(kFALSE);
+                                       isselcuts = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kCandidate,aodEvent);
+                                       //fCuts->SetRemoveDaughtersFromPrim(origFlag);
+                                       fCuts->SetUsePID(iscutusingpid); // restoring usage of the PID from the cuts object
+                                       if (isselcuts == 3 || isselcuts == isD0D0bar){
                                                AliDebug(2,"Particle passed PPR cuts (actually cuts for D0 analysis!)");
-                                               fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR) ;   
+                                               fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR,fWeight) ;   
                                                icountRecoPPR++;
                                                
                                                if(!fAcceptanceUnf){ // unfolding
@@ -686,32 +706,12 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                                        fCorrelation->Fill(fill);
                                                        
                                                }
-                                       }
-                                       else{
-                                               AliDebug(2,"Particle skipped due to PPR cuts");
-                                               if (dca*1E4 > cuts[0]){
-                                                       AliDebug(2,"Problems with dca");
-                                               }
-                                               if ( TMath::Abs(cosThetaStar) > cuts[1]){
-                                                       AliDebug(2,"Problems with cosThetaStar");
-                                               }
-                                               if (pTpi < cuts[2]){
-                                                       AliDebug(2,"Problems with pTpi");
-                                               }
-                                               if (pTK < cuts[2]){
-                                                       AliDebug(2,"Problems with pTK");
-                                               }
-                                               if (TMath::Abs(d0pi*1E4) > cuts[3]){
-                                                       AliDebug(2,"Problems with d0pi");
-                                               }
-                                               if (TMath::Abs(d0K*1E4) > cuts[3]){
-                                                       AliDebug(2,"Problems with d0K");
-                                               }
-                                               if (d0xd0*1E8 > cuts[4]){
-                                                       AliDebug(2,"Problems with d0xd0");
-                                               }
-                                               if (cosPointingAngle < cuts[5]){
-                                                       AliDebug(2,"Problems with cosPointingAngle");
+
+                                               isselpid = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kPID);
+                                               if((fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==isD0D0bar)||(fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==3)){
+                                                       AliDebug(2,"Particle passed PID cuts");
+                                                       fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPID,fWeight) ;   
+                                                       icountRecoPID++;
                                                }
                                        }
                                }
@@ -728,12 +728,13 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
        fCountRecoAcc+= icountRecoAcc;
        fCountRecoITSClusters+= icountRecoITSClusters;
        fCountRecoPPR+= icountRecoPPR;
+       fCountRecoPID+= icountRecoPID;
        
        fHistEventsProcessed->Fill(0);
        /* PostData(0) is taken care of by AliAnalysisTaskSE */
-       PostData(1,fHistEventsProcessed) ;
-       PostData(2,fCFManager->GetParticleContainer()) ;
-        PostData(3,fCorrelation) ;
+       //PostData(1,fHistEventsProcessed) ;
+       //PostData(2,fCFManager->GetParticleContainer()) ;
+        //PostData(3,fCorrelation) ;
 }
 
 
@@ -754,9 +755,11 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
        AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
        AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents));
        AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
+       AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PID cuts, in %d events",fCountRecoPID,fEvents));
        
        // draw some example plots....
        
+       //      AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
        AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
        if(!cont) {
          printf("CONTAINER NOT FOUND\n");
@@ -1146,7 +1149,11 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
         corr2->Draw("text");
       
 
-       TFile* file_projection = new TFile("CFtaskHFprojection.root","RECREATE");
+       TString projectionFilename="CFtaskHFprojection";
+       if(fKeepD0fromBOnly) projectionFilename+="_KeepD0fromBOnly";
+       else if(fKeepD0fromB) projectionFilename+="_KeepD0fromB";
+       projectionFilename+=".root";
+       TFile* fileProjection = new TFile(projectionFilename.Data(),"RECREATE");
 
         corr1->Write();
         corr2->Write();
@@ -1189,7 +1196,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
        h104->Write("cosPointingAngle_step2");
        h114->Write("phi_step2");
 
-       file_projection->Close();
+       fileProjection->Close();
 
     
 
@@ -1219,7 +1226,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
 }
 
 //___________________________________________________________________________
-Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
+Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(const AliAODMCParticle* mcPart, const AliAODMCParticle* mcPartDaughter0, const AliAODMCParticle* mcPartDaughter1) const {
 
        //
        // to calculate cos(ThetaStar) for generated particle
@@ -1251,7 +1258,7 @@ Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle*
        }
        if (pdgprong0 == 211){
                AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
-               AliAODMCParticle* mcPartdummy = mcPartDaughter0;
+               const AliAODMCParticle* mcPartdummy = mcPartDaughter0;
                mcPartDaughter0 = mcPartDaughter1;
                mcPartDaughter1 = mcPartdummy;
                pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
@@ -1284,7 +1291,7 @@ Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle*
        return cts;
 }
 //___________________________________________________________________________
-Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
+Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(const AliAODMCParticle* mcPart, const AliAODMCParticle* mcPartDaughter0, const AliAODMCParticle* mcPartDaughter1) const {
 
        //
        // to calculate cT for generated particle
@@ -1333,7 +1340,7 @@ Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, Al
        return cT;
 }
 //________________________________________________________________________________
-Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
+Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, const TClonesArray* mcArray, Double_t* vectorMC) const {
 
        // 
        // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
@@ -1471,7 +1478,7 @@ Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(
        return isOk;
 }
 //_________________________________________________________________________________________________
-Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPart, TClonesArray* mcArray)const{
+Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(const AliAODMCParticle* mcPart, const TClonesArray* mcArray)const{
 
        //
        // checking whether the very mother of the D0 is a charm or a bottom quark
@@ -1485,6 +1492,7 @@ Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPa
                istep++;
                AliDebug(2,Form("mother at step %d = %d", istep, mother));
                AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
+               if(!mcGranma) break;
                pdgGranma = mcGranma->GetPdgCode();
                AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
                Int_t abspdgGranma = TMath::Abs(pdgGranma);
@@ -1495,4 +1503,43 @@ Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPa
        }
        return pdgGranma;
 }
+//__________________________________________________________________________________________________
+Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetWeight(Float_t pt){
 
+       //
+       // calculating the weight to fill the container
+       //
+
+       // FNOLL central:
+       // p0 = 1.63297e-01 --> 0.322643
+       // p1 = 2.96275e+00
+       // p2 = 2.30301e+00
+       // p3 = 2.50000e+00
+
+       // PYTHIA
+       // p0 = 1.85906e-01 --> 0.36609
+       // p1 = 1.94635e+00
+       // p2 = 1.40463e+00
+       // p3 = 2.50000e+00
+
+       Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
+       Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
+
+       Double_t dndptFunc1 = DodNdptFit(pt,func1);
+       Double_t dndptFunc2 = DodNdptFit(pt,func2);
+       AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndptFunc1,dndptFunc2,dndptFunc1/dndptFunc2));
+       return dndptFunc1/dndptFunc2;
+}
+
+//__________________________________________________________________________________________________
+Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::DodNdptFit(Float_t pt, const Double_t* par)const{
+
+       // 
+       // calculating dNdpt
+       //
+
+       Double_t denom =  TMath::Power((pt/par[1]), par[3] );
+       Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
+       
+       return dNdpt;
+}