[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / AddTaskDvsMultiplicity.C

AliAnalysisTaskSEDvsMultiplicity *AddTaskDvsMultiplicity(Int_t system=0,
							 Bool_t readMC=kFALSE,
							 Int_t MCOption=0,
							 Int_t pdgMeson=411,
							 TString finDirname="Loose",
							 TString filename="",
							 TString finAnObjname="AnalysisCuts", 
							 TString estimatorFilename="",
							 Double_t refMult=9.26,
							 Bool_t subtractDau=kFALSE,
							 Bool_t NchWeight=kFALSE,
							 Int_t recoEstimator = AliAnalysisTaskSEDvsMultiplicity::kNtrk10,
							 Int_t MCEstimator = AliAnalysisTaskSEDvsMultiplicity::kEta10,
							 Bool_t isPPbData=kFALSE)
{
  //
  // Test macro for the AliAnalysisTaskSE for D+ candidates
  //Invariant mass histogram and    
  // association with MC truth (using MC info in AOD)
  //  R. Bala, bala@to.infn.it
  // Get the pointer to the existing analysis manager via the static access method. 
  //============================================================================== 
  
  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  if (!mgr) {
    ::Error("AddTaskDvsMultiplicity", "No analysis manager to connect to.");
  }
  
  Bool_t stdcuts=kFALSE;
  TFile* filecuts;
  if( filename.EqualTo("") ) {
    stdcuts=kTRUE; 
  } else {
    filecuts=TFile::Open(filename.Data());
    if(!filecuts ||(filecuts&& !filecuts->IsOpen())){
      AliFatal("Input file not found : check your cut object");
    }
  }

  
  //Analysis Task  
  AliRDHFCuts *analysiscuts=0x0;
  
  TString Name="";
  if(pdgMeson==411){
    if(stdcuts) {
      analysiscuts = new AliRDHFCutsDplustoKpipi();
      if (system == 0) analysiscuts->SetStandardCutsPP2010();
      else analysiscuts->SetStandardCutsPbPb2011();
    }
    else analysiscuts = (AliRDHFCutsDplustoKpipi*)filecuts->Get(finAnObjname);
    Name="Dplus";
   }else if(pdgMeson==421){
    if(stdcuts) {
      analysiscuts = new AliRDHFCutsD0toKpi();
      if (system == 0) analysiscuts->SetStandardCutsPP2010();
      else analysiscuts->SetStandardCutsPbPb2011();
    }
    else analysiscuts = (AliRDHFCutsD0toKpi*)filecuts->Get(finAnObjname);
    Name="D0";
  }else if(pdgMeson==413){
    if(stdcuts) {
      analysiscuts = new AliRDHFCutsDStartoKpipi();
      if (system == 0) analysiscuts->SetStandardCutsPP2010();
      else analysiscuts->SetStandardCutsPbPb2011();
    }
    else analysiscuts = (AliRDHFCutsDStartoKpipi*)filecuts->Get(finAnObjname);
    Name="DStar";
  }

  AliAnalysisTaskSEDvsMultiplicity *dMultTask = new AliAnalysisTaskSEDvsMultiplicity("dMultAnalysis",pdgMeson,analysiscuts,isPPbData);
  dMultTask->SetReadMC(readMC);  
  dMultTask->SetDebugLevel(0);
  dMultTask->SetUseBit(kTRUE);
  dMultTask->SetDoImpactParameterHistos(kFALSE);
  dMultTask->SetSubtractTrackletsFromDaughters(subtractDau);
  dMultTask->SetMultiplicityEstimator(recoEstimator);
  dMultTask->SetMCPrimariesEstimator(MCEstimator);
  dMultTask->SetMCOption(MCOption);
  if(isPPbData) dMultTask->SetIsPPbData();

  if(NchWeight){
    TH1F *hNchPrimaries = (TH1F*)filecuts->Get("hGenPrimaryParticlesInelGt0");
    if(hNchPrimaries) {
      dMultTask->UseMCNchWeight(true);
      dMultTask->SetHistoNchWeight(hNchPrimaries);
    } else {
      AliFatal("Histogram for multiplicity weights not found");
      return 0x0;
    }
  }

  if(pdgMeson==421) { 
    dMultTask->SetMassLimits(1.5648,2.1648);
    dMultTask->SetNMassBins(200);
  }else if(pdgMeson==411)dMultTask->SetMassLimits(pdgMeson,0.2);
  
  if(estimatorFilename.EqualTo("") ) {
    printf("Estimator file not provided, multiplcity corrected histograms will not be filled\n");
  } else{
         if ( isPPbData) {    //Only use two profiles if pPb
            const Char_t* periodNames[2] = {"LHC13b", "LHC13c"};
            TProfile* multEstimatorAvg[2];                       
            TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
            if(!fileEstimator)  {
               AliFatal("File with multiplicity estimator not found\n"); 
               return;
            }
            
            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])));  
            }
            dMultTask->SetMultiplVsZProfileLHC13b(multEstimatorAvg[0]);
            dMultTask->SetMultiplVsZProfileLHC13c(multEstimatorAvg[1]);
            dMultTask->SetReferenceMultiplcity(refMult);
            
         }
         else {
            const Char_t* periodNames[4] = {"LHC10b", "LHC10c", "LHC10d", "LHC10e"};
            TProfile* multEstimatorAvg[4];                       
            TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
            if(!fileEstimator)  {
               AliFatal("File with multiplicity estimator not found\n"); 
               return;
            }
            
            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])));  
            }
            dMultTask->SetMultiplVsZProfileLHC10b(multEstimatorAvg[0]);
            dMultTask->SetMultiplVsZProfileLHC10c(multEstimatorAvg[1]);
            dMultTask->SetMultiplVsZProfileLHC10d(multEstimatorAvg[2]);
            dMultTask->SetMultiplVsZProfileLHC10e(multEstimatorAvg[3]);
            dMultTask->SetReferenceMultiplcity(refMult);
   }
  }
  mgr->AddTask(dMultTask);
  
  // Create containers for input/output 
  
  TString inname = "cinput";
  TString outname = "coutput";
  TString cutsname = "coutputCuts";
  TString normname = "coutputNorm";
  TString profname = "coutputProf";
  
  inname += Name.Data();
  outname += Name.Data();
  cutsname += Name.Data();
  normname += Name.Data();
  profname += Name.Data();
  inname += finDirname.Data();
  outname += finDirname.Data();
  cutsname += finDirname.Data();
  normname += finDirname.Data();
  profname += finDirname.Data();

  AliAnalysisDataContainer *cinput = mgr->CreateContainer(inname,TChain::Class(),AliAnalysisManager::kInputContainer);

  TString outputfile = AliAnalysisManager::GetCommonFileName();
  outputfile += ":PWG3_D2H_DMult_";
  outputfile += Name.Data(); 
  outputfile += finDirname.Data(); 
    
  AliAnalysisDataContainer *coutputCuts = mgr->CreateContainer(cutsname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
  AliAnalysisDataContainer *coutput = mgr->CreateContainer(outname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
  AliAnalysisDataContainer *coutputNorm = mgr->CreateContainer(normname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
  AliAnalysisDataContainer *coutputProf = mgr->CreateContainer(profname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
  
  mgr->ConnectInput(dMultTask,0,mgr->GetCommonInputContainer());
  
  mgr->ConnectOutput(dMultTask,1,coutput);
  
  mgr->ConnectOutput(dMultTask,2,coutputCuts);

  mgr->ConnectOutput(dMultTask,3,coutputNorm);  

  mgr->ConnectOutput(dMultTask,4,coutputProf);

  return dMultTask;
}
Commit	Line	Data
b7fae104	1	AliAnalysisTaskSEDvsMultiplicity *AddTaskDvsMultiplicity(Int_t system=0,
	2	Bool_t readMC=kFALSE,
	3	Int_t MCOption=0,
3d27f7b3	4	Int_t pdgMeson=411,
b7fae104	5	TString finDirname="Loose",
	6	TString filename="",
	7	TString finAnObjname="AnalysisCuts",
38d3a7e8	8	TString estimatorFilename="",
0d1e32a6	9	Double_t refMult=9.26,
18d726fa	10	Bool_t subtractDau=kFALSE,
18d726fa	11	Bool_t NchWeight=kFALSE,
20550d50	12	Int_t recoEstimator = AliAnalysisTaskSEDvsMultiplicity::kNtrk10,
bba1a026	13	Int_t MCEstimator = AliAnalysisTaskSEDvsMultiplicity::kEta10,
bba1a026	14	Bool_t isPPbData=kFALSE)
b7fae104	15	{
3d27f7b3	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
b7fae104	24	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	25	if (!mgr) {
	26	::Error("AddTaskDvsMultiplicity", "No analysis manager to connect to.");
	27	}
3d27f7b3	28
b7fae104	29	Bool_t stdcuts=kFALSE;
2d11aff8	30	TFile* filecuts;
b7fae104	31	if( filename.EqualTo("") ) {
	32	stdcuts=kTRUE;
	33	} else {
3d27f7b3	34	filecuts=TFile::Open(filename.Data());
	35	if(!filecuts \|\|(filecuts&& !filecuts->IsOpen())){
	36	AliFatal("Input file not found : check your cut object");
	37	}
b7fae104	38	}
b7fae104	39
b7fae104	40
3d27f7b3	41	//Analysis Task
3d27f7b3	42	AliRDHFCuts *analysiscuts=0x0;
b7fae104	43
3d27f7b3	44	TString Name="";
3d27f7b3	45	if(pdgMeson==411){
b7fae104	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";
3d27f7b3	53	}else if(pdgMeson==421){
b7fae104	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);
3d27f7b3	60	Name="D0";
	61	}else if(pdgMeson==413){
	62	if(stdcuts) {
b7fae104	63	analysiscuts = new AliRDHFCutsDStartoKpipi();
	64	if (system == 0) analysiscuts->SetStandardCutsPP2010();
	65	else analysiscuts->SetStandardCutsPbPb2011();
	66	}
	67	else analysiscuts = (AliRDHFCutsDStartoKpipi*)filecuts->Get(finAnObjname);
3d27f7b3	68	Name="DStar";
3d27f7b3	69	}
b7fae104	70
b97c5dc0	71	AliAnalysisTaskSEDvsMultiplicity *dMultTask = new AliAnalysisTaskSEDvsMultiplicity("dMultAnalysis",pdgMeson,analysiscuts,isPPbData);
3d27f7b3	72	dMultTask->SetReadMC(readMC);
b7fae104	73	dMultTask->SetDebugLevel(0);
b7fae104	74	dMultTask->SetUseBit(kTRUE);
b7fae104	75	dMultTask->SetDoImpactParameterHistos(kFALSE);
0d1e32a6	76	dMultTask->SetSubtractTrackletsFromDaughters(subtractDau);
20550d50	77	dMultTask->SetMultiplicityEstimator(recoEstimator);
	78	dMultTask->SetMCPrimariesEstimator(MCEstimator);
	79	dMultTask->SetMCOption(MCOption);
bba1a026	80	if(isPPbData) dMultTask->SetIsPPbData();
18d726fa	81
	82	if(NchWeight){
	83	TH1F 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	}
14eaa841	92
	93	if(pdgMeson==421) {
	94	dMultTask->SetMassLimits(1.5648,2.1648);
8934c74b	95	dMultTask->SetNMassBins(200);
c3f04bec	96	}else if(pdgMeson==411)dMultTask->SetMassLimits(pdgMeson,0.2);
3d27f7b3	97
2d11aff8	98	if(estimatorFilename.EqualTo("") ) {
f09fd3fe	99	printf("Estimator file not provided, multiplcity corrected histograms will not be filled\n");
3d27f7b3	100	} else{
b97c5dc0	101	if ( isPPbData) { //Only use two profiles if pPb
	102	const Char_t* periodNames[2] = {"LHC13b", "LHC13c"};
	103	TProfile* multEstimatorAvg[2];
	104	TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
	105	if(!fileEstimator) {
	106	AliFatal("File with multiplicity estimator not found\n");
	107	return;
	108	}
	109
	110	for(Int_t ip=0; ip<2; ip++) {
	111	multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));
	112	}
	113	dMultTask->SetMultiplVsZProfileLHC13b(multEstimatorAvg[0]);
	114	dMultTask->SetMultiplVsZProfileLHC13c(multEstimatorAvg[1]);
	115	dMultTask->SetReferenceMultiplcity(refMult);
	116
	117	}
	118	else {
	119	const Char_t* periodNames[4] = {"LHC10b", "LHC10c", "LHC10d", "LHC10e"};
	120	TProfile* multEstimatorAvg[4];
	121	TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
	122	if(!fileEstimator) {
	123	AliFatal("File with multiplicity estimator not found\n");
	124	return;
	125	}
	126
	127	for(Int_t ip=0; ip<4; ip++) {
	128	multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));
	129	}
	130	dMultTask->SetMultiplVsZProfileLHC10b(multEstimatorAvg[0]);
	131	dMultTask->SetMultiplVsZProfileLHC10c(multEstimatorAvg[1]);
	132	dMultTask->SetMultiplVsZProfileLHC10d(multEstimatorAvg[2]);
	133	dMultTask->SetMultiplVsZProfileLHC10e(multEstimatorAvg[3]);
	134	dMultTask->SetReferenceMultiplcity(refMult);
	135	}
f09fd3fe	136	}
b7fae104	137	mgr->AddTask(dMultTask);
	138
	139	// Create containers for input/output
3d27f7b3	140
b7fae104	141	TString inname = "cinput";
	142	TString outname = "coutput";
	143	TString cutsname = "coutputCuts";
	144	TString normname = "coutputNorm";
6c06d0fe	145	TString profname = "coutputProf";
3d27f7b3	146
	147	inname += Name.Data();
	148	outname += Name.Data();
	149	cutsname += Name.Data();
	150	normname += Name.Data();
6c06d0fe	151	profname += Name.Data();
b7fae104	152	inname += finDirname.Data();
	153	outname += finDirname.Data();
	154	cutsname += finDirname.Data();
	155	normname += finDirname.Data();
6c06d0fe	156	profname += finDirname.Data();
b7fae104	157
3d27f7b3	158	AliAnalysisDataContainer *cinput = mgr->CreateContainer(inname,TChain::Class(),AliAnalysisManager::kInputContainer);
3d27f7b3	159
b7fae104	160	TString outputfile = AliAnalysisManager::GetCommonFileName();
3d27f7b3	161	outputfile += ":PWG3_D2H_DMult_";
	162	outputfile += Name.Data();
	163	outputfile += finDirname.Data();
	164
	165	AliAnalysisDataContainer *coutputCuts = mgr->CreateContainer(cutsname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
	166	AliAnalysisDataContainer *coutput = mgr->CreateContainer(outname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
	167	AliAnalysisDataContainer *coutputNorm = mgr->CreateContainer(normname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
6c06d0fe	168	AliAnalysisDataContainer *coutputProf = mgr->CreateContainer(profname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
b7fae104	169
	170	mgr->ConnectInput(dMultTask,0,mgr->GetCommonInputContainer());
	171
	172	mgr->ConnectOutput(dMultTask,1,coutput);
	173
	174	mgr->ConnectOutput(dMultTask,2,coutputCuts);
	175
	176	mgr->ConnectOutput(dMultTask,3,coutputNorm);
6c06d0fe	177
	178	mgr->ConnectOutput(dMultTask,4,coutputProf);
	179
b7fae104	180	return dMultTask;
b7fae104	181	}