PWGHF/vertexingHF/macros/AddTaskDvsMultiplicity.C

   1 AliAnalysisTaskSEDvsMultiplicity *AddTaskDvsMultiplicity(Int_t system=0,
   2                                                          Bool_t readMC=kFALSE,
   3                                                          Int_t MCOption=0,
   4                                                          Int_t pdgMeson=411,
   5                                                          TString finDirname="Loose",
   6                                                          TString filename="",
   7                                                          TString finAnObjname="AnalysisCuts",
   8                                                          TString estimatorFilename="",
   9                                                          Double_t refMult=9.26,
  10                                                          Bool_t subtractDau=kFALSE,
  11                                                          Bool_t NchWeight=kFALSE,
  12                                                          Int_t recoEstimator = AliAnalysisTaskSEDvsMultiplicity::kNtrk10,
  13                                                          Int_t MCEstimator = AliAnalysisTaskSEDvsMultiplicity::kEta10,
  14                                                          Bool_t isPPbData=kFALSE)
  15 {
  16   //
  17   // Test macro for the AliAnalysisTaskSE for D+ candidates
  18   //Invariant mass histogram and
  19   // association with MC truth (using MC info in AOD)
  20   //  R. Bala, bala@to.infn.it
  21   // Get the pointer to the existing analysis manager via the static access method.
  22   //==============================================================================
  23
  24   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  25   if (!mgr) {
  26     ::Error("AddTaskDvsMultiplicity", "No analysis manager to connect to.");
  27   }
  28
  29   Bool_t stdcuts=kFALSE;
  30   TFile* filecuts;
  31   if( filename.EqualTo("") ) {
  32     stdcuts=kTRUE;
  33   } else {
  34     filecuts=TFile::Open(filename.Data());
  35     if(!filecuts ||(filecuts&& !filecuts->IsOpen())){
  36       AliFatal("Input file not found : check your cut object");
  37     }
  38   }
  39
  40
  41   //Analysis Task
  42   AliRDHFCuts *analysiscuts=0x0;
  43
  44   TString Name="";
  45   if(pdgMeson==411){
  46     if(stdcuts) {
  47       analysiscuts = new AliRDHFCutsDplustoKpipi();
  48       if (system == 0) analysiscuts->SetStandardCutsPP2010();
  49       else analysiscuts->SetStandardCutsPbPb2011();
  50     }
  51     else analysiscuts = (AliRDHFCutsDplustoKpipi*)filecuts->Get(finAnObjname);
  52     Name="Dplus";
  53    }else if(pdgMeson==421){
  54     if(stdcuts) {
  55       analysiscuts = new AliRDHFCutsD0toKpi();
  56       if (system == 0) analysiscuts->SetStandardCutsPP2010();
  57       else analysiscuts->SetStandardCutsPbPb2011();
  58     }
  59     else analysiscuts = (AliRDHFCutsD0toKpi*)filecuts->Get(finAnObjname);
  60     Name="D0";
  61   }else if(pdgMeson==413){
  62     if(stdcuts) {
  63       analysiscuts = new AliRDHFCutsDStartoKpipi();
  64       if (system == 0) analysiscuts->SetStandardCutsPP2010();
  65       else analysiscuts->SetStandardCutsPbPb2011();
  66     }
  67     else analysiscuts = (AliRDHFCutsDStartoKpipi*)filecuts->Get(finAnObjname);
  68     Name="DStar";
  69   }
  70
  71   AliAnalysisTaskSEDvsMultiplicity *dMultTask = new AliAnalysisTaskSEDvsMultiplicity("dMultAnalysis",pdgMeson,analysiscuts,isPPbData);
  72   dMultTask->SetReadMC(readMC);
  73   dMultTask->SetDebugLevel(0);
  74   dMultTask->SetUseBit(kTRUE);
  75   dMultTask->SetDoImpactParameterHistos(kFALSE);
  76   dMultTask->SetSubtractTrackletsFromDaughters(subtractDau);
  77   dMultTask->SetMultiplicityEstimator(recoEstimator);
  78   dMultTask->SetMCPrimariesEstimator(MCEstimator);
  79   dMultTask->SetMCOption(MCOption);
  80   if(isPPbData) dMultTask->SetIsPPbData();
  81
  82   if(NchWeight){
  83     TH1F *hNchPrimaries = NULL;
  84     if(isPPbData) hNchPrimaries = (TH1F*)filecuts->Get("hNtrUnCorrEvWithDWeight");
  85     else hNchPrimaries = (TH1F*)filecuts->Get("hGenPrimaryParticlesInelGt0");
  86     if(hNchPrimaries) {
  87       dMultTask->UseMCNchWeight(true);
  88       dMultTask->SetHistoNchWeight(hNchPrimaries);
  89     } else {
  90       AliFatal("Histogram for multiplicity weights not found");
  91       return 0x0;
  92     }
  93   }
  94
  95   if(pdgMeson==421) {
  96     dMultTask->SetMassLimits(1.5648,2.1648);
  97     dMultTask->SetNMassBins(200);
  98   }else if(pdgMeson==411)dMultTask->SetMassLimits(pdgMeson,0.2);
  99
 100   if(estimatorFilename.EqualTo("") ) {
 101     printf("Estimator file not provided, multiplcity corrected histograms will not be filled\n");
 102   } else{
 103
 104     TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
 105     if(!fileEstimator)  {
 106       AliFatal("File with multiplicity estimator not found\n");
 107       return;
 108     }
 109
 110     dMultTask->SetReferenceMultiplcity(refMult);
 111
 112     if (isPPbData) {    //Only use two profiles if pPb
 113       const Char_t* periodNames[2] = {"LHC13b", "LHC13c"};
 114       TProfile* multEstimatorAvg[2];
 115       for(Int_t ip=0; ip<2; ip++) {
 116         multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));
 117         if (!multEstimatorAvg[ip]) {
 118           AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
 119           return;
 120         }
 121       }
 122       dMultTask->SetMultiplVsZProfileLHC13b(multEstimatorAvg[0]);
 123       dMultTask->SetMultiplVsZProfileLHC13c(multEstimatorAvg[1]);
 124     }
 125     else {
 126       const Char_t* periodNames[4] = {"LHC10b", "LHC10c", "LHC10d", "LHC10e"};
 127       TProfile* multEstimatorAvg[4];
 128       for(Int_t ip=0; ip<4; ip++) {
 129         multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));
 130         if (!multEstimatorAvg[ip]) {
 131           AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
 132           return;
 133         }
 134       }
 135       dMultTask->SetMultiplVsZProfileLHC10b(multEstimatorAvg[0]);
 136       dMultTask->SetMultiplVsZProfileLHC10c(multEstimatorAvg[1]);
 137       dMultTask->SetMultiplVsZProfileLHC10d(multEstimatorAvg[2]);
 138       dMultTask->SetMultiplVsZProfileLHC10e(multEstimatorAvg[3]);
 139     }
 140   }
 141   mgr->AddTask(dMultTask);
 142
 143   // Create containers for input/output
 144
 145   TString inname = "cinput";
 146   TString outname = "coutput";
 147   TString cutsname = "coutputCuts";
 148   TString normname = "coutputNorm";
 149   TString profname = "coutputProf";
 150
 151   inname += Name.Data();
 152   outname += Name.Data();
 153   cutsname += Name.Data();
 154   normname += Name.Data();
 155   profname += Name.Data();
 156   inname += finDirname.Data();
 157   outname += finDirname.Data();
 158   cutsname += finDirname.Data();
 159   normname += finDirname.Data();
 160   profname += finDirname.Data();
 161
 162   AliAnalysisDataContainer *cinput = mgr->CreateContainer(inname,TChain::Class(),AliAnalysisManager::kInputContainer);
 163
 164   TString outputfile = AliAnalysisManager::GetCommonFileName();
 165   outputfile += ":PWG3_D2H_DMult_";
 166   outputfile += Name.Data();
 167   outputfile += finDirname.Data();
 168
 169   AliAnalysisDataContainer *coutputCuts = mgr->CreateContainer(cutsname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
 170   AliAnalysisDataContainer *coutput = mgr->CreateContainer(outname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
 171   AliAnalysisDataContainer *coutputNorm = mgr->CreateContainer(normname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
 172   AliAnalysisDataContainer *coutputProf = mgr->CreateContainer(profname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
 173
 174   mgr->ConnectInput(dMultTask,0,mgr->GetCommonInputContainer());
 175
 176   mgr->ConnectOutput(dMultTask,1,coutput);
 177
 178   mgr->ConnectOutput(dMultTask,2,coutputCuts);
 179
 180   mgr->ConnectOutput(dMultTask,3,coutputNorm);
 181
 182   mgr->ConnectOutput(dMultTask,4,coutputProf);
 183
 184   return dMultTask;
 185 }