[u/mrichter/AliRoot.git] / PWGHF / hfe / macros / configs / pp / ConfigHFEpid2SYS.C

TF1* GetEtaCorrection(TString listname){
 
  TString etaMap="$ALICE_ROOT/PWGDQ/dielectron/files/EtaCorrMaps.root";
  
  if (gSystem->AccessPathName(gSystem->ExpandPathName(etaMap.Data()))){
    Error("ConfigHFEpid2SYS","Eta map not found: %s",etaMap.Data());
    return 0;
  }

  TFile f(etaMap.Data());
  if (!f.IsOpen()) return 0;
  gROOT->cd();
  TList *keys=f.GetListOfKeys();

  for (Int_t i=0; i<keys->GetEntries(); ++i){
    TString kName=keys->At(i)->GetName();
    TPRegexp reg(kName);
    if (reg.MatchB(listname)){
      printf("Using Eta Correction Function: %s\n",kName.Data());
      return (TF1*)f.Get(kName.Data());
    }
  }
  return 0;
}

AliAnalysisTaskHFE* ConfigHFEpid2SYS(Bool_t useMC, 
				     TString appendix, 
				     UChar_t TPCcl=70, UChar_t TPCclPID = 80, UChar_t ITScl=3, 
				     Double_t DCAxy=1000., Double_t DCAz=1000.,
				     Double_t TPCs=0., Double_t TPCu=3.09, Double_t TOFs=3.,
				     Double_t IpSig=3., Bool_t prodcut = kFALSE, Bool_t IPAbs = kTRUE, Int_t itshitpixel = 0, 
				     Bool_t withetacorrection = kTRUE, TString listname="",
				     Int_t ptbin=0,
				     Bool_t kAnalyseTaggedTracks=kFALSE, Bool_t kMCQA=kFALSE, Bool_t kDEStep=kFALSE,
				     Long_t aodfilter=-1){
  //
  // HFE task configuration PID2 (TOF-TPC only!)
  //

  // Name
  printf("appendix %s\n", appendix.Data());

  // hfecuts
  AliHFEcuts *hfecuts = new AliHFEcuts(Form("hfeCutsPID2_%s",appendix.Data()),"HFE cuts TOF TPC");
  hfecuts->CreateStandardCuts();
  hfecuts->SetMinNClustersTPC(TPCcl);
  hfecuts->SetMinNClustersTPCPID(TPCclPID);
  hfecuts->SetMinNClustersITS(ITScl);
  hfecuts->SetMinRatioTPCclusters(0.6);
  hfecuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
  hfecuts->SetCutITSpixel(itshitpixel);
  hfecuts->SetCheckITSLayerStatus(kFALSE);
  hfecuts->SetIPcutParam(0,0,0,IpSig,kTRUE,IPAbs);
  if(IpSig>100&&IpSig<300){
    hfecuts->SetIPcutParam(0.0064,0.078,-0.56,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077*exp(-0.65*pt))
  }
  else if(IpSig>300&&IpSig<320){
    hfecuts->SetIPcutParam(0.0044,0.078,-0.56,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077*exp(-0.65*pt))
  }
  else if(IpSig>320&&IpSig<350){
    hfecuts->SetIPcutParam(0.0054,0.078,-0.56,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077*exp(-0.65*pt))
  }
  else if(IpSig>350&&IpSig<500){
    hfecuts->SetIPcutParam(0.011,0.077,-0.65,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077*exp(-0.65*pt))
  }
  else if(IpSig>500&&IpSig<700){
    hfecuts->SetIPcutParam(0.012,0.077,-0.65,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077*exp(-0.65*pt))
  }
  else if(IpSig>700&&IpSig<900){
    hfecuts->SetIPcutParam(0.013,0.077,-0.65,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077*exp(-0.65*pt))
  }

  if(prodcut) hfecuts->SetProductionVertex(0,100,0,100);
  else {
    if((itshitpixel==AliHFEextraCuts::kAny) || (itshitpixel==AliHFEextraCuts::kSecond)) hfecuts->SetProductionVertex(0,7,0,7);
  }
  
  hfecuts->SetMaxImpactParam(DCAxy,DCAz);
  hfecuts->SetTOFPIDStep(kTRUE);
  hfecuts->SetUseMixedVertex(kTRUE);
  hfecuts->SetVertexRange(10.);

  // analysis task
  AliAnalysisTaskHFE *task = new AliAnalysisTaskHFE(Form("HFEanalysisPID2_%s",appendix.Data()));
  task->SetHFECuts(hfecuts);
  task->SetRemovePileUp(kTRUE);
  task->GetPIDQAManager()->SetHighResolutionHistos();
  if(useMC) task->SetHasMCData(kTRUE); // necessary for AOD
  printf("AOD filter %d On/OFF?\n",aodfilter);
  if(aodfilter > 0) {
    printf("ON AOD filter %d\n",aodfilter);
    task->SetUseFilterAOD(kTRUE);
    task->SetFilter(aodfilter);
  }


  // Define Variables
  if(ptbin==1){
    Double_t ptbinning[19] = {0., 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 2., 2.5, 3., 4., 5., 6., 8., 12., 16., 20.};
  }
  else{
    Double_t ptbinning[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
  }
  Double_t etabinning[17] = {-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};

  Int_t sizept=(sizeof(ptbinning)/sizeof(double))-1;
  Int_t sizeeta=(sizeof(etabinning)/sizeof(double))-1;

  AliHFEvarManager *vm = task->GetVarManager();
  vm->AddVariable("pt", sizept, ptbinning);
  vm->AddVariable("eta", sizeeta, -0.8,0.8);
  vm->AddVariable("phi",21, -0, 2*TMath::Pi());
  vm->AddVariable("charge");
  vm->AddVariable("source");
  
  // Contamination
  if(!useMC){
    // New background model (LHC10d pass2)
    TF1 *hBackground = new TF1("hadronicBackgroundFunction", "TMath::Exp(([0]/(x**1.5))+[1])", 0., 20.);
    // These settings assume that the default is a cut on .ge.120 TPC clusters (Sep 27, 2011)
    hBackground->SetParameter(0, -55.18);
    hBackground->SetParameter(1, -0.0026);
    if (TPCcl == 100){
      hBackground->SetParameter(0, -39.5);
      hBackground->SetParameter(1, -0.438);
    } elseif (TPCcl == 140){
      hBackground->SetParameter(0, -82.11);
      hBackground->SetParameter(1, 1.138);
    } 

    task->SetBackGroundFactorsFunction(hBackground);
  }

  // Define PID
  AliHFEpid *pid = task->GetPID();
  if(useMC) pid->SetHasMCData(kTRUE);
  pid->AddDetector("TOF", 0);
  pid->AddDetector("TPC", 1);
  // HERE PUT THE STRAIGHT CUT
  Double_t params[4];
  char *cutmodel;
  if(useMC){
    // Monte-Carlo needs modelling of the falling mean with momentum at low momentum
    // for high momentum it is consistent with a flat -0.94
    cutmodel = "expo(0)+pol1(2)";//[0]*TMath::Exp([1]*x) + [2] + [3]*x";
    Double_t paramsMC[4] = {-1.00625e-01, -2.09446e+00, -4.71247e-01, 1.80301e-02};
    for(int ipar = 0; ipar < 4; ipar++) params[ipar] = paramsMC[ipar];
  } else {
    // Data is consistent with a flat constant: (Sep 27, 2011)
    // 100 clusters: mean = -0.076, width = 1.035
    // 120 clusters: mean = -0.113, width = 1.03
    // 140 clusters: mean = -0.093, width = 1.004
    cutmodel = "pol0(0)";
    params[0] = TPCs;

    if(withetacorrection) {
      // Apply eta correction
      AliHFEpidTPC *tpcpid = pid->GetDetPID(AliHFEpid::kTPCpid);
      TF1 *etacorrection = GetEtaCorrection(listname);
      if(etacorrection) tpcpid->SetEtaCorrection(etacorrection);
    }
  }
  pid->ConfigureTOF(TOFs);
  pid->ConfigureTPCdefaultCut(cutmodel, params, TPCu);

  if(kAnalyseTaggedTracks){
    AliHFEcuts *v0trackCuts = new AliHFEcuts("V0trackCuts", "Track Cuts for tagged track Analysis");
    v0trackCuts->CreateStandardCuts();
    v0trackCuts->SetMinNClustersTPC(TPCcl);  
    v0trackCuts->SetMinRatioTPCclusters(0.6);
    v0trackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
    v0trackCuts->SetMinNClustersITS(1);
    v0trackCuts->SetCutITSpixel(AliHFEextraCuts::kFirst);
    v0trackCuts->SetCheckITSLayerStatus(kFALSE);
    v0trackCuts->UnsetVertexRequirement();
    v0trackCuts->SetTOFPIDStep(kTRUE);
    v0trackCuts->SetQAOn();

    task->SwitchOnPlugin(AliAnalysisTaskHFE::kTaggedTrackAnalysis);
    task->SetTaggedTrackCuts(v0trackCuts);
    task->SetCleanTaggedTrack(kTRUE);
  }

  // QA
  printf("task %p\n", task);
  task->SetQAOn(AliAnalysisTaskHFE::kPIDqa);
  if(kMCQA) task->SetQAOn(AliAnalysisTaskHFE::kMCqa);    
  //task->SwitchOnPlugin(AliAnalysisTaskHFE::kIsElecBackGround);
  //task->SwitchOnPlugin(AliAnalysisTaskHFE::kSecVtx);
  if(kDEStep) task->SwitchOnPlugin(AliAnalysisTaskHFE::kDEstep);
  
  printf("*************************************\n");
  printf("Configuring standard Task:\n");
  task->PrintStatus();
  pid->PrintStatus();
  printf("*************************************\n"); 
  return task;
}
Commit	Line	Data
eeb1e952	1	TF1* GetEtaCorrection(TString listname){
	2
	3	TString etaMap="$ALICE_ROOT/PWGDQ/dielectron/files/EtaCorrMaps.root";
	4
	5	if (gSystem->AccessPathName(gSystem->ExpandPathName(etaMap.Data()))){
	6	Error("ConfigHFEpid2SYS","Eta map not found: %s",etaMap.Data());
	7	return 0;
	8	}
	9
	10	TFile f(etaMap.Data());
	11	if (!f.IsOpen()) return 0;
	12	gROOT->cd();
	13	TList *keys=f.GetListOfKeys();
	14
	15	for (Int_t i=0; i<keys->GetEntries(); ++i){
	16	TString kName=keys->At(i)->GetName();
	17	TPRegexp reg(kName);
	18	if (reg.MatchB(listname)){
	19	printf("Using Eta Correction Function: %s\n",kName.Data());
	20	return (TF1*)f.Get(kName.Data());
	21	}
	22	}
	23	return 0;
	24	}
	25
	26	AliAnalysisTaskHFE* ConfigHFEpid2SYS(Bool_t useMC,
	27	TString appendix,
	28	UChar_t TPCcl=70, UChar_t TPCclPID = 80, UChar_t ITScl=3,
d84c68e9	29	Double_t DCAxy=1000., Double_t DCAz=1000.,
d84c68e9	30	Double_t TPCs=0., Double_t TPCu=3.09, Double_t TOFs=3.,
eeb1e952	31	Double_t IpSig=3., Bool_t prodcut = kFALSE, Bool_t IPAbs = kTRUE, Int_t itshitpixel = 0,
31ae1571	32	Bool_t withetacorrection = kTRUE, TString listname="",
eeb1e952	33	Int_t ptbin=0,
31ae1571	34	Bool_t kAnalyseTaggedTracks=kFALSE, Bool_t kMCQA=kFALSE, Bool_t kDEStep=kFALSE,
31ae1571	35	Long_t aodfilter=-1){
d84c68e9	36	//
	37	// HFE task configuration PID2 (TOF-TPC only!)
	38	//
	39
eeb1e952	40	// Name
	41	printf("appendix %s\n", appendix.Data());
	42
	43	// hfecuts
	44	AliHFEcuts *hfecuts = new AliHFEcuts(Form("hfeCutsPID2_%s",appendix.Data()),"HFE cuts TOF TPC");
d84c68e9	45	hfecuts->CreateStandardCuts();
	46	hfecuts->SetMinNClustersTPC(TPCcl);
	47	hfecuts->SetMinNClustersTPCPID(TPCclPID);
	48	hfecuts->SetMinNClustersITS(ITScl);
	49	hfecuts->SetMinRatioTPCclusters(0.6);
	50	hfecuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
eeb1e952	51	hfecuts->SetCutITSpixel(itshitpixel);
d84c68e9	52	hfecuts->SetCheckITSLayerStatus(kFALSE);
eeb1e952	53	hfecuts->SetIPcutParam(0,0,0,IpSig,kTRUE,IPAbs);
	54	if(IpSig>100&&IpSig<300){
	55	hfecuts->SetIPcutParam(0.0064,0.078,-0.56,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077exp(-0.65pt))
	56	}
	57	else if(IpSig>300&&IpSig<320){
	58	hfecuts->SetIPcutParam(0.0044,0.078,-0.56,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077exp(-0.65pt))
	59	}
	60	else if(IpSig>320&&IpSig<350){
	61	hfecuts->SetIPcutParam(0.0054,0.078,-0.56,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077exp(-0.65pt))
	62	}
	63	else if(IpSig>350&&IpSig<500){
	64	hfecuts->SetIPcutParam(0.011,0.077,-0.65,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077exp(-0.65pt))
	65	}
	66	else if(IpSig>500&&IpSig<700){
	67	hfecuts->SetIPcutParam(0.012,0.077,-0.65,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077exp(-0.65pt))
	68	}
	69	else if(IpSig>700&&IpSig<900){
	70	hfecuts->SetIPcutParam(0.013,0.077,-0.65,0,kFALSE,IPAbs); // used Carlo's old parameter (new: 0.011+0.077exp(-0.65pt))
	71	}
d84c68e9	72
eeb1e952	73	if(prodcut) hfecuts->SetProductionVertex(0,100,0,100);
	74	else {
	75	if((itshitpixel==AliHFEextraCuts::kAny) \|\| (itshitpixel==AliHFEextraCuts::kSecond)) hfecuts->SetProductionVertex(0,7,0,7);
	76	}
	77
d84c68e9	78	hfecuts->SetMaxImpactParam(DCAxy,DCAz);
eeb1e952	79	hfecuts->SetTOFPIDStep(kTRUE);
d84c68e9	80	hfecuts->SetUseMixedVertex(kTRUE);
	81	hfecuts->SetVertexRange(10.);
	82
eeb1e952	83	// analysis task
eeb1e952	84	AliAnalysisTaskHFE *task = new AliAnalysisTaskHFE(Form("HFEanalysisPID2_%s",appendix.Data()));
d84c68e9	85	task->SetHFECuts(hfecuts);
	86	task->SetRemovePileUp(kTRUE);
	87	task->GetPIDQAManager()->SetHighResolutionHistos();
31ae1571	88	if(useMC) task->SetHasMCData(kTRUE); // necessary for AOD
a86a8fda	89	printf("AOD filter %d On/OFF?\n",aodfilter);
31ae1571	90	if(aodfilter > 0) {
a86a8fda	91	printf("ON AOD filter %d\n",aodfilter);
	92	task->SetUseFilterAOD(kTRUE);
	93	task->SetFilter(aodfilter);
31ae1571	94	}
31ae1571	95
d84c68e9	96
d84c68e9	97	// Define Variables
eeb1e952	98	if(ptbin==1){
	99	Double_t ptbinning[19] = {0., 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 2., 2.5, 3., 4., 5., 6., 8., 12., 16., 20.};
	100	}
	101	else{
	102	Double_t ptbinning[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
	103	}
d84c68e9	104	Double_t etabinning[17] = {-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};
d84c68e9	105
eeb1e952	106	Int_t sizept=(sizeof(ptbinning)/sizeof(double))-1;
	107	Int_t sizeeta=(sizeof(etabinning)/sizeof(double))-1;
	108
d84c68e9	109	AliHFEvarManager *vm = task->GetVarManager();
eeb1e952	110	vm->AddVariable("pt", sizept, ptbinning);
	111	vm->AddVariable("eta", sizeeta, -0.8,0.8);
	112	vm->AddVariable("phi",21, -0, 2*TMath::Pi());
d84c68e9	113	vm->AddVariable("charge");
d84c68e9	114	vm->AddVariable("source");
eeb1e952	115
eeb1e952	116	// Contamination
d84c68e9	117	if(!useMC){
	118	// New background model (LHC10d pass2)
	119	TF1 hBackground = new TF1("hadronicBackgroundFunction", "TMath::Exp(([0]/(x*1.5))+[1])", 0., 20.);
	120	// These settings assume that the default is a cut on .ge.120 TPC clusters (Sep 27, 2011)
	121	hBackground->SetParameter(0, -55.18);
	122	hBackground->SetParameter(1, -0.0026);
	123	if (TPCcl == 100){
	124	hBackground->SetParameter(0, -39.5);
	125	hBackground->SetParameter(1, -0.438);
	126	} elseif (TPCcl == 140){
	127	hBackground->SetParameter(0, -82.11);
	128	hBackground->SetParameter(1, 1.138);
	129	}
	130
	131	task->SetBackGroundFactorsFunction(hBackground);
	132	}
	133
	134	// Define PID
	135	AliHFEpid *pid = task->GetPID();
	136	if(useMC) pid->SetHasMCData(kTRUE);
	137	pid->AddDetector("TOF", 0);
	138	pid->AddDetector("TPC", 1);
	139	// HERE PUT THE STRAIGHT CUT
	140	Double_t params[4];
	141	char *cutmodel;
	142	if(useMC){
	143	// Monte-Carlo needs modelling of the falling mean with momentum at low momentum
	144	// for high momentum it is consistent with a flat -0.94
	145	cutmodel = "expo(0)+pol1(2)";//[0]TMath::Exp([1]x) + [2] + [3]*x";
	146	Double_t paramsMC[4] = {-1.00625e-01, -2.09446e+00, -4.71247e-01, 1.80301e-02};
	147	for(int ipar = 0; ipar < 4; ipar++) params[ipar] = paramsMC[ipar];
	148	} else {
	149	// Data is consistent with a flat constant: (Sep 27, 2011)
	150	// 100 clusters: mean = -0.076, width = 1.035
	151	// 120 clusters: mean = -0.113, width = 1.03
	152	// 140 clusters: mean = -0.093, width = 1.004
	153	cutmodel = "pol0(0)";
	154	params[0] = TPCs;
eeb1e952	155
	156	if(withetacorrection) {
	157	// Apply eta correction
	158	AliHFEpidTPC *tpcpid = pid->GetDetPID(AliHFEpid::kTPCpid);
	159	TF1 *etacorrection = GetEtaCorrection(listname);
	160	if(etacorrection) tpcpid->SetEtaCorrection(etacorrection);
	161	}
d84c68e9	162	}
	163	pid->ConfigureTOF(TOFs);
	164	pid->ConfigureTPCdefaultCut(cutmodel, params, TPCu);
	165
	166	if(kAnalyseTaggedTracks){
	167	AliHFEcuts *v0trackCuts = new AliHFEcuts("V0trackCuts", "Track Cuts for tagged track Analysis");
	168	v0trackCuts->CreateStandardCuts();
	169	v0trackCuts->SetMinNClustersTPC(TPCcl);
	170	v0trackCuts->SetMinRatioTPCclusters(0.6);
	171	v0trackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
	172	v0trackCuts->SetMinNClustersITS(1);
eeb1e952	173	v0trackCuts->SetCutITSpixel(AliHFEextraCuts::kFirst);
d84c68e9	174	v0trackCuts->SetCheckITSLayerStatus(kFALSE);
d84c68e9	175	v0trackCuts->UnsetVertexRequirement();
d84c68e9	176	v0trackCuts->SetTOFPIDStep(kTRUE);
	177	v0trackCuts->SetQAOn();
	178
	179	task->SwitchOnPlugin(AliAnalysisTaskHFE::kTaggedTrackAnalysis);
	180	task->SetTaggedTrackCuts(v0trackCuts);
	181	task->SetCleanTaggedTrack(kTRUE);
	182	}
	183
	184	// QA
	185	printf("task %p\n", task);
	186	task->SetQAOn(AliAnalysisTaskHFE::kPIDqa);
eeb1e952	187	if(kMCQA) task->SetQAOn(AliAnalysisTaskHFE::kMCqa);
d84c68e9	188	//task->SwitchOnPlugin(AliAnalysisTaskHFE::kIsElecBackGround);
d84c68e9	189	//task->SwitchOnPlugin(AliAnalysisTaskHFE::kSecVtx);
eeb1e952	190	if(kDEStep) task->SwitchOnPlugin(AliAnalysisTaskHFE::kDEstep);
eeb1e952	191
d84c68e9	192	printf("*************************************\n");
	193	printf("Configuring standard Task:\n");
	194	task->PrintStatus();
	195	pid->PrintStatus();
	196	printf("*************************************\n");
	197	return task;
	198	}