[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / macros / makeTFile4CutsLctopKpi.C

#include <Riostream.h>
#include <TFile.h>
#include <AliRDHFCutsLctopKpi.h>
#include <AliAODPidHF.h>
#include <TClonesArray.h>
#include <TParameter.h>


//Use:
//Set hard coded commentet with //set this!!
// root[] .L makeInput...C++
// root[] makeInputAliAnalysisTaskSE...()
//similar macros for the other D mesons

//Author: Rosa Romita, r.romita@gsi.de


//macro to make a .root file which contains an AliRDHFCutsLctopKpi for AliAnalysisTaskSELambdac task

void makeInputAliAnalysisTaskSELctopKpi(){

  AliRDHFCutsLctopKpi* RDHFLctopKpiProd=new AliRDHFCutsLctopKpi();
  RDHFLctopKpiProd->SetName("LctopKpiProdCuts");
  RDHFLctopKpiProd->SetTitle("Production cuts for Lc analysis");

  AliRDHFCutsLctopKpi* RDHFLctopKpiAn=new AliRDHFCutsLctopKpi();
  RDHFLctopKpiAn->SetName("LctopKpiAnalysisCuts");
  RDHFLctopKpiAn->SetTitle("Analysis cuts for Lc analysis");

  AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
  //default
  esdTrackCuts->SetRequireTPCRefit(kTRUE);
  esdTrackCuts->SetRequireITSRefit(kTRUE);
  esdTrackCuts->SetMinNClustersITS(4); // default is 5
  esdTrackCuts->SetMinNClustersTPC(70);
  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
					 AliESDtrackCuts::kAny); 
 // default is kBoth, otherwise kAny
  esdTrackCuts->SetMinDCAToVertexXY(0.);
  esdTrackCuts->SetPtRange(0.3,1.e10);


  RDHFLctopKpiProd->AddTrackCuts(esdTrackCuts);
  RDHFLctopKpiAn->AddTrackCuts(esdTrackCuts);

  const Int_t nvars=13;

  const Int_t nptbins=4;
  Float_t* ptbins;
  ptbins=new Float_t[nptbins+1];
  ptbins[0]=0.;
  ptbins[1]=2.;
  ptbins[2]=3.;
  ptbins[3]=4.;
  ptbins[4]=12.;
  

  Float_t** prodcutsval;
  prodcutsval=new Float_t*[nvars];
  for(Int_t iv=0;iv<nvars;iv++){
    prodcutsval[iv]=new Float_t[nptbins];
  }

  //  inv. mass [GeV]
  // pTK [GeV/c]
  // pTPi [GeV/c]
  // d0K [cm]   lower limit!
  // d0Pi [cm]  lower limit!
  // dist12 (cm)
  // sigmavert (cm)
  // dist prim-sec (cm)
  // pM=Max{pT1,pT2,pT3} (GeV/c)
  // cosThetaPoint
  // Sum d0^2 (cm^2)
  // dca cut (cm)
  for(Int_t ipt=0;ipt<nptbins;ipt++){
    prodcutsval[0][ipt]=0.18;
    prodcutsval[1][ipt]=0.4;
    prodcutsval[2][ipt]=0.5;
    prodcutsval[3][ipt]=0.;
    prodcutsval[4][ipt]=0.;
    prodcutsval[5][ipt]=0.01;
    prodcutsval[6][ipt]=0.06;
    prodcutsval[7][ipt]=0.005;
    prodcutsval[8][ipt]=0.7;
    prodcutsval[9][ipt]=0.;
    prodcutsval[10][ipt]=0.;
    prodcutsval[11][ipt]=0.05;
    prodcutsval[12][ipt]=0.4;
  }

  RDHFLctopKpiProd->SetPtBins(nptbins+1,ptbins);
  RDHFLctopKpiProd->SetCuts(nvars,nptbins,prodcutsval);

  Float_t** anacutsval;
  anacutsval=new Float_t*[nvars];
  for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
  for(Int_t ipt2=0;ipt2<nptbins;ipt2++){
   anacutsval[0][ipt2]=0.18;
   anacutsval[1][ipt2]=0.4;
   anacutsval[2][ipt2]=0.5;
   anacutsval[3][ipt2]=0.;
   anacutsval[4][ipt2]=0.;
   anacutsval[5][ipt2]=0.01;
   anacutsval[6][ipt2]=0.06;
   anacutsval[7][ipt2]=0.005;
   anacutsval[8][ipt2]=0.7;
   anacutsval[9][ipt2]=0.;
   anacutsval[10][ipt2]=0.;
   anacutsval[11][ipt2]=0.05;
   anacutsval[12][ipt2]=0.4;
  }

  // pt kaon
  anacutsval[1][0]=0.5;
  anacutsval[1][1]=0.85;
  anacutsval[1][2]=0.9;
  anacutsval[1][3]=0.4;
  //pt proton
  anacutsval[2][0]=0.5;
  anacutsval[2][1]=0.6;
  anacutsval[2][2]=0.9;
  anacutsval[2][3]=0.9;

  //pt pion
  anacutsval[12][0]=0.475;
  anacutsval[12][1]=0.75;
  anacutsval[12][2]=0.75;
  anacutsval[12][3]=0.7;

  anacutsval[5][0]=0.02;
  anacutsval[5][1]=0.025;
  anacutsval[5][2]=0.02;
  anacutsval[5][3]=0.01;

  anacutsval[7][0]=0.00625;
  anacutsval[7][1]=0.0125;
  anacutsval[7][2]=0.005;
  anacutsval[7][3]=0.007;

  anacutsval[9][0]=0.5;
  anacutsval[9][1]=0.2;
  anacutsval[9][2]=0.6;
  anacutsval[9][3]=0.;

  anacutsval[10][0]=0.00125;

  RDHFLctopKpiAn->SetPtBins(nptbins+1,ptbins);
  RDHFLctopKpiAn->SetCuts(nvars,nptbins,anacutsval);

 //  RDHFLc->SetRecoKF(); //set this if you want to recompute the secondary vertex with the KF package
    //pid settings
    //1. kaon: default one
  AliAODPidHF* pidObjK=new AliAODPidHF();
  Double_t sigmasK[5]={3.,1.,1.,3.,2.};
  pidObjK->SetSigma(sigmasK);
  pidObjK->SetAsym(kTRUE);
  pidObjK->SetMatch(1);
  pidObjK->SetTPC(kTRUE);
  pidObjK->SetTOF(kTRUE);
  pidObjK->SetITS(kTRUE);
  Double_t plimK[2]={0.5,0.8};
  pidObjK->SetPLimit(plimK,2);
  pidObjK->SetTOFdecide(kTRUE);
  
  RDHFLctopKpiProd->SetPidHF(pidObjK);
  RDHFLctopKpiAn->SetPidHF(pidObjK);

    //2. pion 
  AliAODPidHF* pidObjpi=new AliAODPidHF();
  pidObjpi->SetTPC(kTRUE);
  Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
  pidObjpi->SetSigma(sigmaspi);
//  pidObjpi->SetTOFdecide(kTRUE);

  RDHFLctopKpiProd->SetPidpion(pidObjpi);
  RDHFLctopKpiAn->SetPidpion(pidObjpi);

  // 3. proton
  AliAODPidHF* pidObjp=new AliAODPidHF();
  Double_t sigmasp[5]={3.,1.,1.,3.,2.};
  pidObjp->SetSigma(sigmasp);
  pidObjp->SetAsym(kTRUE);
  pidObjp->SetMatch(1);
  pidObjp->SetTPC(kTRUE);
  pidObjp->SetTOF(kTRUE);
  pidObjp->SetITS(kTRUE);
  Double_t plimp[2]={1.,2.};
  pidObjp->SetPLimit(plimp,2);
  pidObjp->SetTOFdecide(kTRUE);

  RDHFLctopKpiProd->SetPidprot(pidObjp);
  RDHFLctopKpiAn->SetPidprot(pidObjp);

  Bool_t pidflag=kTRUE;
  RDHFLctopKpiAn->SetUsePID(pidflag);
  RDHFLctopKpiProd->SetUsePID(pidflag);
  if(pidflag) cout<<"PID is used"<<endl;
  else cout<<"PID is not used"<<endl;

  cout<<"This is the object I'm going to save:"<<endl;
  RDHFLctopKpiProd->PrintAll();
  RDHFLctopKpiAn->PrintAll();
  TFile* fout=new TFile("cuts4LctopKpi.root","RECREATE"); 
  fout->cd();
  RDHFLctopKpiProd->Write();
  RDHFLctopKpiAn->Write();
  fout->Close();
  delete fout;
  delete pidObjp;
  delete pidObjpi;
  delete pidObjK;
  delete RDHFLctopKpiProd;
  delete RDHFLctopKpiAn;

}


//macro to make a .root file (for significance maximization) which contains an AliRDHFCutsLctopKpi with loose set of cuts  and TParameter with the tighest value of these cuts

void makeInputAliAnalysisTaskSESignificanceMaximization(){

  AliRDHFCutsLctopKpi* RDHFLctopKpi=new AliRDHFCutsLctopKpi();
  RDHFLctopKpi->SetName("loosercuts");
  RDHFLctopKpi->SetTitle("Cuts for significance maximization");

  AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
  //default
  esdTrackCuts->SetRequireTPCRefit(kTRUE);
  esdTrackCuts->SetMinNClustersTPC(70);
  esdTrackCuts->SetRequireITSRefit(kTRUE);
  esdTrackCuts->SetMinNClustersITS(4);

  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
  esdTrackCuts->SetMinDCAToVertexXY(0.);
  esdTrackCuts->SetEtaRange(-0.8,0.8);
  esdTrackCuts->SetPtRange(0.3,1.e10);
  
  RDHFLctopKpi->AddTrackCuts(esdTrackCuts);

  const Int_t nvars=13;

  const Int_t nptbins=4; //change this when adding pt bins!
  Float_t ptbins[nptbins+1];
  ptbins[0]=0.;
  ptbins[1]=2.;
  ptbins[2]=3.;
  ptbins[3]=4.;
  ptbins[4]=12.;

  RDHFLctopKpi->SetPtBins(nptbins+1,ptbins);

  Float_t** rdcutsvalmine;
  rdcutsvalmine=new Float_t*[nvars];
  for(Int_t iv=0;iv<nvars;iv++){
    rdcutsvalmine[iv]=new Float_t[nptbins];
  }

  //setting my cut values
  //  inv. mass [GeV]
  // pTK [GeV/c]
  // pTP [GeV/c]
  // d0K [cm]   lower limit!
  // d0Pi [cm]  lower limit!
  // dist12 (cm)
  // sigmavert (cm)
  // dist prim-sec (cm)
  // pM=Max{pT1,pT2,pT3} (GeV/c)
  // cosThetaPoint
  // Sum d0^2 (cm^2)
  // dca cut (cm)
  // pt pion
  Float_t cutsMatrixLctopKpiStand[nptbins][nvars]=
  {{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4},
   {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4},
   {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4},
   {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4}};

  //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
  Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
  for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
  for (Int_t ibin=0;ibin<nptbins;ibin++){
    for (Int_t ivar = 0; ivar<nvars; ivar++){
      cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixLctopKpiStand[ibin][ivar];
    }
  }
  RDHFLctopKpi->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);


  Int_t nvarsforopt=RDHFLctopKpi->GetNVarsForOpt();
  Int_t dim=4; //set this!!
  Bool_t *boolforopt;
  boolforopt=new Bool_t[nvars];
  if(dim>nvarsforopt){
    cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
    return;
  } else {
    if(dim==nvarsforopt){
      boolforopt=RDHFLctopKpi->GetVarsForOpt();
    }else{
      TString *names;
      names=new TString[nvars];
      TString answer="";
      Int_t checktrue=0;
      names=RDHFLctopKpi->GetVarNames();
      for(Int_t i=0;i<nvars;i++){
	cout<<names[i]<<" for opt? (y/n)"<<endl;
	cin>>answer;
	if(answer=="y") {
	  boolforopt[i]=kTRUE;
	  checktrue++;
	}
	else boolforopt[i]=kFALSE;
      }
      if (checktrue!=dim) {
	cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
	return;
      }
      RDHFLctopKpi->SetVarsForOpt(dim,boolforopt);
    }
  }


  Float_t tighterval[dim][nptbins];
  // 5 dist12 (cm)        <---
  // 7 dist prim-sec (cm) <---
  // 9 cosThetaPoint      <---
  // 10 Sum d0^2 (cm^2)   <---
  // 11 dca cut (cm)      
  // {0.18,0.4,0.5,0.,0.,0.01 <--- ,0.06,0.005 <--- ,0.7,0. <--- ,0. <--- ,0.05 }

  //number of steps for each variable is set in the AddTask arguments (default=8)
  //set this!!
  for(Int_t ipt=0;ipt<nptbins;ipt++){
    tighterval[0][ipt]=0.03;
    tighterval[1][ipt]=0.02;
    tighterval[2][ipt]=1.;
    tighterval[3][ipt]=0.01;
  }


  TString name=""; 
  Int_t arrdim=dim*nptbins;
  cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
  TClonesArray max("TParameter<float>",arrdim);
  for(Int_t ival=0;ival<dim;ival++){
    for(Int_t jpt=0;jpt<nptbins;jpt++){
      name=Form("par%dptbin%d",ival,jpt);
      cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
      new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
    }
  }

  Bool_t flagPID=kFALSE;
  RDHFLctopKpi->SetUsePID(flagPID);

  RDHFLctopKpi->PrintAll();
  printf("Use PID? %s\n",flagPID ? "yes" : "no");

    //pid settings
    //1. kaon: default one
  AliAODPidHF* pidObjK=new AliAODPidHF();
  Double_t sigmasK[5]={3.,1.,1.,3.,2.};
  pidObjK->SetSigma(sigmasK);
  pidObjK->SetAsym(kTRUE);
  pidObjK->SetMatch(1);
  pidObjK->SetTPC(kTRUE);
  pidObjK->SetTOF(kTRUE);
  pidObjK->SetITS(kTRUE);
  Double_t plimK[2]={0.5,0.8};
  pidObjK->SetPLimit(plimK,2);
  pidObjK->SetTOFdecide(kTRUE);
  
  RDHFLctopKpi->SetPidHF(pidObjK);

    //2. pion 
  AliAODPidHF* pidObjpi=new AliAODPidHF();
  pidObjpi->SetTPC(kTRUE);
  Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
  pidObjpi->SetSigma(sigmaspi);
  pidObjpi->SetTOFdecide(kTRUE);

  RDHFLctopKpi->SetPidpion(pidObjpi);

  // 3. proton
  AliAODPidHF* pidObjp=new AliAODPidHF();
  Double_t sigmasp[5]={3.,1.,1.,3.,2.};
  pidObjp->SetSigma(sigmasp);
  pidObjp->SetAsym(kTRUE);
  pidObjp->SetMatch(1);
  pidObjp->SetTPC(kTRUE);
  pidObjp->SetTOF(kTRUE);
  pidObjp->SetITS(kTRUE);
  Double_t plimp[2]={1.,2.};
  pidObjp->SetPLimit(plimp,2);
  pidObjp->SetTOFdecide(kTRUE);

  RDHFLctopKpi->SetPidprot(pidObjp);

  //activate pileup rejection (for pp)
  //RDHFLctopKpi->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);

  //Do not recalculate the vertex
  RDHFLctopKpi->SetRemoveDaughtersFromPrim(kFALSE); //activate for pp

  TString cent="";
  //centrality selection (Pb-Pb)
  Float_t minc=20,maxc=80;
  RDHFLctopKpi->SetMinCentrality(minc);
  RDHFLctopKpi->SetMaxCentrality(maxc);
  cent=Form("%.0f%.0f",minc,maxc);
  //  RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
  RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid

  //temporary
  //  RDHFLctopKpi->SetFixRefs();

  TFile* fout=new TFile(Form("cuts4SignifMaxim%s%s%sRecVtx%sPileupRej.root", RDHFLctopKpi->GetUseCentrality()==0 ? "pp" : "PbPb",cent.Data(),RDHFLctopKpi->GetIsPrimaryWithoutDaughters() ? "" : "No",RDHFLctopKpi->GetOptPileUp() ? "" : "No"),"recreate");   //set this!! 

  fout->cd();
  RDHFLctopKpi->Write();
  max.Write();
  fout->Close();
 
}
Commit	Line	Data
7ad4b782	1	#include <Riostream.h>
	2	#include <TFile.h>
	3	#include <AliRDHFCutsLctopKpi.h>
	4	#include <AliAODPidHF.h>
	5	#include <TClonesArray.h>
	6	#include <TParameter.h>
	7
	8
	9	//Use:
	10	//Set hard coded commentet with //set this!!
	11	// root[] .L makeInput...C++
	12	// root[] makeInputAliAnalysisTaskSE...()
	13	//similar macros for the other D mesons
	14
	15	//Author: Rosa Romita, r.romita@gsi.de
	16
	17
	18	//macro to make a .root file which contains an AliRDHFCutsLctopKpi for AliAnalysisTaskSELambdac task
	19
	20	void makeInputAliAnalysisTaskSELctopKpi(){
	21
	22	AliRDHFCutsLctopKpi* RDHFLctopKpiProd=new AliRDHFCutsLctopKpi();
	23	RDHFLctopKpiProd->SetName("LctopKpiProdCuts");
	24	RDHFLctopKpiProd->SetTitle("Production cuts for Lc analysis");
	25
	26	AliRDHFCutsLctopKpi* RDHFLctopKpiAn=new AliRDHFCutsLctopKpi();
	27	RDHFLctopKpiAn->SetName("LctopKpiAnalysisCuts");
	28	RDHFLctopKpiAn->SetTitle("Analysis cuts for Lc analysis");
	29
	30	AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
	31	esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
	32	//default
	33	esdTrackCuts->SetRequireTPCRefit(kTRUE);
	34	esdTrackCuts->SetRequireITSRefit(kTRUE);
	35	esdTrackCuts->SetMinNClustersITS(4); // default is 5
	36	esdTrackCuts->SetMinNClustersTPC(70);
	37	esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
	38	AliESDtrackCuts::kAny);
	39	// default is kBoth, otherwise kAny
	40	esdTrackCuts->SetMinDCAToVertexXY(0.);
	41	esdTrackCuts->SetPtRange(0.3,1.e10);
	42
	43
	44	RDHFLctopKpiProd->AddTrackCuts(esdTrackCuts);
	45	RDHFLctopKpiAn->AddTrackCuts(esdTrackCuts);
	46
a7b533b7	47	const Int_t nvars=13;
7ad4b782	48
	49	const Int_t nptbins=4;
	50	Float_t* ptbins;
	51	ptbins=new Float_t[nptbins+1];
	52	ptbins[0]=0.;
	53	ptbins[1]=2.;
	54	ptbins[2]=3.;
	55	ptbins[3]=4.;
a7b533b7	56	ptbins[4]=12.;
7ad4b782	57
	58
	59	Float_t** prodcutsval;
	60	prodcutsval=new Float_t*[nvars];
	61	for(Int_t iv=0;iv<nvars;iv++){
	62	prodcutsval[iv]=new Float_t[nptbins];
	63	}
	64
5fd38310	65	// inv. mass [GeV]
	66	// pTK [GeV/c]
	67	// pTPi [GeV/c]
	68	// d0K [cm] lower limit!
	69	// d0Pi [cm] lower limit!
	70	// dist12 (cm)
	71	// sigmavert (cm)
	72	// dist prim-sec (cm)
	73	// pM=Max{pT1,pT2,pT3} (GeV/c)
	74	// cosThetaPoint
	75	// Sum d0^2 (cm^2)
	76	// dca cut (cm)
7ad4b782	77	for(Int_t ipt=0;ipt<nptbins;ipt++){
	78	prodcutsval[0][ipt]=0.18;
	79	prodcutsval[1][ipt]=0.4;
	80	prodcutsval[2][ipt]=0.5;
	81	prodcutsval[3][ipt]=0.;
	82	prodcutsval[4][ipt]=0.;
	83	prodcutsval[5][ipt]=0.01;
	84	prodcutsval[6][ipt]=0.06;
	85	prodcutsval[7][ipt]=0.005;
5fd38310	86	prodcutsval[8][ipt]=0.7;
7ad4b782	87	prodcutsval[9][ipt]=0.;
	88	prodcutsval[10][ipt]=0.;
	89	prodcutsval[11][ipt]=0.05;
a7b533b7	90	prodcutsval[12][ipt]=0.4;
7ad4b782	91	}
	92
	93	RDHFLctopKpiProd->SetPtBins(nptbins+1,ptbins);
	94	RDHFLctopKpiProd->SetCuts(nvars,nptbins,prodcutsval);
	95
	96	Float_t** anacutsval;
	97	anacutsval=new Float_t*[nvars];
	98	for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
	99	for(Int_t ipt2=0;ipt2<nptbins;ipt2++){
	100	anacutsval[0][ipt2]=0.18;
029ed93a	101	anacutsval[1][ipt2]=0.4;
029ed93a	102	anacutsval[2][ipt2]=0.5;
7ad4b782	103	anacutsval[3][ipt2]=0.;
	104	anacutsval[4][ipt2]=0.;
	105	anacutsval[5][ipt2]=0.01;
029ed93a	106	anacutsval[6][ipt2]=0.06;
	107	anacutsval[7][ipt2]=0.005;
	108	anacutsval[8][ipt2]=0.7;
7ad4b782	109	anacutsval[9][ipt2]=0.;
	110	anacutsval[10][ipt2]=0.;
	111	anacutsval[11][ipt2]=0.05;
a7b533b7	112	anacutsval[12][ipt2]=0.4;
7ad4b782	113	}
7ad4b782	114
a7b533b7	115	// pt kaon
	116	anacutsval[1][0]=0.5;
	117	anacutsval[1][1]=0.85;
	118	anacutsval[1][2]=0.9;
	119	anacutsval[1][3]=0.4;
	120	//pt proton
	121	anacutsval[2][0]=0.5;
	122	anacutsval[2][1]=0.6;
	123	anacutsval[2][2]=0.9;
	124	anacutsval[2][3]=0.9;
	125
	126	//pt pion
	127	anacutsval[12][0]=0.475;
	128	anacutsval[12][1]=0.75;
	129	anacutsval[12][2]=0.75;
	130	anacutsval[12][3]=0.7;
	131
	132	anacutsval[5][0]=0.02;
	133	anacutsval[5][1]=0.025;
	134	anacutsval[5][2]=0.02;
	135	anacutsval[5][3]=0.01;
	136
	137	anacutsval[7][0]=0.00625;
	138	anacutsval[7][1]=0.0125;
	139	anacutsval[7][2]=0.005;
	140	anacutsval[7][3]=0.007;
	141
	142	anacutsval[9][0]=0.5;
	143	anacutsval[9][1]=0.2;
	144	anacutsval[9][2]=0.6;
	145	anacutsval[9][3]=0.;
	146
	147	anacutsval[10][0]=0.00125;
7ad4b782	148
	149	RDHFLctopKpiAn->SetPtBins(nptbins+1,ptbins);
	150	RDHFLctopKpiAn->SetCuts(nvars,nptbins,anacutsval);
	151
	152	// RDHFLc->SetRecoKF(); //set this if you want to recompute the secondary vertex with the KF package
	153	//pid settings
	154	//1. kaon: default one
	155	AliAODPidHF* pidObjK=new AliAODPidHF();
	156	Double_t sigmasK[5]={3.,1.,1.,3.,2.};
	157	pidObjK->SetSigma(sigmasK);
	158	pidObjK->SetAsym(kTRUE);
	159	pidObjK->SetMatch(1);
	160	pidObjK->SetTPC(kTRUE);
	161	pidObjK->SetTOF(kTRUE);
	162	pidObjK->SetITS(kTRUE);
	163	Double_t plimK[2]={0.5,0.8};
	164	pidObjK->SetPLimit(plimK,2);
a7b533b7	165	pidObjK->SetTOFdecide(kTRUE);
7ad4b782	166
	167	RDHFLctopKpiProd->SetPidHF(pidObjK);
	168	RDHFLctopKpiAn->SetPidHF(pidObjK);
	169
	170	//2. pion
	171	AliAODPidHF* pidObjpi=new AliAODPidHF();
	172	pidObjpi->SetTPC(kTRUE);
	173	Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
	174	pidObjpi->SetSigma(sigmaspi);
029ed93a	175	// pidObjpi->SetTOFdecide(kTRUE);
7ad4b782	176
	177	RDHFLctopKpiProd->SetPidpion(pidObjpi);
	178	RDHFLctopKpiAn->SetPidpion(pidObjpi);
	179
	180	// 3. proton
	181	AliAODPidHF* pidObjp=new AliAODPidHF();
	182	Double_t sigmasp[5]={3.,1.,1.,3.,2.};
	183	pidObjp->SetSigma(sigmasp);
	184	pidObjp->SetAsym(kTRUE);
	185	pidObjp->SetMatch(1);
	186	pidObjp->SetTPC(kTRUE);
	187	pidObjp->SetTOF(kTRUE);
	188	pidObjp->SetITS(kTRUE);
	189	Double_t plimp[2]={1.,2.};
	190	pidObjp->SetPLimit(plimp,2);
a7b533b7	191	pidObjp->SetTOFdecide(kTRUE);
7ad4b782	192
	193	RDHFLctopKpiProd->SetPidprot(pidObjp);
	194	RDHFLctopKpiAn->SetPidprot(pidObjp);
	195
	196	Bool_t pidflag=kTRUE;
	197	RDHFLctopKpiAn->SetUsePID(pidflag);
	198	RDHFLctopKpiProd->SetUsePID(pidflag);
	199	if(pidflag) cout<<"PID is used"<<endl;
	200	else cout<<"PID is not used"<<endl;
	201
	202	cout<<"This is the object I'm going to save:"<<endl;
	203	RDHFLctopKpiProd->PrintAll();
	204	RDHFLctopKpiAn->PrintAll();
	205	TFile* fout=new TFile("cuts4LctopKpi.root","RECREATE");
	206	fout->cd();
	207	RDHFLctopKpiProd->Write();
	208	RDHFLctopKpiAn->Write();
	209	fout->Close();
	210	delete fout;
	211	delete pidObjp;
	212	delete pidObjpi;
	213	delete pidObjK;
	214	delete RDHFLctopKpiProd;
	215	delete RDHFLctopKpiAn;
	216
	217	}
5fd38310	218
	219
	220	//macro to make a .root file (for significance maximization) which contains an AliRDHFCutsLctopKpi with loose set of cuts and TParameter with the tighest value of these cuts
	221
	222	void makeInputAliAnalysisTaskSESignificanceMaximization(){
	223
	224	AliRDHFCutsLctopKpi* RDHFLctopKpi=new AliRDHFCutsLctopKpi();
	225	RDHFLctopKpi->SetName("loosercuts");
	226	RDHFLctopKpi->SetTitle("Cuts for significance maximization");
	227
	228	AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
	229	esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
	230	//default
	231	esdTrackCuts->SetRequireTPCRefit(kTRUE);
	232	esdTrackCuts->SetMinNClustersTPC(70);
	233	esdTrackCuts->SetRequireITSRefit(kTRUE);
	234	esdTrackCuts->SetMinNClustersITS(4);
	235
	236	esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
	237	esdTrackCuts->SetMinDCAToVertexXY(0.);
	238	esdTrackCuts->SetEtaRange(-0.8,0.8);
	239	esdTrackCuts->SetPtRange(0.3,1.e10);
	240
	241	RDHFLctopKpi->AddTrackCuts(esdTrackCuts);
	242
a7b533b7	243	const Int_t nvars=13;
5fd38310	244
	245	const Int_t nptbins=4; //change this when adding pt bins!
	246	Float_t ptbins[nptbins+1];
	247	ptbins[0]=0.;
	248	ptbins[1]=2.;
	249	ptbins[2]=3.;
	250	ptbins[3]=4.;
a7b533b7	251	ptbins[4]=12.;
5fd38310	252
	253	RDHFLctopKpi->SetPtBins(nptbins+1,ptbins);
	254
	255	Float_t** rdcutsvalmine;
	256	rdcutsvalmine=new Float_t*[nvars];
	257	for(Int_t iv=0;iv<nvars;iv++){
	258	rdcutsvalmine[iv]=new Float_t[nptbins];
	259	}
	260
	261	//setting my cut values
	262	// inv. mass [GeV]
	263	// pTK [GeV/c]
a7b533b7	264	// pTP [GeV/c]
5fd38310	265	// d0K [cm] lower limit!
	266	// d0Pi [cm] lower limit!
	267	// dist12 (cm)
	268	// sigmavert (cm)
	269	// dist prim-sec (cm)
	270	// pM=Max{pT1,pT2,pT3} (GeV/c)
	271	// cosThetaPoint
	272	// Sum d0^2 (cm^2)
	273	// dca cut (cm)
a7b533b7	274	// pt pion
5fd38310	275	Float_t cutsMatrixLctopKpiStand[nptbins][nvars]=
a7b533b7	276	{{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4},
	277	{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4},
	278	{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4},
	279	{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4}};
5fd38310	280
	281	//CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
	282	Float_t *cutsMatrixTransposeStand=new Float_t[nvars];
	283	for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
	284	for (Int_t ibin=0;ibin<nptbins;ibin++){
	285	for (Int_t ivar = 0; ivar<nvars; ivar++){
	286	cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixLctopKpiStand[ibin][ivar];
	287	}
	288	}
	289	RDHFLctopKpi->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
	290
	291
	292	Int_t nvarsforopt=RDHFLctopKpi->GetNVarsForOpt();
	293	Int_t dim=4; //set this!!
	294	Bool_t *boolforopt;
	295	boolforopt=new Bool_t[nvars];
	296	if(dim>nvarsforopt){
	297	cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
	298	return;
	299	} else {
	300	if(dim==nvarsforopt){
	301	boolforopt=RDHFLctopKpi->GetVarsForOpt();
	302	}else{
	303	TString *names;
	304	names=new TString[nvars];
	305	TString answer="";
	306	Int_t checktrue=0;
	307	names=RDHFLctopKpi->GetVarNames();
	308	for(Int_t i=0;i<nvars;i++){
	309	cout<<names[i]<<" for opt? (y/n)"<<endl;
	310	cin>>answer;
	311	if(answer=="y") {
	312	boolforopt[i]=kTRUE;
	313	checktrue++;
	314	}
	315	else boolforopt[i]=kFALSE;
	316	}
	317	if (checktrue!=dim) {
	318	cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
	319	return;
	320	}
	321	RDHFLctopKpi->SetVarsForOpt(dim,boolforopt);
	322	}
	323	}
	324
	325
	326	Float_t tighterval[dim][nptbins];
	327	// 5 dist12 (cm) <---
	328	// 7 dist prim-sec (cm) <---
	329	// 9 cosThetaPoint <---
	330	// 10 Sum d0^2 (cm^2) <---
	331	// 11 dca cut (cm)
	332	// {0.18,0.4,0.5,0.,0.,0.01 <--- ,0.06,0.005 <--- ,0.7,0. <--- ,0. <--- ,0.05 }
	333
	334	//number of steps for each variable is set in the AddTask arguments (default=8)
	335	//set this!!
	336	for(Int_t ipt=0;ipt<nptbins;ipt++){
	337	tighterval[0][ipt]=0.03;
	338	tighterval[1][ipt]=0.02;
	339	tighterval[2][ipt]=1.;
	340	tighterval[3][ipt]=0.01;
	341	}
	342
	343
344
345	TString name="";
346	Int_t arrdim=dim*nptbins;
347	cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
348	TClonesArray max("TParameter<float>",arrdim);
349	for(Int_t ival=0;ival<dim;ival++){
350	for(Int_t jpt=0;jpt<nptbins;jpt++){
351	name=Form("par%dptbin%d",ival,jpt);
352	cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
353	new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
354	}
355	}
356
357	Bool_t flagPID=kFALSE;
358	RDHFLctopKpi->SetUsePID(flagPID);
359
360	RDHFLctopKpi->PrintAll();
361	printf("Use PID? %s\n",flagPID ? "yes" : "no");
362
363	//pid settings
364	//1. kaon: default one
365	AliAODPidHF* pidObjK=new AliAODPidHF();
366	Double_t sigmasK[5]={3.,1.,1.,3.,2.};
367	pidObjK->SetSigma(sigmasK);
368	pidObjK->SetAsym(kTRUE);
369	pidObjK->SetMatch(1);
370	pidObjK->SetTPC(kTRUE);
371	pidObjK->SetTOF(kTRUE);
372	pidObjK->SetITS(kTRUE);
373	Double_t plimK[2]={0.5,0.8};
374	pidObjK->SetPLimit(plimK,2);
a7b533b7	375	pidObjK->SetTOFdecide(kTRUE);
5fd38310	376
	377	RDHFLctopKpi->SetPidHF(pidObjK);
	378
	379	//2. pion
	380	AliAODPidHF* pidObjpi=new AliAODPidHF();
	381	pidObjpi->SetTPC(kTRUE);
	382	Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
	383	pidObjpi->SetSigma(sigmaspi);
a7b533b7	384	pidObjpi->SetTOFdecide(kTRUE);
5fd38310	385
	386	RDHFLctopKpi->SetPidpion(pidObjpi);
	387
	388	// 3. proton
	389	AliAODPidHF* pidObjp=new AliAODPidHF();
	390	Double_t sigmasp[5]={3.,1.,1.,3.,2.};
	391	pidObjp->SetSigma(sigmasp);
	392	pidObjp->SetAsym(kTRUE);
	393	pidObjp->SetMatch(1);
	394	pidObjp->SetTPC(kTRUE);
	395	pidObjp->SetTOF(kTRUE);
	396	pidObjp->SetITS(kTRUE);
	397	Double_t plimp[2]={1.,2.};
	398	pidObjp->SetPLimit(plimp,2);
a7b533b7	399	pidObjp->SetTOFdecide(kTRUE);
5fd38310	400
	401	RDHFLctopKpi->SetPidprot(pidObjp);
	402
	403	//activate pileup rejection (for pp)
	404	//RDHFLctopKpi->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
	405
	406	//Do not recalculate the vertex
	407	RDHFLctopKpi->SetRemoveDaughtersFromPrim(kFALSE); //activate for pp
	408
	409	TString cent="";
	410	//centrality selection (Pb-Pb)
	411	Float_t minc=20,maxc=80;
	412	RDHFLctopKpi->SetMinCentrality(minc);
	413	RDHFLctopKpi->SetMaxCentrality(maxc);
	414	cent=Form("%.0f%.0f",minc,maxc);
	415	// RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
	416	RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
	417
	418	//temporary
	419	// RDHFLctopKpi->SetFixRefs();
	420
	421	TFile* fout=new TFile(Form("cuts4SignifMaxim%s%s%sRecVtx%sPileupRej.root", RDHFLctopKpi->GetUseCentrality()==0 ? "pp" : "PbPb",cent.Data(),RDHFLctopKpi->GetIsPrimaryWithoutDaughters() ? "" : "No",RDHFLctopKpi->GetOptPileUp() ? "" : "No"),"recreate"); //set this!!
	422
	423	fout->cd();
	424	RDHFLctopKpi->Write();
	425	max.Write();
	426	fout->Close();
	427
	428	}
	429