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

#include <Riostream.h>
#include <TFile.h>
//#include <AliRDHFCutsDstoKKpi.h>
#include <TClonesArray.h>
#include <TParameter.h>

//macro to make a .root file which contains an AliRDHFCutsDstoKKpi with loose set of cuts (for significance maximization) and TParameter with the tighest value of these cuts
//Needed for AliAnalysisTaskSEDs, AliCFTaskVertexingHF3Prong, AliAnalysisTaskSESignificance

//Use:
//Set hard coded commented with //set this!!

//.L makeTFile4CutsDstoKKpi.C
// makeInputAliAnalysisTaskSEDs()
// makeInputAliAnalysisTaskSESignificanceMaximization()

//similar macros for the other D mesons

//Author: Chiara Bianchin, cbianchi@pd.infn.it
//        Giacomo Ortona, ortona@to.infn.it
//        Renu Bala [Dplus Analysis and CF]

//Modified for Ds meson: G.M. Innocenti, innocent@to.infn.it

void makeInputAliAnalysisTaskSEDsPP(){

 //  gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/PWG3 -I$ALICE_ROOT/PWG3/vertexingHF -I$ALICE_ROOT/PWG3/vertexingH/macros -g"); 

    
    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);
    
    
    const Int_t nptbins=4;
    Float_t* ptbins;
    ptbins=new Float_t[nptbins+1];
    ptbins[0]=2.;
    ptbins[1]=4.;
    ptbins[2]=6.;
    ptbins[3]=8.;
    ptbins[4]=12.;
    
    
    const Int_t nvars=20;
    
    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 ipt=0;ipt<nptbins;ipt++){
      
      anacutsval[0][ipt]=0.35;
      anacutsval[1][ipt]=0.3;
      anacutsval[2][ipt]=0.3;
      anacutsval[3][ipt]=0.;
      anacutsval[4][ipt]=0.;
      anacutsval[5][ipt]=0.005;
      anacutsval[6][ipt]=0.06;
      anacutsval[7][ipt]=0.0;
      anacutsval[8][ipt]=0.;
      anacutsval[9][ipt]=0.9;
      anacutsval[10][ipt]=0.;
      anacutsval[11][ipt]=1000.0;
      anacutsval[12][ipt]=0.015;
      anacutsval[13][ipt]=0.1;
      anacutsval[14][ipt]=0.;
      anacutsval[15][ipt]=1.;
      anacutsval[16][ipt]=0.;
      anacutsval[17][ipt]=0.;
      anacutsval[18][ipt]=0.;
      anacutsval[19][ipt]=-1.;
      
      
    }
    
    
    /*    
   
    Cut list                                                rejection condition
    0           "inv. mass [GeV]",                          invmassDS-massDspdg>fCutsRD
    1			"pTK [GeV/c]",                              pTK<fCutsRd
    2			"pTPi [GeV/c]",                             pTPi<fCutsRd
    3			"d0K [cm]",                                 d0K<fCutsRd
    4			"d0Pi [cm]",                                d0Pi<fCutsRd
    5			"dist12 [cm]",                              dist12<fCutsRd
    6			"sigmavert [cm]",                           sigmavert>fCutsRd
    7			"decLen [cm]",                              decLen<fCutsRD
    8			"ptMax [GeV/c]",                            ptMax<fCutsRD
    9			"cosThetaPoint",                            CosThetaPoint<fCutsRD
    10			"Sum d0^2 (cm^2)",                          sumd0<fCutsRD
    11			"dca [cm]",                                 dca(i)>fCutsRD
    12			"inv. mass (Mphi-MKK) [GeV]",               invmass-pdg>fCutsRD
    13			"inv. mass (MKo*-MKpi) [GeV]"};             invmass-pdg>fCutsRD
    14    		"Abs(CosineKpiPhiRFrame)^3",
	15  		"CosPiDsLabFrame"};
    */
 
    
    AliRDHFCutsDstoKKpi* analysiscuts=new AliRDHFCutsDstoKKpi();
    analysiscuts->SetName("AnalysisCuts");
    analysiscuts->SetTitle("Cuts for Ds Analysis and CF");
    analysiscuts->SetPtBins(nptbins+1,ptbins);
    analysiscuts->SetCuts(nvars,nptbins,anacutsval);
    analysiscuts->AddTrackCuts(esdTrackCuts);
    analysiscuts->SetUsePID(kTRUE);
    analysiscuts->SetPidOption(1);
    analysiscuts->SetOptPileup(kTRUE);
    analysiscuts->SetRemoveDaughtersFromPrim(kTRUE);
    
    // To be set only in case of strong pid
    //analysiscuts->SetMaxPtStrongPid(9999.);
    //analysiscuts->SetMaxPStrongPidK(1.5);
    //analysiscuts->SetMaxPStrongPidpi(1.5);
    cout<<"This is the odject I'm going to save:"<<nptbins<<endl;
    
    analysiscuts->PrintAll();
    TFile* fout=new TFile("DstoKKpiCuts.root","recreate");   
    fout->cd();
    analysiscuts->Write();
    fout->Close();
    
    
}


void makeInputAliAnalysisTaskSEDsPbPb(){

 //  gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/PWG3 -I$ALICE_ROOT/PWG3/vertexingHF -I$ALICE_ROOT/PWG3/vertexingH/macros -g"); 

    
    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.7,1.e10);
    
    
    const Int_t nptbins=4;
    Float_t* ptbins;
    ptbins=new Float_t[nptbins+1];
    ptbins[0]=2.;
    ptbins[1]=4.;
    ptbins[2]=6.;
    ptbins[3]=8.;
    ptbins[4]=12.;
    
    
    const Int_t nvars=20;
    
    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 ipt=0;ipt<nptbins;ipt++){
      
      anacutsval[0][ipt]=0.35;
      anacutsval[1][ipt]=0.3;
      anacutsval[2][ipt]=0.3;
      anacutsval[3][ipt]=0.;
      anacutsval[4][ipt]=0.;
      anacutsval[5][ipt]=0.005;
      anacutsval[6][ipt]=0.06;
      anacutsval[7][ipt]=0.0;
      anacutsval[8][ipt]=0.;
      anacutsval[9][ipt]=0.7;
      anacutsval[10][ipt]=0.;
      anacutsval[11][ipt]=1000.0;
      anacutsval[12][ipt]=0.1;
      anacutsval[13][ipt]=0.1;
      anacutsval[14][ipt]=0.;
      anacutsval[15][ipt]=1.;
      anacutsval[16][ipt]=0.;
      anacutsval[17][ipt]=0.;
      anacutsval[18][ipt]=0.;
      anacutsval[19][ipt]=-1.;
      
   
    }
    
    
    /*    
   
    Cut list                                                rejection condition
    0           "inv. mass [GeV]",                          invmassDS-massDspdg>fCutsRD
    1			"pTK [GeV/c]",                              pTK<fCutsRd
    2			"pTPi [GeV/c]",                             pTPi<fCutsRd
    3			"d0K [cm]",                                 d0K<fCutsRd
    4			"d0Pi [cm]",                                d0Pi<fCutsRd
    5			"dist12 [cm]",                              dist12<fCutsRd
    6			"sigmavert [cm]",                           sigmavert>fCutsRd
    7			"decLen [cm]",                              decLen<fCutsRD
    8			"ptMax [GeV/c]",                            ptMax<fCutsRD
    9			"cosThetaPoint",                            CosThetaPoint<fCutsRD
    10			"Sum d0^2 (cm^2)",                          sumd0<fCutsRD
    11			"dca [cm]",                                 dca(i)>fCutsRD
    12			"inv. mass (Mphi-MKK) [GeV]",               invmass-pdg>fCutsRD
    13			"inv. mass (MKo*-MKpi) [GeV]"};             invmass-pdg>fCutsRD
    14    		"Abs(CosineKpiPhiRFrame)^3",
	15  		"CosPiDsLabFrame"};
    */
    
    AliRDHFCutsDstoKKpi* analysiscuts=new AliRDHFCutsDstoKKpi();
    analysiscuts->SetName("AnalysisCuts");
    analysiscuts->SetTitle("Cuts for Ds Analysis and CF");
    analysiscuts->SetPtBins(nptbins+1,ptbins);
    analysiscuts->SetCuts(nvars,nptbins,anacutsval);
    analysiscuts->AddTrackCuts(esdTrackCuts);
    
    TString cent="";
    Float_t mincen=20.;
    Float_t maxcen=40.;
    
    analysiscuts->SetUsePID(kTRUE);
    analysiscuts->SetPidOption(0);
    //analysiscuts->SetUseImpParProdCorrCut(kFALSE);
    analysiscuts->SetOptPileup(kFALSE);
    //analysiscuts->SetUseAOD049(kTRUE);
    analysiscuts->SetMinCentrality(mincen);
    analysiscuts->SetMaxCentrality(maxcen);
    cent=Form("%.0f%.0f",mincen,maxcen);
    analysiscuts->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid              
    analysiscuts->SetMinPtCandidate(3.);
    analysiscuts->SetMaxPtCandidate(1000.);
    analysiscuts->SetRemoveDaughtersFromPrim(kFALSE);
    // To be set only in case of strong pid
    //analysiscuts->SetMaxPtStrongPid(9999.);
    //analysiscuts->SetMaxPStrongPidK(1.5);
    //analysiscuts->SetMaxPStrongPidpi(1.5);
    
    cout<<"This is the odject I'm going to save:"<<nptbins<<endl;
      
    analysiscuts->PrintAll();
    TFile* fout=new TFile("DstoKKpiCuts.root","recreate");   
    fout->cd();
    analysiscuts->Write();
    fout->Close();
    
    
}


void makeInputAliAnalysisTaskSESignificanceMaximization(){
  
  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);
  
  
  AliRDHFCutsDstoKKpi* RDHFDstoKKpi=new AliRDHFCutsDstoKKpi();
  RDHFDstoKKpi->SetName("loosercuts");
  RDHFDstoKKpi->SetTitle("Cuts for significance maximization");
  RDHFDstoKKpi->AddTrackCuts(esdTrackCuts);
  RDHFDstoKKpi->SetUsePID(kTRUE);
  RDHFDstoKKpi->SetPidOption(1);
  RDHFDstoKKpi->SetOptPileup(kTRUE);
  RDHFDstoKKpi->SetRemoveDaughtersFromPrim(kTRUE);

  const Int_t nvars=20;

  const Int_t nptbins=4;
  Float_t* ptbins;
  ptbins=new Float_t[nptbins+1];
  ptbins[0]=2.;
  ptbins[1]=4.;
  ptbins[2]=6.;
  ptbins[3]=8.;
  ptbins[4]=12.;
  
  RDHFDstoKKpi->SetPtBins(nptbins+1,ptbins);
    
  //setting my cut values
   
  Float_t** prodcutsval;
  prodcutsval=new Float_t*[nvars];
  for(Int_t ic=0;ic<nvars;ic++){prodcutsval[ic]=new Float_t[nptbins];}  
  for(Int_t ipt=0;ipt<nptbins;ipt++){
    
      prodcutsval[0][ipt]=0.35;
      prodcutsval[1][ipt]=0.3;
      prodcutsval[2][ipt]=0.3;
      prodcutsval[3][ipt]=0.;
      prodcutsval[4][ipt]=0.;
      prodcutsval[5][ipt]=0.005;
      prodcutsval[6][ipt]=0.06;
      prodcutsval[7][ipt]=0.0;
      prodcutsval[8][ipt]=0.;
      prodcutsval[9][ipt]=0.9;
      prodcutsval[10][ipt]=0.;
      prodcutsval[11][ipt]=1000.0;
      prodcutsval[12][ipt]=0.015;
      prodcutsval[13][ipt]=0.1;
      prodcutsval[14][ipt]=0.;
      prodcutsval[15][ipt]=1.;
      prodcutsval[16][ipt]=0.;
      prodcutsval[17][ipt]=0.;
      prodcutsval[18][ipt]=0.;
      prodcutsval[19][ipt]=-1.;
      
  }

  RDHFDstoKKpi->SetCuts(nvars,nptbins,prodcutsval);

  Int_t nvarsforopt=RDHFDstoKKpi->GetNVarsForOpt();
  //Int_t nvarsforopt=2;
  const 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=RDHFDstoKKpi->GetVarsForOpt();
    }else{
      TString *names;
      names=new TString[nvars];
      TString answer="";
      Int_t checktrue=0;
      names=RDHFDstoKKpi->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;
      }
      RDHFDstoKKpi->SetVarsForOpt(dim,boolforopt);
    }
  }


  Float_t tighterval[dim][nptbins];
 
  //number of steps for each variable is 4 now
  //set this!!
  
  /*    
   Cut list                                                rejection condition

   0           "inv. mass [GeV]",                          invmassDS-massDspdg>fCutsRD
   1			"pTK [GeV/c]",                              pTK<fCutsRd
   2			"pTPi [GeV/c]",                             pTPi<fCutsRd
   3			"d0K [cm]",                                 d0K<fCutsRd
   4			"d0Pi [cm]",                                d0Pi<fCutsRd
   5			"dist12 [cm]",                              dist12<fCutsRd
   6			"sigmavert [cm]",                           sigmavert>fCutsRd
   7			"decLen [cm]",                              decLen<fCutsRD
   8			"ptMax [GeV/c]",                            ptMax<fCutsRD
   9			"cosThetaPoint",                            CosThetaPoint<fCutsRD
   10			"Sum d0^2 (cm^2)",                          sumd0<fCutsRD
   11			"dca [cm]",                                 dca(i)>fCutsRD
   12			"inv. mass (Mphi-MKK) [GeV]",               invmass-pdg>fCutsRD
   13			"inv. mass (MKo*-MKpi) [GeV]"};             invmass-pdg>fCutsRD
   14    		"Abs(CosineKpiPhiRFrame)^3",
   15  	    	"CosPiDsLabFrame"}
   */  
  
  //sigmavert [cm]
  
   
  tighterval[0][0]=0.0;
  tighterval[0][1]=0.0;
  tighterval[0][2]=0.0;
  tighterval[0][3]=0.0;
 
  //decay length
  
  tighterval[1][0]=0.05;
  tighterval[1][1]=0.05;
  tighterval[1][2]=0.05;
  tighterval[1][3]=0.05;
 
  //costhetapoint
  
  tighterval[2][0]=1.;
  tighterval[2][1]=1.;
  tighterval[2][2]=1.;
  tighterval[2][3]=1.;
  
  
  //inv mass phi meson
  
  tighterval[3][0]=0.;
  tighterval[3][1]=0.;
  tighterval[3][2]=0.;
  tighterval[3][3]=0.;

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

  TFile* fout=new TFile("cuts4SignifMaximDs.root","recreate");   //set this!! 
  fout->cd();
  RDHFDstoKKpi->Write();
  max.Write();
  fout->Close();

}
Commit	Line	Data
1f796bc2	1	#include <Riostream.h>
	2	#include <TFile.h>
	3	//#include <AliRDHFCutsDstoKKpi.h>
	4	#include <TClonesArray.h>
	5	#include <TParameter.h>
	6
	7	//macro to make a .root file which contains an AliRDHFCutsDstoKKpi with loose set of cuts (for significance maximization) and TParameter with the tighest value of these cuts
	8	//Needed for AliAnalysisTaskSEDs, AliCFTaskVertexingHF3Prong, AliAnalysisTaskSESignificance
	9
	10	//Use:
	11	//Set hard coded commented with //set this!!
	12
	13	//.L makeTFile4CutsDstoKKpi.C
	14	// makeInputAliAnalysisTaskSEDs()
	15	// makeInputAliAnalysisTaskSESignificanceMaximization()
	16
	17	//similar macros for the other D mesons
	18
	19	//Author: Chiara Bianchin, cbianchi@pd.infn.it
	20	// Giacomo Ortona, ortona@to.infn.it
	21	// Renu Bala [Dplus Analysis and CF]
	22
	23	//Modified for Ds meson: G.M. Innocenti, innocent@to.infn.it
	24
	25	void makeInputAliAnalysisTaskSEDsPP(){
	26
	27	// gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/PWG3 -I$ALICE_ROOT/PWG3/vertexingHF -I$ALICE_ROOT/PWG3/vertexingH/macros -g");
	28
	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	const Int_t nptbins=4;
	45	Float_t* ptbins;
	46	ptbins=new Float_t[nptbins+1];
	47	ptbins[0]=2.;
	48	ptbins[1]=4.;
	49	ptbins[2]=6.;
	50	ptbins[3]=8.;
	51	ptbins[4]=12.;
	52
	53
ac8cab2d	54	const Int_t nvars=20;
1f796bc2	55
	56	Float_t** anacutsval;
	57	anacutsval=new Float_t*[nvars];
	58
	59	for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
	60	for(Int_t ipt=0;ipt<nptbins;ipt++){
	61
	62	anacutsval[0][ipt]=0.35;
	63	anacutsval[1][ipt]=0.3;
	64	anacutsval[2][ipt]=0.3;
	65	anacutsval[3][ipt]=0.;
	66	anacutsval[4][ipt]=0.;
	67	anacutsval[5][ipt]=0.005;
	68	anacutsval[6][ipt]=0.06;
	69	anacutsval[7][ipt]=0.0;
	70	anacutsval[8][ipt]=0.;
	71	anacutsval[9][ipt]=0.9;
	72	anacutsval[10][ipt]=0.;
	73	anacutsval[11][ipt]=1000.0;
	74	anacutsval[12][ipt]=0.015;
	75	anacutsval[13][ipt]=0.1;
	76	anacutsval[14][ipt]=0.;
	77	anacutsval[15][ipt]=1.;
ac8cab2d	78	anacutsval[16][ipt]=0.;
	79	anacutsval[17][ipt]=0.;
	80	anacutsval[18][ipt]=0.;
	81	anacutsval[19][ipt]=-1.;
	82
1f796bc2	83
	84
	85	}
	86
	87
	88
	89	/*
	90
	91	Cut list rejection condition
	92	0 "inv. mass [GeV]", invmassDS-massDspdg>fCutsRD
	93	1 "pTK [GeV/c]", pTK<fCutsRd
	94	2 "pTPi [GeV/c]", pTPi<fCutsRd
	95	3 "d0K [cm]", d0K<fCutsRd
	96	4 "d0Pi [cm]", d0Pi<fCutsRd
	97	5 "dist12 [cm]", dist12<fCutsRd
	98	6 "sigmavert [cm]", sigmavert>fCutsRd
	99	7 "decLen [cm]", decLen<fCutsRD
	100	8 "ptMax [GeV/c]", ptMax<fCutsRD
	101	9 "cosThetaPoint", CosThetaPoint<fCutsRD
	102	10 "Sum d0^2 (cm^2)", sumd0<fCutsRD
	103	11 "dca [cm]", dca(i)>fCutsRD
	104	12 "inv. mass (Mphi-MKK) [GeV]", invmass-pdg>fCutsRD
	105	13 "inv. mass (MKo*-MKpi) [GeV]"}; invmass-pdg>fCutsRD
	106	14 "Abs(CosineKpiPhiRFrame)^3",
	107	15 "CosPiDsLabFrame"};
	108	*/
	109
	110
	111	AliRDHFCutsDstoKKpi* analysiscuts=new AliRDHFCutsDstoKKpi();
	112	analysiscuts->SetName("AnalysisCuts");
	113	analysiscuts->SetTitle("Cuts for Ds Analysis and CF");
	114	analysiscuts->SetPtBins(nptbins+1,ptbins);
	115	analysiscuts->SetCuts(nvars,nptbins,anacutsval);
	116	analysiscuts->AddTrackCuts(esdTrackCuts);
	117	analysiscuts->SetUsePID(kTRUE);
	118	analysiscuts->SetPidOption(1);
	119	analysiscuts->SetOptPileup(kTRUE);
	120	analysiscuts->SetRemoveDaughtersFromPrim(kTRUE);
ac8cab2d	121
	122	// To be set only in case of strong pid
	123	//analysiscuts->SetMaxPtStrongPid(9999.);
	124	//analysiscuts->SetMaxPStrongPidK(1.5);
	125	//analysiscuts->SetMaxPStrongPidpi(1.5);
1f796bc2	126	cout<<"This is the odject I'm going to save:"<<nptbins<<endl;
	127
	128	analysiscuts->PrintAll();
	129	TFile* fout=new TFile("DstoKKpiCuts.root","recreate");
	130	fout->cd();
	131	analysiscuts->Write();
	132	fout->Close();
	133
	134
	135	}
	136
	137
	138	void makeInputAliAnalysisTaskSEDsPbPb(){
	139
	140	// gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/PWG3 -I$ALICE_ROOT/PWG3/vertexingHF -I$ALICE_ROOT/PWG3/vertexingH/macros -g");
	141
	142
	143	AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
	144	esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
	145	//default
	146	esdTrackCuts->SetRequireTPCRefit(kTRUE);
	147	esdTrackCuts->SetRequireITSRefit(kTRUE);
	148	//esdTrackCuts->SetMinNClustersITS(4); // default is 5
	149	esdTrackCuts->SetMinNClustersTPC(70);
	150	esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
	151	AliESDtrackCuts::kAny);
	152	// default is kBoth, otherwise kAny
	153	esdTrackCuts->SetMinDCAToVertexXY(0.);
	154	esdTrackCuts->SetPtRange(0.7,1.e10);
	155
	156
	157	const Int_t nptbins=4;
	158	Float_t* ptbins;
	159	ptbins=new Float_t[nptbins+1];
	160	ptbins[0]=2.;
	161	ptbins[1]=4.;
	162	ptbins[2]=6.;
	163	ptbins[3]=8.;
	164	ptbins[4]=12.;
	165
	166
ac8cab2d	167	const Int_t nvars=20;
1f796bc2	168
	169	Float_t** anacutsval;
	170	anacutsval=new Float_t*[nvars];
	171
	172	for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
	173	for(Int_t ipt=0;ipt<nptbins;ipt++){
	174
	175	anacutsval[0][ipt]=0.35;
	176	anacutsval[1][ipt]=0.3;
	177	anacutsval[2][ipt]=0.3;
	178	anacutsval[3][ipt]=0.;
	179	anacutsval[4][ipt]=0.;
	180	anacutsval[5][ipt]=0.005;
	181	anacutsval[6][ipt]=0.06;
	182	anacutsval[7][ipt]=0.0;
	183	anacutsval[8][ipt]=0.;
	184	anacutsval[9][ipt]=0.7;
	185	anacutsval[10][ipt]=0.;
	186	anacutsval[11][ipt]=1000.0;
	187	anacutsval[12][ipt]=0.1;
	188	anacutsval[13][ipt]=0.1;
	189	anacutsval[14][ipt]=0.;
	190	anacutsval[15][ipt]=1.;
ac8cab2d	191	anacutsval[16][ipt]=0.;
	192	anacutsval[17][ipt]=0.;
	193	anacutsval[18][ipt]=0.;
	194	anacutsval[19][ipt]=-1.;
1f796bc2	195
	196
	197	}
	198
	199
	200
	201	/*
	202
	203	Cut list rejection condition
	204	0 "inv. mass [GeV]", invmassDS-massDspdg>fCutsRD
	205	1 "pTK [GeV/c]", pTK<fCutsRd
	206	2 "pTPi [GeV/c]", pTPi<fCutsRd
	207	3 "d0K [cm]", d0K<fCutsRd
	208	4 "d0Pi [cm]", d0Pi<fCutsRd
	209	5 "dist12 [cm]", dist12<fCutsRd
	210	6 "sigmavert [cm]", sigmavert>fCutsRd
	211	7 "decLen [cm]", decLen<fCutsRD
	212	8 "ptMax [GeV/c]", ptMax<fCutsRD
	213	9 "cosThetaPoint", CosThetaPoint<fCutsRD
	214	10 "Sum d0^2 (cm^2)", sumd0<fCutsRD
	215	11 "dca [cm]", dca(i)>fCutsRD
	216	12 "inv. mass (Mphi-MKK) [GeV]", invmass-pdg>fCutsRD
	217	13 "inv. mass (MKo*-MKpi) [GeV]"}; invmass-pdg>fCutsRD
	218	14 "Abs(CosineKpiPhiRFrame)^3",
	219	15 "CosPiDsLabFrame"};
	220	*/
	221
	222	AliRDHFCutsDstoKKpi* analysiscuts=new AliRDHFCutsDstoKKpi();
	223	analysiscuts->SetName("AnalysisCuts");
	224	analysiscuts->SetTitle("Cuts for Ds Analysis and CF");
	225	analysiscuts->SetPtBins(nptbins+1,ptbins);
	226	analysiscuts->SetCuts(nvars,nptbins,anacutsval);
	227	analysiscuts->AddTrackCuts(esdTrackCuts);
	228
	229	TString cent="";
	230	Float_t mincen=20.;
	231	Float_t maxcen=40.;
	232
	233	analysiscuts->SetUsePID(kTRUE);
	234	analysiscuts->SetPidOption(0);
	235	//analysiscuts->SetUseImpParProdCorrCut(kFALSE);
	236	analysiscuts->SetOptPileup(kFALSE);
	237	//analysiscuts->SetUseAOD049(kTRUE);
	238	analysiscuts->SetMinCentrality(mincen);
	239	analysiscuts->SetMaxCentrality(maxcen);
	240	cent=Form("%.0f%.0f",mincen,maxcen);
	241	analysiscuts->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
	242	analysiscuts->SetMinPtCandidate(3.);
	243	analysiscuts->SetMaxPtCandidate(1000.);
	244	analysiscuts->SetRemoveDaughtersFromPrim(kFALSE);
ac8cab2d	245	// To be set only in case of strong pid
	246	//analysiscuts->SetMaxPtStrongPid(9999.);
	247	//analysiscuts->SetMaxPStrongPidK(1.5);
	248	//analysiscuts->SetMaxPStrongPidpi(1.5);
1f796bc2	249
	250	cout<<"This is the odject I'm going to save:"<<nptbins<<endl;
	251
	252	analysiscuts->PrintAll();
	253	TFile* fout=new TFile("DstoKKpiCuts.root","recreate");
	254	fout->cd();
	255	analysiscuts->Write();
	256	fout->Close();
	257
	258
	259	}
	260
	261
	262	void makeInputAliAnalysisTaskSESignificanceMaximization(){
	263
	264	AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
	265	esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
	266	//default
	267	esdTrackCuts->SetRequireTPCRefit(kTRUE);
	268	esdTrackCuts->SetRequireITSRefit(kTRUE);
	269	//esdTrackCuts->SetMinNClustersITS(4); // default is 5
	270	esdTrackCuts->SetMinNClustersTPC(70);
	271	esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
	272	AliESDtrackCuts::kAny);
	273	// default is kBoth, otherwise kAny
	274	esdTrackCuts->SetMinDCAToVertexXY(0.);
	275	esdTrackCuts->SetPtRange(0.3,1.e10);
	276
	277
	278
	279	AliRDHFCutsDstoKKpi* RDHFDstoKKpi=new AliRDHFCutsDstoKKpi();
	280	RDHFDstoKKpi->SetName("loosercuts");
	281	RDHFDstoKKpi->SetTitle("Cuts for significance maximization");
	282	RDHFDstoKKpi->AddTrackCuts(esdTrackCuts);
	283	RDHFDstoKKpi->SetUsePID(kTRUE);
	284	RDHFDstoKKpi->SetPidOption(1);
	285	RDHFDstoKKpi->SetOptPileup(kTRUE);
	286	RDHFDstoKKpi->SetRemoveDaughtersFromPrim(kTRUE);
	287
ac8cab2d	288	const Int_t nvars=20;
1f796bc2	289
	290	const Int_t nptbins=4;
	291	Float_t* ptbins;
	292	ptbins=new Float_t[nptbins+1];
	293	ptbins[0]=2.;
	294	ptbins[1]=4.;
	295	ptbins[2]=6.;
	296	ptbins[3]=8.;
	297	ptbins[4]=12.;
	298
	299	RDHFDstoKKpi->SetPtBins(nptbins+1,ptbins);
	300
	301	//setting my cut values
	302
	303	Float_t** prodcutsval;
	304	prodcutsval=new Float_t*[nvars];
	305	for(Int_t ic=0;ic<nvars;ic++){prodcutsval[ic]=new Float_t[nptbins];}
	306	for(Int_t ipt=0;ipt<nptbins;ipt++){
	307
	308	prodcutsval[0][ipt]=0.35;
	309	prodcutsval[1][ipt]=0.3;
	310	prodcutsval[2][ipt]=0.3;
	311	prodcutsval[3][ipt]=0.;
	312	prodcutsval[4][ipt]=0.;
	313	prodcutsval[5][ipt]=0.005;
	314	prodcutsval[6][ipt]=0.06;
	315	prodcutsval[7][ipt]=0.0;
	316	prodcutsval[8][ipt]=0.;
	317	prodcutsval[9][ipt]=0.9;
	318	prodcutsval[10][ipt]=0.;
	319	prodcutsval[11][ipt]=1000.0;
	320	prodcutsval[12][ipt]=0.015;
	321	prodcutsval[13][ipt]=0.1;
	322	prodcutsval[14][ipt]=0.;
	323	prodcutsval[15][ipt]=1.;
ac8cab2d	324	prodcutsval[16][ipt]=0.;
	325	prodcutsval[17][ipt]=0.;
	326	prodcutsval[18][ipt]=0.;
	327	prodcutsval[19][ipt]=-1.;
	328
1f796bc2	329	}
	330
	331	RDHFDstoKKpi->SetCuts(nvars,nptbins,prodcutsval);
	332
	333	Int_t nvarsforopt=RDHFDstoKKpi->GetNVarsForOpt();
	334	//Int_t nvarsforopt=2;
	335	const Int_t dim=4; //set this!!
	336	Bool_t *boolforopt;
	337	boolforopt=new Bool_t[nvars];
	338	if(dim>nvarsforopt){
	339	cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
	340	return;
	341	} else {
	342	if(dim==nvarsforopt){
	343	boolforopt=RDHFDstoKKpi->GetVarsForOpt();
	344	}else{
	345	TString *names;
	346	names=new TString[nvars];
	347	TString answer="";
	348	Int_t checktrue=0;
	349	names=RDHFDstoKKpi->GetVarNames();
	350	for(Int_t i=0;i<nvars;i++){
	351	cout<<names[i]<<" for opt? (y/n)"<<endl;
	352	cin>>answer;
	353	if(answer=="y") {
	354	boolforopt[i]=kTRUE;
	355	checktrue++;
	356	}
	357	else boolforopt[i]=kFALSE;
	358	}
	359	if (checktrue!=dim) {
	360	cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
	361	return;
	362	}
	363	RDHFDstoKKpi->SetVarsForOpt(dim,boolforopt);
	364	}
	365	}
	366
	367
	368	Float_t tighterval[dim][nptbins];
	369
	370	//number of steps for each variable is 4 now
	371	//set this!!
	372
	373	/*
	374	Cut list rejection condition
	375
	376	0 "inv. mass [GeV]", invmassDS-massDspdg>fCutsRD
	377	1 "pTK [GeV/c]", pTK<fCutsRd
	378	2 "pTPi [GeV/c]", pTPi<fCutsRd
	379	3 "d0K [cm]", d0K<fCutsRd
	380	4 "d0Pi [cm]", d0Pi<fCutsRd
	381	5 "dist12 [cm]", dist12<fCutsRd
	382	6 "sigmavert [cm]", sigmavert>fCutsRd
	383	7 "decLen [cm]", decLen<fCutsRD
	384	8 "ptMax [GeV/c]", ptMax<fCutsRD
	385	9 "cosThetaPoint", CosThetaPoint<fCutsRD
	386	10 "Sum d0^2 (cm^2)", sumd0<fCutsRD
	387	11 "dca [cm]", dca(i)>fCutsRD
	388	12 "inv. mass (Mphi-MKK) [GeV]", invmass-pdg>fCutsRD
	389	13 "inv. mass (MKo*-MKpi) [GeV]"}; invmass-pdg>fCutsRD
	390	14 "Abs(CosineKpiPhiRFrame)^3",
	391	15 "CosPiDsLabFrame"}
	392	*/
393
394	//sigmavert [cm]
395
396
397	tighterval[0][0]=0.0;
398	tighterval[0][1]=0.0;
399	tighterval[0][2]=0.0;
400	tighterval[0][3]=0.0;
401
402	//decay length
403
404	tighterval[1][0]=0.05;
405	tighterval[1][1]=0.05;
406	tighterval[1][2]=0.05;
407	tighterval[1][3]=0.05;
408
409	//costhetapoint
410
411	tighterval[2][0]=1.;
412	tighterval[2][1]=1.;
413	tighterval[2][2]=1.;
414	tighterval[2][3]=1.;
415
416
417	//inv mass phi meson
418
419	tighterval[3][0]=0.;
420	tighterval[3][1]=0.;
421	tighterval[3][2]=0.;
422	tighterval[3][3]=0.;
423
424	TString name="";
425	Int_t arrdim=dim*nptbins;
426	TClonesArray max("TParameter<float>",arrdim);
427	for (Int_t ipt=0;ipt<nptbins;ipt++){
428	for(Int_t ival=0;ival<dim;ival++){
429	name=Form("par%dptbin%d",ival,ipt);
430	new(max[ipt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][ipt]);
431	}
432	}
433
434	TFile* fout=new TFile("cuts4SignifMaximDs.root","recreate"); //set this!!
435	fout->cd();
436	RDHFDstoKKpi->Write();
437	max.Write();
438	fout->Close();
439
440	}