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