]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/vertexingHF/macros/AddTaskCFVertexingHFCascade.C
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / AddTaskCFVertexingHFCascade.C
index 433d65c5368daec7c4b404385dcad856c7746011..20fb794933beb12677ba60b546f9245937256a18 100644 (file)
@@ -45,11 +45,13 @@ const Float_t multmin_50_80 = 50;
 const Float_t multmax_50_80 = 80;
 const Float_t multmin_80_100 = 80;
 const Float_t multmax_80_100 = 100;
+const Float_t multmin_100_400 = 100;
+const Float_t multmax_100_400 = 400;
 
 
 //----------------------------------------------------
 
-AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStartoKpipiCuts010.root", TString cutObjectName="DStartoKpipiCuts", TString suffix="", Int_t configuration = AliCFTaskVertexingHF::kCheetah, Bool_t isKeepDfromB=kFALSE, Bool_t isKeepDfromBOnly=kFALSE, Int_t pdgCode = 413, Char_t isSign = 2, Bool_t useWeight=kFALSE, Bool_t useFlatPtWeight=kFALSE, Bool_t useZWeight=kFALSE, Bool_t useNchWeight=kFALSE, Bool_t isPPData=kFALSE)
+AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStartoKpipiCuts010.root", TString cutObjectName="DStartoKpipiCuts", TString suffix="suf", Int_t configuration = AliCFTaskVertexingHF::kCheetah, Bool_t isKeepDfromB=kFALSE, Bool_t isKeepDfromBOnly=kFALSE, Int_t pdgCode = 413, Char_t isSign = 2, Bool_t useWeight=kTRUE, Bool_t useFlatPtWeight=kFALSE, Bool_t useZWeight=kFALSE, Bool_t useNchWeight=kFALSE, Bool_t useNtrkWeight=kFALSE, TString estimatorFilename="", Int_t multiplicityEstimator = AliCFTaskVertexingHF::kNtrk10, Bool_t isPPData=kFALSE, Bool_t isPPbData=kFALSE, Double_t refMult=9.26, Bool_t isFineNtrkBin=kFALSE)
 {
        printf("Adding CF task using cuts from file %s\n",cutFile);
        if (configuration == AliCFTaskVertexingHF::kSnail){
@@ -89,7 +91,7 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
          return 0x0;
        }
 
-       AliRDHFCutsD0toKpi *cutsD0toKpi = (AliRDHFCutsD0toKpi*)fileCuts->Get(cutObjectName.Data());
+       AliRDHFCutsDStartoKpipi *cutsDStartoKpi = (AliRDHFCutsDStartoKpipi*)fileCuts->Get(cutObjectName.Data());
        
        // check that the fKeepD0fromB flag is set to true when the fKeepD0fromBOnly flag is true
        //  for now the binning is the same than for all D's
@@ -138,8 +140,10 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
        const Int_t nbinmult = 49;  //bins in multiplicity (total number)
         const Int_t nbinmult_0_20 = 20; //bins in multiplicity between 0 and 20
         const Int_t nbinmult_20_50 = 15; //bins in multiplicity between 20 and 50
-        const Int_t nbinmult_50_80 = 10; //bins in multiplicity between 50 and 102
-        const Int_t nbinmult_80_100 = 4; //bins in multiplicity between 50 and 102
+        const Int_t nbinmult_50_80 = 10; //bins in multiplicity between 50 and 80
+        const Int_t nbinmult_80_100 = 4; //bins in multiplicity between 80 and 100
+        const Int_t nbinmult_100_400 = 3; //bins in multiplicity between 100 and 400
+       if(isPPbData) nbinmult += nbinmult_100_400;
 
 
        //the sensitive variables, their indices
@@ -169,7 +173,7 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
        //arrays for the number of bins in each dimension
        Int_t iBin[nvarTot];
                
-       const Int_t nbinpt = cutsD0toKpi->GetNPtBins(); // bins in pT
+       const Int_t nbinpt = cutsDStartoKpi->GetNPtBins(); // bins in pT
        iBin[ipT]=nbinpt;
        iBin[ipTpi]=nbinpt;
        iBin[ipTk]=nbinpt;
@@ -177,7 +181,7 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
        Double_t *binLimpTpi=new Double_t[iBin[ipTpi]+1];
        Double_t *binLimpTk=new Double_t[iBin[ipTk]+1];
        // values for bin lower bounds
-       Float_t* floatbinLimpT = cutsD0toKpi->GetPtBinLimits();
+       Float_t* floatbinLimpT = cutsDStartoKpi->GetPtBinLimits();
        for (Int_t ibin0 = 0 ; ibin0<iBin[ipT]+1; ibin0++){
                binLimpT[ibin0] = (Double_t)floatbinLimpT[ibin0];
                binLimpTpi[ibin0] = (Double_t)floatbinLimpT[ibin0];
@@ -186,7 +190,29 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
        for(Int_t i=0; i<=nbinpt; i++) printf("binLimpT[%d]=%f\n",i,binLimpT[i]);  
        
        printf("pT: nbin (from cuts file) = %d\n",nbinpt);
+
+
        
+
+
+
+
+
+       // Fine Ntrk binning setting
+       Double_t *binLimmultFine;
+       Int_t nbinmultTmp=nbinmult;
+       if (isFineNtrkBin) {
+               Int_t nbinLimmultFine = 100;
+               if (isPPbData) nbinLimmultFine = 200;
+               const UInt_t nbinMultFine = nbinLimmultFine;
+               binLimmultFine = new Double_t[nbinMultFine+1];
+               for (Int_t ibin0 = 0; ibin0 < nbinMultFine+1; ibin0++) {
+                       binLimmultFine[ibin0] = ibin0;
+               }
+               nbinmultTmp = nbinLimmultFine;
+       }
+       const Int_t nbinmultTot = nbinmultTmp;
+
        // defining now the binning for the other variables:
        
        AliLog::SetClassDebugLevel("AliCFManager",AliLog::kInfo);
@@ -203,7 +229,7 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
        iBin[ifake]=nbinfake;
        iBin[ipointingXY]=nbinpointingXY;
        iBin[inormDecayLXY]=nbinnormDecayLXY;
-       iBin[imult]=nbinmult;
+       iBin[imult]=nbinmultTot;
        
        //arrays for lower bounds :
        Double_t *binLimy=new Double_t[iBin[iy]+1];
@@ -281,9 +307,15 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
         }
         for(Int_t i=0; i<=nbinmult_50_80; i++) binLimmult[i+nbinmult_0_20+nbinmult_20_50]=(Double_t)multmin_50_80 + (multmax_50_80-multmin_50_80)/nbinmult_50_80*(Double_t)i ; 
         if (binLimmult[nbinmult_0_20+nbinmult_20_50+nbinmult_50_80] != multmin_80_100)  {
-                Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for mult - 2nd range - differs from expected!\n");
+                Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for mult - 3rd range - differs from expected!\n");
         }
         for(Int_t i=0; i<=nbinmult_80_100; i++) binLimmult[i+nbinmult_0_20+nbinmult_20_50+nbinmult_50_80]=(Double_t)multmin_80_100 + (multmax_80_100-multmin_80_100)/nbinmult_80_100*(Double_t)i ; 
+        if (binLimmult[nbinmult_0_20+nbinmult_20_50+nbinmult_50_80+nbinmult_80_100] != multmin_100_400) {
+                Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for mult - 4th range - differs from expected!\n");
+        }
+       if(isPPbData){
+         for (Int_t i = 0; i<=nbinmult_100_400; i++) binLimmult[i+nbinmult_0_20+nbinmult_20_50+nbinmult_50_80+nbinmult_80_100]= (Double_t)multmin_100_400 + (multmax_100_400-multmin_100_400)/nbinmult_100_400*(Double_t)i ;
+       }
 
        //one "container" for MC
        TString nameContainer="";
@@ -335,7 +367,9 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
                printf("normDecayLXY\n");
                container -> SetBinLimits(inormDecayLXY,binLimnormDecayLXY);
                printf("multiplicity\n");
-               container -> SetBinLimits(imult,binLimmult);
+               
+               if (isFineNtrkBin) container->SetBinLimits(imult,binLimmultFine);
+               else               container -> SetBinLimits(imult,binLimmult);
 
                container -> SetVarTitle(ipT,"pt");
                container -> SetVarTitle(iy,"y");
@@ -393,7 +427,8 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
                printf("fake\n");
                container -> SetBinLimits(ifakeFast,binLimfake);
                printf("multiplicity\n");
-               container -> SetBinLimits(imultFast,binLimmult);
+               if(isFineNtrkBin) container -> SetBinLimits(imultFast,binLimmultFine);
+               else              container -> SetBinLimits(imultFast,binLimmult);
 
                container -> SetVarTitle(ipTFast,"pt");
                container -> SetVarTitle(iyFast,"y");
@@ -488,7 +523,7 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
        printf("CREATE TASK\n");
 
        // create the task
-       AliCFTaskVertexingHF *task = new AliCFTaskVertexingHF("AliCFTaskVertexingHF",cutsD0toKpi);
+       AliCFTaskVertexingHF *task = new AliCFTaskVertexingHF("AliCFTaskVertexingHF",cutsDStartoKpi);
        task->SetConfiguration(configuration);
        task->SetFillFromGenerated(kFALSE);
        task->SetCFManager(man); //here is set the CF manager
@@ -516,40 +551,90 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
                else {
                        task->SetWeightFunction(funcWeight);
                        Printf("User-defined Weights will be used. The function being:");
-                       task->GetWeightFunction(funcWeight)->Print();
+                       task->GetWeightFunction()->Print();
                }
        }
 
-       if(useNchWeight){
-          TH1F *hNchPrimaries = (TH1F*)fileCuts->Get("hGenPrimaryParticlesInelGt0");
+       if(useNchWeight || useNtrkWeight){
+          TH1F *hNchPrimaries;
+         TH1F *hNchMeasured;           
+         if (isPPbData) hNchPrimaries = (TH1F*)fileCuts->Get("hNtrUnCorrEvWithCandWeight");
+         else hNchPrimaries = (TH1F*)fileCuts->Get("hGenPrimaryParticlesInelGt0");
+         hNchMeasured = (TH1F*)fileCuts->Get("hNchMeasured");
           if(hNchPrimaries) {
             task->SetUseNchWeight(kTRUE);
             task->SetMCNchHisto(hNchPrimaries);
+           if(isPPbData) task->SetUseNchTrackletsWeight();
           } else {
             AliFatal("Histogram for multiplicity weights not found");
             return 0x0;
           }
+         if(hNchMeasured) task->SetMeasuredNchHisto(hNchMeasured);
+         if(useNtrkWeight) task->SetUseNchTrackletsWeight();
         }
-       
-       Printf("***************** CONTAINER SETTINGS *****************");       
-       Printf("decay channel = %d",(Int_t)task->GetDecayChannel());
-       Printf("FillFromGenerated = %d",(Int_t)task->GetFillFromGenerated());
-       Printf("Dselection = %d",(Int_t)task->GetDselection());
-       Printf("UseWeight = %d",(Int_t)task->GetUseWeight());
-       if (task->GetUseWeight()) {
-               Printf("User-defined Weight function:");
-               task->GetWeightFunction(funcWeight)->Print();
-       }
-       else{
-               Printf("FONLL will be used for the weights");
+
+       task->SetMultiplicityEstimator(multiplicityEstimator);
+       if(estimatorFilename.EqualTo("") ) {
+         printf("Estimator file not provided, multiplicity corrected histograms will not be filled\n");
+         task->SetUseZvtxCorrectedNtrkEstimator(kFALSE);
+       } else{
+         TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
+         if(!fileEstimator)  {
+           AliFatal("File with multiplicity estimator not found"); 
+           return;
+         }      
+         task->SetUseZvtxCorrectedNtrkEstimator(kTRUE);
+         task->SetReferenceMultiplcity(refMult);
+         
+         if (isPPbData) {  //load multiplicity estimators for pPb
+           task->SetIsPPbData(kTRUE);
+           const Char_t* periodNames[2] = {"LHC13b", "LHC13c"};
+           TProfile *multEstimatorAvg[2];
+           for (Int_t ip=0; ip < 2; ip++) {
+             multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));    
+             if (!multEstimatorAvg[ip]) {
+               AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
+               return;               
+             }
+           }
+           task->SetMultiplVsZProfileLHC13b(multEstimatorAvg[0]);
+           task->SetMultiplVsZProfileLHC13c(multEstimatorAvg[1]);
+           
+         } else {    //load multiplicity estimators for pp
+           const Char_t* periodNames[4] = {"LHC10b", "LHC10c", "LHC10d", "LHC10e"};
+           TProfile* multEstimatorAvg[4];                       
+           
+           for(Int_t ip=0; ip<4; ip++) {
+             multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));  
+             if (!multEstimatorAvg[ip]) {
+               AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
+               return;               
+             }  
+           }
+           task->SetMultiplVsZProfileLHC10b(multEstimatorAvg[0]);
+           task->SetMultiplVsZProfileLHC10c(multEstimatorAvg[1]);
+           task->SetMultiplVsZProfileLHC10d(multEstimatorAvg[2]);
+           task->SetMultiplVsZProfileLHC10e(multEstimatorAvg[3]);
+         }
        }
-       Printf("Use Nch weight = %d",(Int_t)task->GetUseNchWeight());
-       Printf("Sign = %d",(Int_t)task->GetSign());
-       Printf("Centrality selection = %d",(Int_t)task->GetCentralitySelection());
-       Printf("Fake selection = %d",(Int_t)task->GetFakeSelection());
-       Printf("RejectCandidateIfNotFromQuark selection = %d",(Int_t)task->GetRejectCandidateIfNotFromQuark());
-       Printf("UseMCVertex selection = %d",(Int_t)task->GetUseMCVertex());
-       Printf("***************END CONTAINER SETTINGS *****************\n");
+       Printf("***************** CONTAINER SETTINGS *****************");       
+        Printf("decay channel = %d",(Int_t)task->GetDecayChannel());
+        Printf("FillFromGenerated = %d",(Int_t)task->GetFillFromGenerated());
+        Printf("Dselection = %d",(Int_t)task->GetDselection());
+        Printf("UseWeight = %d",(Int_t)task->GetUseWeight());
+        if (task->GetUseWeight()) {
+          if(funcWeight) Printf("User-defined Weight function");
+          else Printf("FONLL will be used for the weights");
+        }
+
+        Printf("Use Nch weight = %d",(Int_t)task->GetUseNchWeight());
+        Printf("Sign = %d",(Int_t)task->GetSign());
+        Printf("Centrality selection = %d",(Int_t)task->GetCentralitySelection());
+        Printf("Fake selection = %d",(Int_t)task->GetFakeSelection());
+        Printf("RejectCandidateIfNotFromQuark selection = %d",(Int_t)task->GetRejectCandidateIfNotFromQuark());
+        Printf("UseMCVertex selection = %d",(Int_t)task->GetUseMCVertex());
+        Printf("***************END CONTAINER SETTINGS *****************\n");
+
 
         //-----------------------------------------------------------//
         //   create correlation matrix for unfolding - only eta-pt   //
@@ -609,32 +694,37 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
        // ----- output data -----
        
        TString outputfile = AliAnalysisManager::GetCommonFileName();
-       TString output1name="", output2name="", output3name="",output4name="";
+       TString output1name="", output2name="", output3name="",output4name="", output5name="";
        output2name=nameContainer;
        output3name=nameCorr;
        output4name= "Cuts";
+       output5name= "coutProfDst";
        if(!isKeepDfromB) {
                outputfile += ":PWG3_D2H_CFtaskDstartoKpipi";
                output1name="CFHFchist0";
                output3name+="_cOnly";
                output4name+="_cOnly";
+               output5name+="_cOnly";
        }
        else  if(isKeepDfromBOnly){
                outputfile += ":PWG3_D2H_CFtaskDstartoKpiKeepDfromBOnly";
                output1name="CFHFchist0DfromB";
                output3name+="_bOnly";
                output4name+="_bOnly";
+               output5name+="_bOnly";
        }
        else{
                outputfile += ":PWG3_D2H_CFtaskDstartoKpiKeepDfromB";
                output1name="CFHFchist0allD";
                output3name+="_all";
                output4name+="_all";
+               output5name+="_all";
        }
 
        outputfile += suffix;
        output1name += suffix;
        output4name += suffix;
+       output5name += suffix;
 
        //now comes user's output objects :
        // output TH1I for event counting
@@ -645,6 +735,8 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
         AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(output3name, THnSparseD::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
        // cuts
        AliAnalysisDataContainer *coutput4 = mgr->CreateContainer(output4name, AliRDHFCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());
+       // estimators list
+       AliAnalysisDataContainer *coutput5 = mgr->CreateContainer(output5name, TList::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());
 
        mgr->AddTask(task);
        
@@ -653,8 +745,7 @@ AliCFTaskVertexingHF *AddTaskCFVertexingHFCascade(const char* cutFile = "DStarto
        mgr->ConnectOutput(task,2,coutput2);
         mgr->ConnectOutput(task,3,coutput3);
        mgr->ConnectOutput(task,4,coutput4);
+       mgr->ConnectOutput(task,5,coutput5);
        return task;
        
 }
-
-