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

#include <Riostream.h>
#include <TFile.h>
#include <AliRDHFCutsD0toKpi.h>
#include <TClonesArray.h>
#include <TParameter.h>


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

//Author: Chiara Bianchin, cbianchi@pd.infn.it


//macro to make a .root file which contains an AliRDHFCutsD0toKpi for AliAnalysisTaskSED0Mass task

void makeInputAliAnalysisTaskSED0Mass(){

  AliRDHFCutsD0toKpi* RDHFD0toKpi=new AliRDHFCutsD0toKpi();
  RDHFD0toKpi->SetName("D0toKpiCuts");
  RDHFD0toKpi->SetTitle("Cuts for D0 analysis");

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


  RDHFD0toKpi->AddTrackCuts(esdTrackCuts);

  const Int_t nvars=11;

  const Int_t nptbins=13;
  Float_t* ptbins;
  ptbins=new Float_t[nptbins+1];
  ptbins[0]=0.;
  ptbins[1]=0.5;
  ptbins[2]=1.;
  ptbins[3]=2.;
  ptbins[4]=3.;
  ptbins[5]=4.;
  ptbins[6]=5.;
  ptbins[7]=6.;
  ptbins[8]=8.;
  ptbins[9]=12.;
  ptbins[10]=16.;
  ptbins[11]=20.;
  ptbins[12]=24.;
  ptbins[13]=9999.;

  RDHFD0toKpi->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
    //cuts order
    //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
    //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
    //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
    //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
    //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
    //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
    //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
    //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
    //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
    //     printf("    |cosThetaPointXY| < %f\n",fD0toKpiCuts[9]);
    //     printf("    NormDecayLenghtXY    > %f\n",fD0toKpiCuts[10]);


  Double_t arrcuts[11]={0.3,0.03,0.8,0.8,0.8,0.1,0.1,-0.0004,0.9,0.998,5.}; //put the last 2 values at 0. for pp

  //setting my cut values
  //0-0.5 GeV
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][0]=arrcuts[ic];

  //0.5-1 GeV/c
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][1]=arrcuts[ic];

  //1-2 GeV 
  arrcuts[1]=0.025;  arrcuts[7]=-0.0003;
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][2]=arrcuts[ic];

  //2-3 GeV
  arrcuts[7]=-0.00026;
  arrcuts[9]=0.998;
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][3]=arrcuts[ic];

  //3-4 GeV
  arrcuts[7]=-0.00015; arrcuts[8]=0.85;
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][4]=arrcuts[ic];

  //4-5 GeV 
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][5]=arrcuts[ic];

  //5-6 GeV
  arrcuts[7]=-0.0001; arrcuts[8]=0.85;
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][6]=arrcuts[ic];

  //6-8 GeV
  arrcuts[2]=1.;
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][7]=arrcuts[ic];

  //8-12 GeV
  arrcuts[8]=0.8;
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][8]=arrcuts[ic];

  //12-16 GeV
  arrcuts[1]=0.03;
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][9]=arrcuts[ic];

  //16-20 GeV
  arrcuts[1]=0.035;
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][10]=arrcuts[ic];

  //20-24 GeV
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][11]=arrcuts[ic];

  //24-9999 GeV
  for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][12]=arrcuts[ic];

  RDHFD0toKpi->SetCuts(nvars,nptbins,rdcutsvalmine);

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

    //pid settings
  AliAODPidHF* pidObj=new AliAODPidHF();
  //pidObj->SetName("pid4D0");
  Int_t mode=1;
  const Int_t nlims=2;
  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
  Bool_t compat=kTRUE; //effective only for this mode
  Bool_t asym=kTRUE;
  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
  pidObj->SetMatch(mode);
  pidObj->SetPLimit(plims,nlims);
  pidObj->SetSigma(sigmas);
  pidObj->SetCompat(compat);
  pidObj->SetTPC(kTRUE);
  pidObj->SetTOF(kTRUE);
  RDHFD0toKpi->SetPidHF(pidObj);

  RDHFD0toKpi->SetUseDefaultPID(kFALSE); //to use the AliAODPidHF

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

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

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

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

  cout<<"This is the odject I'm going to save:"<<endl;
  RDHFD0toKpi->PrintAll();
  TFile* fout=new TFile(Form("D0toKpiCuts%s%s%sRecVtx%sPileupRej.root", RDHFD0toKpi->GetUseCentrality()==0 ? "pp" : "PbPb",cent.Data(),RDHFD0toKpi->GetIsPrimaryWithoutDaughters() ? "" : "No",RDHFD0toKpi->GetOptPileUp() ? "" : "No"),"recreate");   //set this!! 

  fout->cd();
  RDHFD0toKpi->Write();
  fout->Close();

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

void makeInputAliAnalysisTaskSESignificanceMaximization(){

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

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

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

  const Int_t nvars=11;

  const Int_t nptbins=14; //change this when adding pt bins!
  Float_t ptbins[nptbins+1];
  ptbins[0]=0.;
  ptbins[1]=0.5;
  ptbins[2]=1.;
  ptbins[3]=2.;
  ptbins[4]=3.;
  ptbins[5]=4.;
  ptbins[6]=5.;
  ptbins[7]=6.;
  ptbins[8]=8.;
  ptbins[9]=10.;
  ptbins[10]=12.;
  ptbins[11]=16.;
  ptbins[12]=20.;
  ptbins[13]=24.;
  ptbins[14]=9999.;

  RDHFD0toKpi->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
    //cuts order
    //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
    //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
    //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
    //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
    //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
    //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
    //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
    //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
    //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
    //     printf("    |cosThetaPointXY| < %f\n",fD0toKpiCuts[9]);
    //     printf("    NormDecayLenghtXY    > %f\n",fD0toKpiCuts[10]);


    Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.75,0.,2.},/* pt<0.5*/
						  {0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.75,0.,2.},/* 0.5<pt<1*/
						  {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-33000.*1E-8,0.75,0.,2.},/* 1<pt<2 */
						  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-30000.*1E-8,0.85,0.994,2.},/* 2<pt<3 */
						  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-26000.*1E-8,0.85,0.994,2.},/* 3<pt<4 */
						  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.85,0.994,2.},/* 4<pt<5 */
						  {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-11000.*1E-8,0.82,0.994,2.},/* 5<pt<6 */
						  {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.78,0.994,2.},/* 6<pt<8 */
						  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.7,0.994,2.},/* 8<pt<10 */
						  {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.7,0.994,2.},/* 10<pt<12 */
						  {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-500.*1E-8,0.7,0.994,2.},/* 12<pt<16 */
						  {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-500.*1E-8,0.7,0.994,2.},/* 16<pt<20 */
						  {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-500.*1E-8,0.7,0.994,2.},/* 20<pt<24 */
						  {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-500.*1E-8,0.7,0.994,2.}};/* pt>24 */

  //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]=cutsMatrixD0toKpiStand[ibin][ivar];
    }
  }
  RDHFD0toKpi->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);


  Bool_t stdvaropt=kFALSE;
  Int_t dim=4; //set this!!
  Bool_t *boolforopt;
  boolforopt=new Bool_t[nvars];
  if(stdvaropt){
    boolforopt=RDHFD0toKpi->GetVarsForOpt();
  }else{
    TString *names;
    names=new TString[nvars];
    TString answer="";
    Int_t checktrue=0;
    names=RDHFD0toKpi->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;
    }
    RDHFD0toKpi->SetVarsForOpt(dim,boolforopt);
  }


  Float_t tighterval[dim][nptbins];
  //dca  
  //costhetastar
  //d0d0 <-this 
  //costhetapoint <-this 
  //cosThetaPointXY <-this 
  //NormDecayLength <-this 
  
  //number of steps for each variable is set in the AddTask arguments (default=8)
  //set this!!
  //0-0.5
  tighterval[0][0]=-0.00065;
  tighterval[1][0]=1.;
  tighterval[2][0]=0.4;
  tighterval[3][0]=6.;

  //0.5-1.
  tighterval[0][1]=-0.00065;
  tighterval[1][1]=1.;
  tighterval[2][1]=0.4;
  tighterval[3][1]=6.;

  //1-2
  tighterval[0][2]=-0.00065;
  tighterval[1][2]=1.;
  tighterval[2][2]=0.4;
  tighterval[3][2]=6.;
 
  //2-3
  tighterval[0][3]=-0.0006;
  tighterval[1][3]=1.;
  tighterval[2][3]=1.;
  tighterval[3][3]=6.;

  //3-4
  tighterval[0][4]=-0.00046;
  tighterval[1][4]=1.;
  tighterval[2][4]=1.;
  tighterval[3][4]=6.;
 
  //4-5
  tighterval[0][5]=-0.00045;
  tighterval[1][5]=1.;
  tighterval[2][5]=1.;
  tighterval[3][5]=6.;
 
  //5-6
  tighterval[0][6]=-0.00031;
  tighterval[1][6]=1.;
  tighterval[2][6]=1.;
  tighterval[3][6]=6.;

  //6-8
  tighterval[0][7]=-0.00021;
  tighterval[1][7]=0.98;
  tighterval[2][7]=1.;
  tighterval[3][7]=6.;

  //8-10
  tighterval[0][8]=-0.0001;
  tighterval[1][8]=0.98;
  tighterval[2][8]=1.;
  tighterval[3][8]=6.;

  //10-12
  tighterval[0][9]=-0.0001;
  tighterval[1][9]=0.9;
  tighterval[2][9]=1.;
  tighterval[3][9]=6.;
 
  //12-16
  tighterval[0][10]=-0.00005;
  tighterval[1][10]=0.9;
  tighterval[2][10]=1.;
  tighterval[3][10]=6.;

  //16-20
  tighterval[0][11]=-0.00005;
  tighterval[1][11]=0.9;
  tighterval[2][11]=1.;
  tighterval[3][11]=6.;

  //20-24
  tighterval[0][12]=-0.00005;
  tighterval[1][12]=0.9;
  tighterval[2][12]=1.;
  tighterval[3][12]=6.;

  //>24
  tighterval[0][13]=-0.00005;
  tighterval[1][13]=0.9;
  tighterval[2][13]=1.;
  tighterval[3][13]=6.;


  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=kTRUE;
  RDHFD0toKpi->SetUsePID(flagPID);

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

  //pid settings
  AliAODPidHF* pidObj=new AliAODPidHF();
  //pidObj->SetName("pid4D0");
  Int_t mode=1;
  const Int_t nlims=2;
  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
  Bool_t compat=kTRUE; //effective only for this mode
  Bool_t asym=kTRUE;
  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
  pidObj->SetMatch(mode);
  pidObj->SetPLimit(plims,nlims);
  pidObj->SetSigma(sigmas);
  pidObj->SetCompat(compat);
  pidObj->SetTPC(kTRUE);
  pidObj->SetTOF(kTRUE);
  RDHFD0toKpi->SetPidHF(pidObj);

  RDHFD0toKpi->SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
  
  //activate pileup rejection (for pp)
  //RDHFD0toKpi->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);

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

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

  //temporary
  RDHFD0toKpi->SetFixRefs();

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

  fout->cd();
  RDHFD0toKpi->Write();
  max.Write();
  fout->Close();
 
}
Commit	Line	Data
d52f7b50	1	#include <Riostream.h>
	2	#include <TFile.h>
	3	#include <AliRDHFCutsD0toKpi.h>
	4	#include <TClonesArray.h>
	5	#include <TParameter.h>
cbedddce	6
cbedddce	7
	8	//Use:
	9	//Set hard coded commentet with //set this!!
3b79985d	10	// root[] .L makeInputD0tasks.C++
	11	// root[] makeInputAliAnalysisTaskSED0Mass()
	12	// root[] makeInputAliAnalysisTaskSESignificanceMaximization()
cbedddce	13	//similar macros for the other D mesons
	14
	15	//Author: Chiara Bianchin, cbianchi@pd.infn.it
	16
3b79985d	17
	18	//macro to make a .root file which contains an AliRDHFCutsD0toKpi for AliAnalysisTaskSED0Mass task
	19
	20	void makeInputAliAnalysisTaskSED0Mass(){
	21
	22	AliRDHFCutsD0toKpi* RDHFD0toKpi=new AliRDHFCutsD0toKpi();
	23	RDHFD0toKpi->SetName("D0toKpiCuts");
	24	RDHFD0toKpi->SetTitle("Cuts for D0 analysis");
	25
	26	AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
	27	esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
	28	//default
	29	esdTrackCuts->SetRequireTPCRefit(kTRUE);
	30	esdTrackCuts->SetRequireITSRefit(kTRUE);
44fb0179	31	esdTrackCuts->SetEtaRange(-0.8,0.8);
0d4398a2	32	//esdTrackCuts->SetMinNClustersITS(4); // default is 5
44fb0179	33	//esdTrackCuts->SetMinNClustersTPC(120);
3b79985d	34	esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
	35	AliESDtrackCuts::kAny);
	36	// default is kBoth, otherwise kAny
	37	esdTrackCuts->SetMinDCAToVertexXY(0.);
0beb4580	38	esdTrackCuts->SetPtRange(0.8,1.e10);
3b79985d	39
	40
	41	RDHFD0toKpi->AddTrackCuts(esdTrackCuts);
	42
29068278	43	const Int_t nvars=11;
3b79985d	44
0beb4580	45	const Int_t nptbins=13;
3b79985d	46	Float_t* ptbins;
	47	ptbins=new Float_t[nptbins+1];
	48	ptbins[0]=0.;
0beb4580	49	ptbins[1]=0.5;
	50	ptbins[2]=1.;
	51	ptbins[3]=2.;
	52	ptbins[4]=3.;
	53	ptbins[5]=4.;
	54	ptbins[6]=5.;
	55	ptbins[7]=6.;
	56	ptbins[8]=8.;
	57	ptbins[9]=12.;
	58	ptbins[10]=16.;
	59	ptbins[11]=20.;
	60	ptbins[12]=24.;
	61	ptbins[13]=9999.;
	62
3b79985d	63	RDHFD0toKpi->SetPtBins(nptbins+1,ptbins);
	64
	65
	66	Float_t** rdcutsvalmine;
	67	rdcutsvalmine=new Float_t*[nvars];
	68	for(Int_t iv=0;iv<nvars;iv++){
	69	rdcutsvalmine[iv]=new Float_t[nptbins];
	70	}
	71
	72	//setting my cut values
	73	//cuts order
	74	// printf(" \|M-MD0\| [GeV] < %f\n",fD0toKpiCuts[0]);
	75	// printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
	76	// printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
	77	// printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
	78	// printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
	79	// printf(" \|d0K\| [cm] < %f\n",fD0toKpiCuts[5]);
	80	// printf(" \|d0pi\| [cm] < %f\n",fD0toKpiCuts[6]);
	81	// printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
	82	// printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
29068278	83	// printf(" \|cosThetaPointXY\| < %f\n",fD0toKpiCuts[9]);
29068278	84	// printf(" NormDecayLenghtXY > %f\n",fD0toKpiCuts[10]);
3b79985d	85
0beb4580	86
29068278	87	Double_t arrcuts[11]={0.3,0.03,0.8,0.8,0.8,0.1,0.1,-0.0004,0.9,0.998,5.}; //put the last 2 values at 0. for pp
0beb4580	88
3b79985d	89	//setting my cut values
0beb4580	90	//0-0.5 GeV
	91	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][0]=arrcuts[ic];
	92
	93	//0.5-1 GeV/c
	94	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][1]=arrcuts[ic];
6b3e3c78	95
0beb4580	96	//1-2 GeV
	97	arrcuts[1]=0.025; arrcuts[7]=-0.0003;
	98	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][2]=arrcuts[ic];
6b3e3c78	99
6b3e3c78	100	//2-3 GeV
0beb4580	101	arrcuts[7]=-0.00026;
29068278	102	arrcuts[9]=0.998;
0beb4580	103	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][3]=arrcuts[ic];
	104
	105	//3-4 GeV
	106	arrcuts[7]=-0.00015; arrcuts[8]=0.85;
	107	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][4]=arrcuts[ic];
	108
	109	//4-5 GeV
	110	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][5]=arrcuts[ic];
	111
	112	//5-6 GeV
	113	arrcuts[7]=-0.0001; arrcuts[8]=0.85;
	114	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][6]=arrcuts[ic];
	115
	116	//6-8 GeV
	117	arrcuts[2]=1.;
	118	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][7]=arrcuts[ic];
	119
	120	//8-12 GeV
	121	arrcuts[8]=0.8;
	122	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][8]=arrcuts[ic];
	123
	124	//12-16 GeV
	125	arrcuts[1]=0.03;
	126	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][9]=arrcuts[ic];
	127
	128	//16-20 GeV
	129	arrcuts[1]=0.035;
	130	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][10]=arrcuts[ic];
	131
	132	//20-24 GeV
	133	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][11]=arrcuts[ic];
	134
	135	//24-9999 GeV
	136	for(Int_t ic=0;ic<nvars;ic++) rdcutsvalmine[ic][12]=arrcuts[ic];
3b79985d	137
3b79985d	138	RDHFD0toKpi->SetCuts(nvars,nptbins,rdcutsvalmine);
6b3e3c78	139
	140	Bool_t pidflag=kTRUE;
	141	RDHFD0toKpi->SetUsePID(pidflag);
	142	if(pidflag) cout<<"PID is used"<<endl;
	143	else cout<<"PID is not used"<<endl;
	144
d52f7b50	145	//pid settings
6b3e3c78	146	AliAODPidHF* pidObj=new AliAODPidHF();
	147	//pidObj->SetName("pid4D0");
	148	Int_t mode=1;
	149	const Int_t nlims=2;
	150	Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
	151	Bool_t compat=kTRUE; //effective only for this mode
	152	Bool_t asym=kTRUE;
	153	Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
	154	pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
	155	pidObj->SetMatch(mode);
	156	pidObj->SetPLimit(plims,nlims);
	157	pidObj->SetSigma(sigmas);
	158	pidObj->SetCompat(compat);
	159	pidObj->SetTPC(kTRUE);
	160	pidObj->SetTOF(kTRUE);
	161	RDHFD0toKpi->SetPidHF(pidObj);
	162
	163	RDHFD0toKpi->SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
	164
0beb4580	165	//activate pileup rejection (for pp)
	166	//RDHFD0toKpi->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
	167
	168	//Do not recalculate the vertex
	169	RDHFD0toKpi->SetRemoveDaughtersFromPrim(kFALSE); //activate for pp
	170
	171	TString cent="";
	172	//centrality selection (Pb-Pb)
	173	Float_t minc=20,maxc=80;
	174	RDHFD0toKpi->SetMinCentrality(minc);
	175	RDHFD0toKpi->SetMaxCentrality(maxc);
	176	cent=Form("%.0f%.0f",minc,maxc);
	177	RDHFD0toKpi->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
	178
	179	//temporary
29068278	180	//RDHFD0toKpi->SetFixRefs();
d52f7b50	181
3b79985d	182	cout<<"This is the odject I'm going to save:"<<endl;
3b79985d	183	RDHFD0toKpi->PrintAll();
0beb4580	184	TFile* fout=new TFile(Form("D0toKpiCuts%s%s%sRecVtx%sPileupRej.root", RDHFD0toKpi->GetUseCentrality()==0 ? "pp" : "PbPb",cent.Data(),RDHFD0toKpi->GetIsPrimaryWithoutDaughters() ? "" : "No",RDHFD0toKpi->GetOptPileUp() ? "" : "No"),"recreate"); //set this!!
0beb4580	185
3b79985d	186	fout->cd();
	187	RDHFD0toKpi->Write();
	188	fout->Close();
	189
	190	}
29068278	191
3b79985d	192	//macro to make a .root file (for significance maximization) which contains an AliRDHFCutsD0toKpi with loose set of cuts and TParameter with the tighest value of these cuts
	193
	194	void makeInputAliAnalysisTaskSESignificanceMaximization(){
cbedddce	195
	196	AliRDHFCutsD0toKpi* RDHFD0toKpi=new AliRDHFCutsD0toKpi();
	197	RDHFD0toKpi->SetName("loosercuts");
	198	RDHFD0toKpi->SetTitle("Cuts for significance maximization");
	199
	200	AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
	201	esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
	202	//default
	203	esdTrackCuts->SetRequireTPCRefit(kTRUE);
	204	esdTrackCuts->SetRequireITSRefit(kTRUE);
0d4398a2	205	//esdTrackCuts->SetMinNClustersITS(4);
44fb0179	206	//esdTrackCuts->SetMinNClustersTPC(120);
44fb0179	207
cbedddce	208	esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
cbedddce	209	esdTrackCuts->SetMinDCAToVertexXY(0.);
44fb0179	210	esdTrackCuts->SetEtaRange(-0.8,0.8);
0beb4580	211	esdTrackCuts->SetPtRange(0.8,1.e10);
cbedddce	212
	213	RDHFD0toKpi->AddTrackCuts(esdTrackCuts);
	214
29068278	215	const Int_t nvars=11;
cbedddce	216
29068278	217	const Int_t nptbins=14; //change this when adding pt bins!
cbedddce	218	Float_t ptbins[nptbins+1];
cbedddce	219	ptbins[0]=0.;
6b3e3c78	220	ptbins[1]=0.5;
	221	ptbins[2]=1.;
	222	ptbins[3]=2.;
	223	ptbins[4]=3.;
	224	ptbins[5]=4.;
	225	ptbins[6]=5.;
	226	ptbins[7]=6.;
	227	ptbins[8]=8.;
29068278	228	ptbins[9]=10.;
	229	ptbins[10]=12.;
	230	ptbins[11]=16.;
	231	ptbins[12]=20.;
	232	ptbins[13]=24.;
	233	ptbins[14]=9999.;
0beb4580	234
cbedddce	235	RDHFD0toKpi->SetPtBins(nptbins+1,ptbins);
cbedddce	236
	237	Float_t** rdcutsvalmine;
	238	rdcutsvalmine=new Float_t*[nvars];
	239	for(Int_t iv=0;iv<nvars;iv++){
	240	rdcutsvalmine[iv]=new Float_t[nptbins];
	241	}
	242
	243	//setting my cut values
	244	//cuts order
	245	// printf(" \|M-MD0\| [GeV] < %f\n",fD0toKpiCuts[0]);
	246	// printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
	247	// printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
	248	// printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
	249	// printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
	250	// printf(" \|d0K\| [cm] < %f\n",fD0toKpiCuts[5]);
	251	// printf(" \|d0pi\| [cm] < %f\n",fD0toKpiCuts[6]);
	252	// printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
	253	// printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
29068278	254	// printf(" \|cosThetaPointXY\| < %f\n",fD0toKpiCuts[9]);
	255	// printf(" NormDecayLenghtXY > %f\n",fD0toKpiCuts[10]);
	256
	257
	258	Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.1E-4,0.8,0.3,0.3,1000.1E-4,1000.1E-4,-40000.1E-8,0.75,0.,2.},/* pt<0.5*/
	259	{0.400,400.1E-4,0.8,0.3,0.3,1000.1E-4,1000.1E-4,-40000.1E-8,0.75,0.,2.},/* 0.5<pt<1*/
	260	{0.400,400.1E-4,0.8,0.4,0.4,1000.1E-4,1000.1E-4,-33000.1E-8,0.75,0.,2.},/* 1<pt<2 */
	261	{0.400,250.1E-4,0.8,0.7,0.7,1000.1E-4,1000.1E-4,-30000.1E-8,0.85,0.994,2.},/* 2<pt<3 */
	262	{0.400,250.1E-4,0.8,0.7,0.7,1000.1E-4,1000.1E-4,-26000.1E-8,0.85,0.994,2.},/* 3<pt<4 */
	263	{0.400,250.1E-4,0.8,0.7,0.7,1000.1E-4,1000.1E-4,-15000.1E-8,0.85,0.994,2.},/* 4<pt<5 */
	264	{0.400,250.1E-4,0.8,0.7,0.7,1000.1E-4,1000.1E-4,-11000.1E-8,0.82,0.994,2.},/* 5<pt<6 */
	265	{0.400,270.1E-4,0.8,0.7,0.7,1000.1E-4,1000.1E-4,-10000.1E-8,0.78,0.994,2.},/* 6<pt<8 */
	266	{0.400,300.1E-4,0.8,0.7,0.7,1000.1E-4,1000.1E-4,-1000.1E-8,0.7,0.994,2.},/* 8<pt<10 */
	267	{0.400,300.1E-4,0.8,0.7,0.7,1000.1E-4,1000.1E-4,-1000.1E-8,0.7,0.994,2.},/* 10<pt<12 */
	268	{0.400,350.1E-4,1.0,0.7,0.7,1000.1E-4,1000.1E-4,-500.1E-8,0.7,0.994,2.},/* 12<pt<16 */
	269	{0.400,400.1E-4,1.0,0.7,0.7,1000.1E-4,1000.1E-4,-500.1E-8,0.7,0.994,2.},/* 16<pt<20 */
	270	{0.400,400.1E-4,1.0,0.7,0.7,1000.1E-4,1000.1E-4,-500.1E-8,0.7,0.994,2.},/* 20<pt<24 */
	271	{0.400,400.1E-4,1.0,0.7,0.7,1000.1E-4,1000.1E-4,-500.1E-8,0.7,0.994,2.}};/* pt>24 */
d52f7b50	272
6b3e3c78	273	//CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
d52f7b50	274	Float_t *cutsMatrixTransposeStand=new Float_t[nvars];
d52f7b50	275	for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
6b3e3c78	276	for (Int_t ibin=0;ibin<nptbins;ibin++){
6b3e3c78	277	for (Int_t ivar = 0; ivar<nvars; ivar++){
d52f7b50	278	cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
6b3e3c78	279	}
6b3e3c78	280	}
d52f7b50	281	RDHFD0toKpi->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
cbedddce	282
cbedddce	283
29068278	284	Bool_t stdvaropt=kFALSE;
29068278	285	Int_t dim=4; //set this!!
cbedddce	286	Bool_t *boolforopt;
cbedddce	287	boolforopt=new Bool_t[nvars];
29068278	288	if(stdvaropt){
	289	boolforopt=RDHFD0toKpi->GetVarsForOpt();
	290	}else{
	291	TString *names;
	292	names=new TString[nvars];
	293	TString answer="";
	294	Int_t checktrue=0;
	295	names=RDHFD0toKpi->GetVarNames();
	296	for(Int_t i=0;i<nvars;i++){
	297	cout<<names[i]<<" for opt? (y/n)"<<endl;
	298	cin>>answer;
	299	if(answer=="y") {
	300	boolforopt[i]=kTRUE;
	301	checktrue++;
cbedddce	302	}
29068278	303	else boolforopt[i]=kFALSE;
cbedddce	304	}
29068278	305	if (checktrue!=dim) {
	306	cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
	307	return;
	308	}
	309	RDHFD0toKpi->SetVarsForOpt(dim,boolforopt);
cbedddce	310	}
	311
	312
	313	Float_t tighterval[dim][nptbins];
44fb0179	314	//dca
cbedddce	315	//costhetastar
	316	//d0d0 <-this
	317	//costhetapoint <-this
29068278	318	//cosThetaPointXY <-this
29068278	319	//NormDecayLength <-this
cbedddce	320
44fb0179	321	//number of steps for each variable is set in the AddTask arguments (default=8)
cbedddce	322	//set this!!
29068278	323	//0-0.5
44fb0179	324	tighterval[0][0]=-0.00065;
44fb0179	325	tighterval[1][0]=1.;
29068278	326	tighterval[2][0]=0.4;
29068278	327	tighterval[3][0]=6.;
cbedddce	328
29068278	329	//0.5-1.
44fb0179	330	tighterval[0][1]=-0.00065;
44fb0179	331	tighterval[1][1]=1.;
29068278	332	tighterval[2][1]=0.4;
29068278	333	tighterval[3][1]=6.;
44fb0179	334
29068278	335	//1-2
44fb0179	336	tighterval[0][2]=-0.00065;
44fb0179	337	tighterval[1][2]=1.;
29068278	338	tighterval[2][2]=0.4;
29068278	339	tighterval[3][2]=6.;
cbedddce	340
29068278	341	//2-3
29068278	342	tighterval[0][3]=-0.0006;
44fb0179	343	tighterval[1][3]=1.;
29068278	344	tighterval[2][3]=1.;
29068278	345	tighterval[3][3]=6.;
cbedddce	346
29068278	347	//3-4
29068278	348	tighterval[0][4]=-0.00046;
44fb0179	349	tighterval[1][4]=1.;
29068278	350	tighterval[2][4]=1.;
29068278	351	tighterval[3][4]=6.;
44fb0179	352
29068278	353	//4-5
29068278	354	tighterval[0][5]=-0.00045;
44fb0179	355	tighterval[1][5]=1.;
29068278	356	tighterval[2][5]=1.;
29068278	357	tighterval[3][5]=6.;
44fb0179	358
29068278	359	//5-6
29068278	360	tighterval[0][6]=-0.00031;
44fb0179	361	tighterval[1][6]=1.;
29068278	362	tighterval[2][6]=1.;
	363	tighterval[3][6]=6.;
	364
	365	//6-8
	366	tighterval[0][7]=-0.00021;
	367	tighterval[1][7]=0.98;
	368	tighterval[2][7]=1.;
	369	tighterval[3][7]=6.;
	370
	371	//8-10
	372	tighterval[0][8]=-0.0001;
	373	tighterval[1][8]=0.98;
	374	tighterval[2][8]=1.;
	375	tighterval[3][8]=6.;
	376
	377	//10-12
	378	tighterval[0][9]=-0.0001;
	379	tighterval[1][9]=0.9;
	380	tighterval[2][9]=1.;
	381	tighterval[3][9]=6.;
44fb0179	382
29068278	383	//12-16
	384	tighterval[0][10]=-0.00005;
	385	tighterval[1][10]=0.9;
	386	tighterval[2][10]=1.;
	387	tighterval[3][10]=6.;
	388
	389	//16-20
	390	tighterval[0][11]=-0.00005;
	391	tighterval[1][11]=0.9;
	392	tighterval[2][11]=1.;
	393	tighterval[3][11]=6.;
	394
	395	//20-24
	396	tighterval[0][12]=-0.00005;
	397	tighterval[1][12]=0.9;
	398	tighterval[2][12]=1.;
	399	tighterval[3][12]=6.;
	400
	401	//>24
	402	tighterval[0][13]=-0.00005;
	403	tighterval[1][13]=0.9;
	404	tighterval[2][13]=1.;
	405	tighterval[3][13]=6.;
44fb0179	406
44fb0179	407
cbedddce	408	TString name="";
	409	Int_t arrdim=dim*nptbins;
	410	cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
	411	TClonesArray max("TParameter<float>",arrdim);
3b79985d	412	for(Int_t ival=0;ival<dim;ival++){
	413	for(Int_t jpt=0;jpt<nptbins;jpt++){
	414	name=Form("par%dptbin%d",ival,jpt);
	415	cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
	416	new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
cbedddce	417	}
cbedddce	418	}
0beb4580	419
6b3e3c78	420	Bool_t flagPID=kTRUE;
6b3e3c78	421	RDHFD0toKpi->SetUsePID(flagPID);
0beb4580	422
d52f7b50	423	RDHFD0toKpi->PrintAll();
6b3e3c78	424	printf("Use PID? %s\n",flagPID ? "yes" : "no");
	425
	426	//pid settings
	427	AliAODPidHF* pidObj=new AliAODPidHF();
	428	//pidObj->SetName("pid4D0");
	429	Int_t mode=1;
	430	const Int_t nlims=2;
	431	Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
	432	Bool_t compat=kTRUE; //effective only for this mode
	433	Bool_t asym=kTRUE;
	434	Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
	435	pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
	436	pidObj->SetMatch(mode);
	437	pidObj->SetPLimit(plims,nlims);
	438	pidObj->SetSigma(sigmas);
	439	pidObj->SetCompat(compat);
	440	pidObj->SetTPC(kTRUE);
	441	pidObj->SetTOF(kTRUE);
	442	RDHFD0toKpi->SetPidHF(pidObj);
	443
6b3e3c78	444	RDHFD0toKpi->SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
d52f7b50	445
0beb4580	446	//activate pileup rejection (for pp)
	447	//RDHFD0toKpi->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
	448
	449	//Do not recalculate the vertex
	450	RDHFD0toKpi->SetRemoveDaughtersFromPrim(kFALSE); //activate for pp
	451
	452	TString cent="";
	453	//centrality selection (Pb-Pb)
	454	Float_t minc=20,maxc=80;
	455	RDHFD0toKpi->SetMinCentrality(minc);
	456	RDHFD0toKpi->SetMaxCentrality(maxc);
	457	cent=Form("%.0f%.0f",minc,maxc);
	458	RDHFD0toKpi->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
d52f7b50	459
0beb4580	460	//temporary
0beb4580	461	RDHFD0toKpi->SetFixRefs();
6b3e3c78	462
0beb4580	463	TFile* fout=new TFile(Form("cuts4SignifMaxim%s%s%sRecVtx%sPileupRej.root", RDHFD0toKpi->GetUseCentrality()==0 ? "pp" : "PbPb",cent.Data(),RDHFD0toKpi->GetIsPrimaryWithoutDaughters() ? "" : "No",RDHFD0toKpi->GetOptPileUp() ? "" : "No"),"recreate"); //set this!!
44fb0179	464
cbedddce	465	fout->cd();
	466	RDHFD0toKpi->Write();
	467	max.Write();
	468	fout->Close();
	469
	470	}
3b79985d	471