PWGHF/vertexingHF/macros/makeTFile4CutsDstoKKpi.C

   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
  54     const Int_t nvars=20;
  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.;
  78       anacutsval[16][ipt]=0.;
  79       anacutsval[17][ipt]=0.;
  80       anacutsval[18][ipt]=0.;
  81       anacutsval[19][ipt]=-1.;
  82
  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);
 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);
 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
 167     const Int_t nvars=20;
 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.;
 191       anacutsval[16][ipt]=0.;
 192       anacutsval[17][ipt]=0.;
 193       anacutsval[18][ipt]=0.;
 194       anacutsval[19][ipt]=-1.;
 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);
 245     // To be set only in case of strong pid
 246     //analysiscuts->SetMaxPtStrongPid(9999.);
 247     //analysiscuts->SetMaxPStrongPidK(1.5);
 248     //analysiscuts->SetMaxPStrongPidpi(1.5);
 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
 288   const Int_t nvars=20;
 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.;
 324       prodcutsval[16][ipt]=0.;
 325       prodcutsval[17][ipt]=0.;
 326       prodcutsval[18][ipt]=0.;
 327       prodcutsval[19][ipt]=-1.;
 328
 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 }