PWGHF/vertexingHF/macros/makeTFile4CutsLctoV0bachelor.C

   1 #include <Riostream.h>
   2 #include <TFile.h>
   3 #include <AliRDHFCutsLctoV0.h>
   4 #include <AliAODPidHF.h>
   5 #include <TClonesArray.h>
   6 #include <TParameter.h>
   7 #include <TF1.h>
   8
   9
  10 //Use:
  11 //Set hard coded commentet with //set this!!
  12 // root[] .L makeInput...C++
  13 // root[] makeInputAliAnalysisTaskSE...()
  14 //similar macros for D mesons as well as for Lc->3prongs
  15
  16 //Author: Annalisa De Caro - decaro@sa.infn.it
  17
  18
  19 //macro to make a .root file which contains an AliRDHFCutsLctoV0 for AliAnalysisTaskSELc2pK0S task
  20
  21 void makeInputAliAnalysisTaskSELctoV0bachelor(){
  22
  23   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
  24   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
  25   //default
  26   esdTrackCuts->SetRequireTPCRefit(kTRUE);
  27   esdTrackCuts->SetRequireITSRefit(kTRUE);
  28   esdTrackCuts->SetMinNClustersITS(4); // default is 5
  29   esdTrackCuts->SetMinNClustersTPC(70);
  30   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
  31                                          AliESDtrackCuts::kAny);
  32   // default is kBoth, otherwise kAny
  33   esdTrackCuts->SetMinDCAToVertexXY(0.);
  34   esdTrackCuts->SetPtRange(0.3,1.e10);
  35
  36   AliRDHFCutsLctoV0* RDHFLctoV0An=new AliRDHFCutsLctoV0();
  37   RDHFLctoV0An->SetName("LctoV0AnalysisCuts");
  38   RDHFLctoV0An->SetTitle("Analysis cuts for Lc analysis");
  39
  40   AliRDHFCutsLctoV0* RDHFLctoV0Prod=new AliRDHFCutsLctoV0();
  41   RDHFLctoV0Prod->SetName("LctoV0ProductionCuts");
  42   RDHFLctoV0Prod->SetTitle("Production cuts for Lc analysis");
  43
  44   RDHFLctoV0Prod->AddTrackCuts(esdTrackCuts);
  45   RDHFLctoV0An->AddTrackCuts(esdTrackCuts);
  46
  47   const Int_t nptbins=1;
  48   Float_t* ptbins;
  49   ptbins=new Float_t[nptbins+1];
  50   ptbins[0]=0.;
  51   ptbins[1]=99999999.;
  52   RDHFLctoV0Prod->SetPtBins(nptbins+1,ptbins);
  53   RDHFLctoV0An->SetPtBins(nptbins+1,ptbins);
  54
  55   const Int_t nvars=9 ;
  56
  57   Float_t** prodcutsval;
  58   prodcutsval=new Float_t*[nvars];
  59   for(Int_t ic=0;ic<nvars;ic++){prodcutsval[ic]=new Float_t[nptbins];}
  60   for(Int_t ipt2=0;ipt2<nptbins;ipt2++){
  61    prodcutsval[0][ipt2]=1.; // inv. mass if K0S [GeV/c2]
  62    prodcutsval[1][ipt2]=1.; // inv. mass if Lambda [GeV/c2]
  63    prodcutsval[2][ipt2]=0.05; // inv. mass V0 if K0S [GeV/c2]
  64    prodcutsval[3][ipt2]=0.05; // inv. mass V0 if Lambda [GeV/c2]
  65    prodcutsval[4][ipt2]=0.3; // pT min bachelor track [GeV/c] // AOD by construction
  66    prodcutsval[5][ipt2]=0.;  // pT min V0-positive track [GeV/c]
  67    prodcutsval[6][ipt2]=0.;  // pT min V0-negative track [GeV/c]
  68    prodcutsval[7][ipt2]=1000.; // dca cascade cut [cm]
  69    prodcutsval[8][ipt2]=1000.; // dca V0 cut [nSigma] // it's 1.5 x offline V0s
  70   }
  71   RDHFLctoV0Prod->SetCuts(nvars,nptbins,prodcutsval);
  72
  73   Float_t** anacutsval;
  74   anacutsval=new Float_t*[nvars];
  75   for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
  76   for(Int_t ipt2=0;ipt2<nptbins;ipt2++){
  77    anacutsval[0][ipt2]=0.25; // inv. mass if K0S [GeV/c2]
  78    anacutsval[1][ipt2]=0.25; // inv. mass if Lambda [GeV/c2]
  79    anacutsval[2][ipt2]=0.0075; // inv. mass V0 if K0S [GeV/c2]
  80    anacutsval[3][ipt2]=0.0030; // inv. mass V0 if Lambda [GeV/c2]
  81    anacutsval[4][ipt2]=0.3; // pT min bachelor track [GeV/c] // AOD by construction
  82    anacutsval[5][ipt2]=0.;  // pT min V0-positive track [GeV/c]
  83    anacutsval[6][ipt2]=0.;  // pT min V0-negative track [GeV/c]
  84    anacutsval[7][ipt2]=1000.; // dca cascade cut [cm]
  85    anacutsval[8][ipt2]=1.5;   // dca V0 cut [nSigma] // it's 1.5 x offline V0s
  86   }
  87   RDHFLctoV0An->SetCuts(nvars,nptbins,anacutsval);
  88
  89
  90   //RDHFLc->SetRecoKF(); //set this if you want to recompute the secondary vertex with the KF package
  91
  92   //pid settings
  93   //1. bachelor: default one
  94   AliAODPidHF* pidObjBachelor = new AliAODPidHF();
  95   Double_t sigmasBac[5]={3.,1.,1.,3.,3.}; // 0, 1(A), 2(A) -> TPC; 3 -> TOF; 4 -> ITS
  96   pidObjBachelor->SetSigma(sigmasBac);
  97   pidObjBachelor->SetAsym(kFALSE);
  98   pidObjBachelor->SetMatch(1);
  99   pidObjBachelor->SetTPC(kTRUE);
 100   pidObjBachelor->SetTOF(kTRUE);
 101   pidObjBachelor->SetTOFdecide(kFALSE);
 102   RDHFLctoV0An->SetPidHF(pidObjBachelor);
 103
 104   //2. V0pos
 105   AliAODPidHF* pidObjV0pos = new AliAODPidHF();
 106   Double_t sigmasV0pos[5]={3.,1.,1.,3.,3.}; // 0, 1(A), 2(A) -> TPC; 3 -> TOF; 4 -> ITS
 107   pidObjV0pos->SetSigma(sigmasV0pos);
 108   pidObjV0pos->SetAsym(kFALSE);
 109   pidObjV0pos->SetMatch(1);
 110   pidObjV0pos->SetTPC(kTRUE);
 111   pidObjV0pos->SetTOF(kTRUE);
 112   pidObjV0pos->SetTOFdecide(kFALSE);
 113   RDHFLctoV0An->SetPidHF(pidObjV0pos);
 114
 115   //2. V0neg
 116   AliAODPidHF* pidObjV0neg = new AliAODPidHF();
 117   Double_t sigmasV0neg[5]={3.,1.,1.,3.,3.}; // 0, 1(A), 2(A) -> TPC; 3 -> TOF; 4 -> ITS
 118   pidObjV0neg->SetSigma(sigmasV0neg);
 119   pidObjV0neg->SetAsym(kFALSE);
 120   pidObjV0neg->SetMatch(1);
 121   pidObjV0neg->SetTPC(kTRUE);
 122   pidObjV0neg->SetTOF(kTRUE);
 123   pidObjV0neg->SetTOFdecide(kFALSE);
 124   RDHFLctoV0An->SetPidHF(pidObjV0neg);
 125
 126
 127   // uncomment these lines for Baysian PID:
 128   // Double_t threshold=0.3;
 129   // SetupCombinedPID(RDHFLctoV0An  ,threshold);
 130   // RDHFLctoV0An  ->SetPIDStrategy(AliRDHFCutsLctoV0::kCombined);
 131   //
 132
 133
 134   //uncomment these lines to apply cuts with the KF package
 135   //RDHFLctoV0An->SetCutsStrategy(AliRDHFCutsLctoV0::kKF);
 136   //for(Int_t ipt2=0;ipt2<nptbins;ipt2++){
 137   //   anacutsval[0][ipt2]=1.;  //if <0., no topological constraint
 138   //   anacutsval[1][ipt2]=2.;  //cut on the Chi2/Ndf
 139   // }
 140
 141   Bool_t pidflag=kTRUE;
 142   RDHFLctoV0An->SetUsePID(pidflag);
 143   if(pidflag) cout<<"PID is used"<<endl;
 144   else cout<<"PID is not used"<<endl;
 145
 146   cout<<"This is the object I'm going to save:"<<endl;
 147   RDHFLctoV0Prod->PrintAll();
 148   cout<<"This is the object I'm going to save:"<<endl;
 149   RDHFLctoV0An->PrintAll();
 150   TFile* fout=new TFile("Lc2pK0SCuts.root","RECREATE");
 151   fout->cd();
 152   RDHFLctoV0Prod->Write();
 153   RDHFLctoV0An->Write();
 154   fout->Close();
 155   delete fout;
 156
 157   delete pidObjBachelor;
 158   delete pidObjV0neg;
 159   delete pidObjV0pos;
 160   delete RDHFLctoV0Prod;
 161   delete RDHFLctoV0An;
 162
 163 }
 164
 165
 166 //macro to make a .root file (for significance maximization) which contains an AliRDHFCutsLctoV0 with loose set of cuts  and TParameter with the tighest value of these cuts
 167
 168 void makeInputAliAnalysisTaskSESignificanceMaximization(){
 169
 170   AliRDHFCutsLctoV0* RDHFLctoV0=new AliRDHFCutsLctoV0();
 171   RDHFLctoV0->SetName("loosercuts");
 172   RDHFLctoV0->SetTitle("Cuts for significance maximization");
 173
 174   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
 175   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 176   //default
 177   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 178   esdTrackCuts->SetMinNClustersTPC(70);
 179   esdTrackCuts->SetRequireITSRefit(kTRUE);
 180   esdTrackCuts->SetMinNClustersITS(4);
 181
 182   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 183   esdTrackCuts->SetMinDCAToVertexXY(0.);
 184   esdTrackCuts->SetEtaRange(-0.8,0.8);
 185   esdTrackCuts->SetPtRange(0.3,1.e10);
 186
 187   RDHFLctoV0->AddTrackCuts(esdTrackCuts);
 188
 189   const Int_t nvars=9;
 190
 191   const Int_t nptbins=1; //change this when adding pt bins!
 192   Float_t ptbins[nptbins+1];
 193   ptbins[0]=0.;
 194   ptbins[1]=999999999.;
 195
 196   RDHFLctoV0->SetPtBins(nptbins+1,ptbins);
 197
 198   Float_t** rdcutsvalmine;
 199   rdcutsvalmine=new Float_t*[nvars];
 200   for(Int_t iv=0;iv<nvars;iv++){
 201     rdcutsvalmine[iv]=new Float_t[nptbins];
 202   }
 203
 204   //setting my cut values
 205   // inv. mass if K0s [GeV]
 206   // inv. mass if Lambda [GeV]
 207   // inv. mass V0 if K0S [GeV]
 208   // inv. mass V0 if Lambda [GeV]
 209   // pT min bachelor track [GeV/c]
 210   // pT min V0-positive track [GeV/c]
 211   // pT min V0-negative track [GeV/c]
 212   // dca cascade cut (cm)
 213   // dca V0 cut (cm)
 214   Float_t cutsMatrixLctoV0Stand[nptbins][nvars]=
 215     { {0.25,0.25,0.0075,0.0075,0.,0.,0.,1.5,1.5} };
 216
 217   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
 218   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
 219   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
 220   for (Int_t ibin=0;ibin<nptbins;ibin++){
 221     for (Int_t ivar=0; ivar<nvars; ivar++){
 222       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixLctoV0Stand[ibin][ivar];
 223     }
 224   }
 225   RDHFLctoV0->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
 226
 227
 228   Int_t nvarsforopt=RDHFLctoV0->GetNVarsForOpt();
 229   Int_t dim=7; //set this!!
 230   Bool_t *boolforopt;
 231   boolforopt=new Bool_t[nvars];
 232   if(dim>nvarsforopt){
 233     cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
 234     return;
 235   } else {
 236     if(dim==nvarsforopt){
 237       boolforopt=RDHFLctoV0->GetVarsForOpt();
 238     }else{
 239       TString *names;
 240       names=new TString[nvars];
 241       TString answer="";
 242       Int_t checktrue=0;
 243       names=RDHFLctoV0->GetVarNames();
 244       for(Int_t i=0;i<nvars;i++){
 245         cout<<names[i]<<" for opt? (y/n)"<<endl;
 246         cin>>answer;
 247         if(answer=="y") {
 248           boolforopt[i]=kTRUE;
 249           checktrue++;
 250         }
 251         else boolforopt[i]=kFALSE;
 252       }
 253       if (checktrue!=dim) {
 254         cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
 255         return;
 256       }
 257       RDHFLctoV0->SetVarsForOpt(dim,boolforopt);
 258     }
 259   }
 260
 261
 262   Float_t tighterval[dim][nptbins];
 263   // 0(2): inv. mass V0 if K0S [GeV]
 264   // 1(3): inv. mass V0 if Lambda [GeV]
 265   // 2(4): pT min bachelor track [GeV/c]
 266   // 3(5): pT min V0-positive track [GeV/c]
 267   // 4(6): pT min V0-negative track [GeV/c]
 268   // 5(7): dca cascade cut (cm)
 269   // 6(8): dca V0 cut (cm)
 270
 271   // number of steps for each variable is set in the AddTask arguments (default=8)
 272   // set this!!
 273   for(Int_t ipt=0;ipt<nptbins;ipt++){
 274     tighterval[0][ipt]=0.075; // inv. mass V0 if K0S [GeV]
 275     tighterval[1][ipt]=0.040; // inv. mass V0 if Lambda [GeV]
 276     tighterval[2][ipt]=0.1; // pT min bachelor track [GeV/c]
 277     tighterval[3][ipt]=0.1; // pT min V0-positive track [GeV/c]
 278     tighterval[4][ipt]=0.1; // pT min V0-negative track [GeV/c]
 279     tighterval[5][ipt]=10.; // dca cascade cut (cm)
 280     tighterval[6][ipt]=10.; // dca v0 cut (cm)
 281   }
 282
 283   TString name="";
 284   Int_t arrdim=dim*nptbins;
 285   cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
 286   TClonesArray max("TParameter<float>",arrdim);
 287   for(Int_t ival=0;ival<dim;ival++){
 288     for(Int_t jpt=0;jpt<nptbins;jpt++){
 289       name=Form("par%dptbin%d",ival,jpt);
 290       cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
 291       new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
 292     }
 293   }
 294
 295   Bool_t flagPID=kFALSE;
 296   RDHFLctoV0->SetUsePID(flagPID);
 297
 298   RDHFLctoV0->PrintAll();
 299   printf("Use PID? %s\n",flagPID ? "yes" : "no");
 300
 301   //pid settings
 302   //1. bachelor: default one
 303   AliAODPidHF* pidObjBachelor = new AliAODPidHF();
 304   Double_t sigmasBac[5]={3.,1.,1.,3.,3.}; // 0, 1(A), 2(A) -> TPC; 3 -> TOF; 4 -> ITS
 305   pidObjBachelor->SetSigma(sigmasBac);
 306   pidObjBachelor->SetAsym(kFALSE);
 307   pidObjBachelor->SetMatch(1);
 308   pidObjBachelor->SetTPC(kTRUE);
 309   pidObjBachelor->SetTOF(kTRUE);
 310   pidObjBachelor->SetTOFdecide(kFALSE);
 311   RDHFLctoV0->SetPidHF(pidObjBachelor);
 312
 313   //2. V0pos
 314   AliAODPidHF* pidObjV0pos = new AliAODPidHF();
 315   Double_t sigmasV0pos[5]={3.,1.,1.,3.,3.}; // 0, 1(A), 2(A) -> TPC; 3 -> TOF; 4 -> ITS
 316   pidObjV0pos->SetSigma(sigmasV0pos);
 317   pidObjV0pos->SetAsym(kFALSE);
 318   pidObjV0pos->SetMatch(1);
 319   pidObjV0pos->SetTPC(kTRUE);
 320   pidObjV0pos->SetTOF(kTRUE);
 321   pidObjV0pos->SetTOFdecide(kFALSE);
 322   RDHFLctoV0->SetPidHF(pidObjV0pos);
 323
 324   //2. V0neg
 325   AliAODPidHF* pidObjV0neg = new AliAODPidHF();
 326   Double_t sigmasV0neg[5]={3.,1.,1.,3.,3.}; // 0, 1(A), 2(A) -> TPC; 3 -> TOF; 4 -> ITS
 327   pidObjV0neg->SetSigma(sigmasV0neg);
 328   pidObjV0neg->SetAsym(kFALSE);
 329   pidObjV0neg->SetMatch(1);
 330   pidObjV0neg->SetTPC(kTRUE);
 331   pidObjV0neg->SetTOF(kTRUE);
 332   pidObjV0neg->SetTOFdecide(kFALSE);
 333   RDHFLctoV0->SetPidHF(pidObjV0neg);
 334
 335   //activate pileup rejection (for pp)
 336   //RDHFLctoV0->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
 337
 338   //Do not recalculate the vertex
 339   RDHFLctoV0->SetRemoveDaughtersFromPrim(kFALSE); //activate for pp
 340
 341   TString cent="";
 342   //centrality selection (Pb-Pb)
 343   Float_t minc=20,maxc=80;
 344   RDHFLctoV0->SetMinCentrality(minc);
 345   RDHFLctoV0->SetMaxCentrality(maxc);
 346   cent=Form("%.0f%.0f",minc,maxc);
 347   //RDHFLctoV0->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
 348   RDHFLctoV0->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
 349
 350   //temporary
 351   //RDHFLctoV0->SetFixRefs();
 352
 353   TFile* fout=new TFile(Form("cuts4SignifMaxim%s%s%sRecVtx%sPileupRej.root",
 354                              RDHFLctoV0->GetUseCentrality()==0 ? "pp" : "PbPb",
 355                              cent.Data(),
 356                              RDHFLctoV0->GetIsPrimaryWithoutDaughters() ? "" : "No",
 357                              RDHFLctoV0->GetOptPileUp() ? "" : "No"),"recreate"); //set this!!
 358
 359   fout->cd();
 360   RDHFLctoV0->Write();
 361   max.Write();
 362   fout->Close();
 363
 364 }