[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=12;

  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]=9999.;
  

  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;
  }

  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.6;
   anacutsval[3][ipt2]=0.;
   anacutsval[4][ipt2]=0.;
   anacutsval[5][ipt2]=0.01;
   anacutsval[6][ipt2]=0.03;
   anacutsval[9][ipt2]=0.;
   anacutsval[10][ipt2]=0.;
   anacutsval[11][ipt2]=0.05;
  }

  anacutsval[2][0]=0.6;
  anacutsval[2][1]=0.8;
  anacutsval[2][2]=1.;
  anacutsval[2][3]=1.2;

  anacutsval[7][0]=0.005;
  anacutsval[7][1]=0.015;
  anacutsval[7][2]=0.018;
  anacutsval[7][3]=0.018;

  anacutsval[8][0]=0.6;
  anacutsval[8][1]=0.8;
  anacutsval[8][2]=1.;
  anacutsval[8][3]=1.2;

  anacutsval[11][0]=0.04;
  anacutsval[11][1]=0.04;
  anacutsval[11][2]=0.03;
  anacutsval[11][3]=0.03;


  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);
  
  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);

  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);

  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=12;

  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]=9999.;

  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]
  // 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)
  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.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05},
   {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05},
   {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05}};

  //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);
  
  RDHFLctopKpi->SetPidHF(pidObjK);

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

  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);

  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
	47	const Int_t nvars=12;
	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.;
	56	ptbins[4]=9999.;
	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;
	90	}
	91
	92	RDHFLctopKpiProd->SetPtBins(nptbins+1,ptbins);
	93	RDHFLctopKpiProd->SetCuts(nvars,nptbins,prodcutsval);
	94
	95	Float_t** anacutsval;
	96	anacutsval=new Float_t*[nvars];
	97	for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
	98	for(Int_t ipt2=0;ipt2<nptbins;ipt2++){
	99	anacutsval[0][ipt2]=0.18;
	100	anacutsval[1][ipt2]=0.6;
	101	anacutsval[3][ipt2]=0.;
	102	anacutsval[4][ipt2]=0.;
	103	anacutsval[5][ipt2]=0.01;
	104	anacutsval[6][ipt2]=0.03;
	105	anacutsval[9][ipt2]=0.;
	106	anacutsval[10][ipt2]=0.;
	107	anacutsval[11][ipt2]=0.05;
	108	}
	109
	110	anacutsval[2][0]=0.6;
	111	anacutsval[2][1]=0.8;
	112	anacutsval[2][2]=1.;
	113	anacutsval[2][3]=1.2;
	114
	115	anacutsval[7][0]=0.005;
	116	anacutsval[7][1]=0.015;
	117	anacutsval[7][2]=0.018;
	118	anacutsval[7][3]=0.018;
	119
	120	anacutsval[8][0]=0.6;
	121	anacutsval[8][1]=0.8;
	122	anacutsval[8][2]=1.;
	123	anacutsval[8][3]=1.2;
	124
	125	anacutsval[11][0]=0.04;
	126	anacutsval[11][1]=0.04;
	127	anacutsval[11][2]=0.03;
	128	anacutsval[11][3]=0.03;
	129
	130
	131	RDHFLctopKpiAn->SetPtBins(nptbins+1,ptbins);
	132	RDHFLctopKpiAn->SetCuts(nvars,nptbins,anacutsval);
	133
	134	// RDHFLc->SetRecoKF(); //set this if you want to recompute the secondary vertex with the KF package
	135	//pid settings
	136	//1. kaon: default one
	137	AliAODPidHF* pidObjK=new AliAODPidHF();
	138	Double_t sigmasK[5]={3.,1.,1.,3.,2.};
	139	pidObjK->SetSigma(sigmasK);
	140	pidObjK->SetAsym(kTRUE);
	141	pidObjK->SetMatch(1);
	142	pidObjK->SetTPC(kTRUE);
	143	pidObjK->SetTOF(kTRUE);
	144	pidObjK->SetITS(kTRUE);
	145	Double_t plimK[2]={0.5,0.8};
	146	pidObjK->SetPLimit(plimK,2);
	147
	148	RDHFLctopKpiProd->SetPidHF(pidObjK);
	149	RDHFLctopKpiAn->SetPidHF(pidObjK);
	150
151	//2. pion
152	AliAODPidHF* pidObjpi=new AliAODPidHF();
153	pidObjpi->SetTPC(kTRUE);
154	Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
155	pidObjpi->SetSigma(sigmaspi);
156
157	RDHFLctopKpiProd->SetPidpion(pidObjpi);
158	RDHFLctopKpiAn->SetPidpion(pidObjpi);
159
160	// 3. proton
161	AliAODPidHF* pidObjp=new AliAODPidHF();
162	Double_t sigmasp[5]={3.,1.,1.,3.,2.};
163	pidObjp->SetSigma(sigmasp);
164	pidObjp->SetAsym(kTRUE);
165	pidObjp->SetMatch(1);
166	pidObjp->SetTPC(kTRUE);
167	pidObjp->SetTOF(kTRUE);
168	pidObjp->SetITS(kTRUE);
169	Double_t plimp[2]={1.,2.};
170	pidObjp->SetPLimit(plimp,2);
171
172	RDHFLctopKpiProd->SetPidprot(pidObjp);
173	RDHFLctopKpiAn->SetPidprot(pidObjp);
174
175	Bool_t pidflag=kTRUE;
176	RDHFLctopKpiAn->SetUsePID(pidflag);
177	RDHFLctopKpiProd->SetUsePID(pidflag);
178	if(pidflag) cout<<"PID is used"<<endl;
179	else cout<<"PID is not used"<<endl;
180
181	cout<<"This is the object I'm going to save:"<<endl;
182	RDHFLctopKpiProd->PrintAll();
183	RDHFLctopKpiAn->PrintAll();
184	TFile* fout=new TFile("cuts4LctopKpi.root","RECREATE");
185	fout->cd();
186	RDHFLctopKpiProd->Write();
187	RDHFLctopKpiAn->Write();
188	fout->Close();
189	delete fout;
190	delete pidObjp;
191	delete pidObjpi;
192	delete pidObjK;
193	delete RDHFLctopKpiProd;
194	delete RDHFLctopKpiAn;
195
196	}
5fd38310	197
	198
	199	//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
	200
	201	void makeInputAliAnalysisTaskSESignificanceMaximization(){
	202
	203	AliRDHFCutsLctopKpi* RDHFLctopKpi=new AliRDHFCutsLctopKpi();
	204	RDHFLctopKpi->SetName("loosercuts");
	205	RDHFLctopKpi->SetTitle("Cuts for significance maximization");
	206
	207	AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
	208	esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
	209	//default
	210	esdTrackCuts->SetRequireTPCRefit(kTRUE);
	211	esdTrackCuts->SetMinNClustersTPC(70);
	212	esdTrackCuts->SetRequireITSRefit(kTRUE);
	213	esdTrackCuts->SetMinNClustersITS(4);
	214
	215	esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
	216	esdTrackCuts->SetMinDCAToVertexXY(0.);
	217	esdTrackCuts->SetEtaRange(-0.8,0.8);
	218	esdTrackCuts->SetPtRange(0.3,1.e10);
	219
	220	RDHFLctopKpi->AddTrackCuts(esdTrackCuts);
	221
	222	const Int_t nvars=12;
	223
	224	const Int_t nptbins=4; //change this when adding pt bins!
	225	Float_t ptbins[nptbins+1];
	226	ptbins[0]=0.;
	227	ptbins[1]=2.;
	228	ptbins[2]=3.;
	229	ptbins[3]=4.;
	230	ptbins[4]=9999.;
	231
	232	RDHFLctopKpi->SetPtBins(nptbins+1,ptbins);
	233
	234	Float_t** rdcutsvalmine;
	235	rdcutsvalmine=new Float_t*[nvars];
	236	for(Int_t iv=0;iv<nvars;iv++){
	237	rdcutsvalmine[iv]=new Float_t[nptbins];
	238	}
	239
	240	//setting my cut values
	241	// inv. mass [GeV]
	242	// pTK [GeV/c]
	243	// pTPi [GeV/c]
	244	// d0K [cm] lower limit!
	245	// d0Pi [cm] lower limit!
	246	// dist12 (cm)
	247	// sigmavert (cm)
	248	// dist prim-sec (cm)
	249	// pM=Max{pT1,pT2,pT3} (GeV/c)
	250	// cosThetaPoint
	251	// Sum d0^2 (cm^2)
	252	// dca cut (cm)
	253	Float_t cutsMatrixLctopKpiStand[nptbins][nvars]=
	254	{{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05},
	255	{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05},
	256	{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05},
	257	{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05}};
	258
	259	//CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
	260	Float_t *cutsMatrixTransposeStand=new Float_t[nvars];
261	for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
262	for (Int_t ibin=0;ibin<nptbins;ibin++){
263	for (Int_t ivar = 0; ivar<nvars; ivar++){
264	cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixLctopKpiStand[ibin][ivar];
265	}
266	}
267	RDHFLctopKpi->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
268
269
270	Int_t nvarsforopt=RDHFLctopKpi->GetNVarsForOpt();
271	Int_t dim=4; //set this!!
272	Bool_t *boolforopt;
273	boolforopt=new Bool_t[nvars];
274	if(dim>nvarsforopt){
275	cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
276	return;
277	} else {
278	if(dim==nvarsforopt){
279	boolforopt=RDHFLctopKpi->GetVarsForOpt();
280	}else{
281	TString *names;
282	names=new TString[nvars];
283	TString answer="";
284	Int_t checktrue=0;
285	names=RDHFLctopKpi->GetVarNames();
286	for(Int_t i=0;i<nvars;i++){
287	cout<<names[i]<<" for opt? (y/n)"<<endl;
288	cin>>answer;
289	if(answer=="y") {
290	boolforopt[i]=kTRUE;
291	checktrue++;
292	}
293	else boolforopt[i]=kFALSE;
294	}
295	if (checktrue!=dim) {
296	cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
297	return;
298	}
299	RDHFLctopKpi->SetVarsForOpt(dim,boolforopt);
300	}
301	}
302
303
304	Float_t tighterval[dim][nptbins];
305	// 5 dist12 (cm) <---
306	// 7 dist prim-sec (cm) <---
307	// 9 cosThetaPoint <---
308	// 10 Sum d0^2 (cm^2) <---
309	// 11 dca cut (cm)
310	// {0.18,0.4,0.5,0.,0.,0.01 <--- ,0.06,0.005 <--- ,0.7,0. <--- ,0. <--- ,0.05 }
311
312	//number of steps for each variable is set in the AddTask arguments (default=8)
313	//set this!!
314	for(Int_t ipt=0;ipt<nptbins;ipt++){
315	tighterval[0][ipt]=0.03;
316	tighterval[1][ipt]=0.02;
317	tighterval[2][ipt]=1.;
318	tighterval[3][ipt]=0.01;
319	}
320
321
322
323	TString name="";
324	Int_t arrdim=dim*nptbins;
325	cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
326	TClonesArray max("TParameter<float>",arrdim);
327	for(Int_t ival=0;ival<dim;ival++){
328	for(Int_t jpt=0;jpt<nptbins;jpt++){
329	name=Form("par%dptbin%d",ival,jpt);
330	cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
331	new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
332	}
333	}
334
335	Bool_t flagPID=kFALSE;
336	RDHFLctopKpi->SetUsePID(flagPID);
337
338	RDHFLctopKpi->PrintAll();
339	printf("Use PID? %s\n",flagPID ? "yes" : "no");
340
341	//pid settings
342	//1. kaon: default one
343	AliAODPidHF* pidObjK=new AliAODPidHF();
344	Double_t sigmasK[5]={3.,1.,1.,3.,2.};
345	pidObjK->SetSigma(sigmasK);
346	pidObjK->SetAsym(kTRUE);
347	pidObjK->SetMatch(1);
348	pidObjK->SetTPC(kTRUE);
349	pidObjK->SetTOF(kTRUE);
350	pidObjK->SetITS(kTRUE);
351	Double_t plimK[2]={0.5,0.8};
352	pidObjK->SetPLimit(plimK,2);
353
354	RDHFLctopKpi->SetPidHF(pidObjK);
355
356	//2. pion
357	AliAODPidHF* pidObjpi=new AliAODPidHF();
358	pidObjpi->SetTPC(kTRUE);
359	Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
360	pidObjpi->SetSigma(sigmaspi);
361
362	RDHFLctopKpi->SetPidpion(pidObjpi);
363
364	// 3. proton
365	AliAODPidHF* pidObjp=new AliAODPidHF();
366	Double_t sigmasp[5]={3.,1.,1.,3.,2.};
367	pidObjp->SetSigma(sigmasp);
368	pidObjp->SetAsym(kTRUE);
369	pidObjp->SetMatch(1);
370	pidObjp->SetTPC(kTRUE);
371	pidObjp->SetTOF(kTRUE);
372	pidObjp->SetITS(kTRUE);
373	Double_t plimp[2]={1.,2.};
374	pidObjp->SetPLimit(plimp,2);
375
376	RDHFLctopKpi->SetPidprot(pidObjp);
377
378	//activate pileup rejection (for pp)
379	//RDHFLctopKpi->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
380
381	//Do not recalculate the vertex
382	RDHFLctopKpi->SetRemoveDaughtersFromPrim(kFALSE); //activate for pp
383
384	TString cent="";
385	//centrality selection (Pb-Pb)
386	Float_t minc=20,maxc=80;
387	RDHFLctopKpi->SetMinCentrality(minc);
388	RDHFLctopKpi->SetMaxCentrality(maxc);
389	cent=Form("%.0f%.0f",minc,maxc);
390	// RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
391	RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
392
393	//temporary
394	// RDHFLctopKpi->SetFixRefs();
395
396	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!!
397
398	fout->cd();
399	RDHFLctopKpi->Write();
400	max.Write();
401	fout->Close();
402
403	}
404