PWGHF/vertexingHF/macros/makeTFile4CutsDplustoKpipi.C

   1 #include <Riostream.h>
   2 #include <TFile.h>
   3 #include <AliRDHFCutsDplustoKpipi.h>
   4 #include <TClonesArray.h>
   5 #include <TParameter.h>
   6
   7 //macro to make a .root file which contains an AliRDHFCutsDplustoKpipi with loose set of cuts (for significance maximization) and TParameter with the tighest value of these cuts
   8 //Needed for AliAnalysisTaskSEDplus, AliCFTaskVertexingHF3Prong, AliAnalysisTaskSESignificance
   9
  10 //Use:
  11 //Set hard coded commented with //set this!!
  12
  13 //.L makeTFile4CutsDplustoKpipi.C
  14 // makeInputAliAnalysisTaskSEDplus()
  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 void makeInputAliAnalysisTaskSEDplusPP(){
  24
  25  //  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");
  26
  27
  28     AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
  29     esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
  30     //default
  31     esdTrackCuts->SetRequireTPCRefit(kTRUE);
  32     esdTrackCuts->SetRequireITSRefit(kTRUE);
  33     //esdTrackCuts->SetMinNClustersITS(4); // default is 5
  34     esdTrackCuts->SetMinNClustersTPC(70);
  35     esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
  36                                            AliESDtrackCuts::kAny);
  37     // default is kBoth, otherwise kAny
  38     esdTrackCuts->SetMinDCAToVertexXY(0.);
  39     esdTrackCuts->SetPtRange(0.3,1.e10);
  40
  41
  42
  43     const Int_t nptbins=14;
  44     Float_t* ptbins;
  45     ptbins=new Float_t[nptbins+1];
  46     ptbins[0]=0.;
  47     ptbins[1]=1.;
  48     ptbins[2]=2.;
  49     ptbins[3]=3.;
  50     ptbins[4]=4.;
  51     ptbins[5]=5.;
  52     ptbins[6]=6.;
  53     ptbins[7]=7.;
  54     ptbins[8]=8.;
  55     ptbins[9]=9.;
  56     ptbins[10]=10.;
  57     ptbins[11]=12.;
  58     ptbins[12]=14.;
  59     ptbins[13]=16.;
  60     ptbins[14]=24.;
  61     const Int_t nvars=14;
  62
  63
  64
  65     Float_t** anacutsval;
  66     anacutsval=new Float_t*[nvars];
  67
  68     for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
  69     //Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
  70     Int_t ic=0;
  71     for(Int_t ipt=0;ipt<nptbins;ipt++){
  72       anacutsval[ic][ipt]=0.2;
  73     }
  74     ic=1;
  75     for(Int_t ipt=3;ipt<nptbins;ipt++){
  76       anacutsval[ic][ipt]=0.4;
  77     }
  78     anacutsval[1][0]=0.3;
  79     anacutsval[1][1]=0.3;
  80     anacutsval[1][2]=0.4;
  81
  82     ic=2;
  83     for(Int_t ipt=3;ipt<nptbins;ipt++){
  84       anacutsval[ic][ipt]=0.4;
  85     }
  86
  87     anacutsval[2][0]=0.3;
  88     anacutsval[2][1]=0.3;
  89     anacutsval[2][2]=0.4;
  90
  91     ic=3;
  92     for(Int_t ipt=0;ipt<nptbins;ipt++){
  93       anacutsval[ic][ipt]=0.;
  94     }
  95
  96     ic=4;
  97     for(Int_t ipt=0;ipt<nptbins;ipt++){
  98       anacutsval[ic][ipt]=0.;
  99     }
 100
 101     ic=5;
 102     for(Int_t ipt=0;ipt<nptbins;ipt++){
 103       anacutsval[ic][ipt]=0.01;
 104     }
 105
 106
 107     ic=6;
 108     for(Int_t ipt=5;ipt<nptbins;ipt++){
 109       anacutsval[ic][ipt]=0.023333;
 110     }
 111
 112     anacutsval[6][0]=0.022100;
 113     anacutsval[6][1]=0.022100;
 114     anacutsval[6][2]=0.034;
 115     anacutsval[6][3]=0.020667;
 116     anacutsval[6][4]=0.020667;
 117
 118
 119     ic=7;
 120     for(Int_t ipt=5;ipt<nptbins;ipt++){
 121       anacutsval[ic][ipt]=0.115;
 122     }
 123
 124     anacutsval[7][0]=0.08;
 125     anacutsval[7][1]=0.08;
 126     anacutsval[7][2]=0.09;
 127     anacutsval[7][3]=0.095;
 128     anacutsval[7][4]=0.095;
 129
 130     ic=8;
 131     for(Int_t ipt=3;ipt<nptbins;ipt++){
 132       anacutsval[ic][ipt]=0.5;
 133     }
 134
 135     anacutsval[8][0]=0.5;
 136     anacutsval[8][1]=0.5;
 137     anacutsval[8][2]=1.0;
 138
 139
 140     ic=9;
 141     for(Int_t ipt=9;ipt<nptbins;ipt++){
 142       anacutsval[ic][ipt]=0.90;
 143     }
 144
 145     anacutsval[9][0]=0.95;
 146     anacutsval[9][1]=0.95;
 147     anacutsval[9][2]=0.95;
 148     anacutsval[9][3]=0.95;
 149     anacutsval[9][4]= 0.95;
 150     anacutsval[9][5]=0.92;
 151     anacutsval[9][6]=0.92;
 152     anacutsval[9][7]=0.92;
 153     anacutsval[9][8]=0.92;
 154
 155
 156     ic=10;
 157     for(Int_t ipt=3;ipt<nptbins;ipt++){
 158       anacutsval[ic][ipt]=0.000883;
 159     }
 160
 161     anacutsval[10][0]=0.0055;
 162     anacutsval[10][1]=0.0055;
 163     anacutsval[10][2]= 0.0028;
 164
 165
 166
 167
 168     ic=11;
 169     for(Int_t ipt=0;ipt<nptbins;ipt++){
 170       anacutsval[ic][ipt]=10000000000.;
 171     }
 172
 173     ic=12;
 174     for(Int_t ipt=0;ipt<nptbins;ipt++){
 175       anacutsval[ic][ipt]=0.;
 176     }
 177
 178     ic=13;
 179     for(Int_t ipt=0;ipt<nptbins;ipt++){
 180       anacutsval[ic][ipt]=0.;
 181     }
 182
 183
 184
 185
 186     AliRDHFCutsDplustoKpipi* analysiscuts=new AliRDHFCutsDplustoKpipi();
 187     analysiscuts->SetName("AnalysisCuts");
 188     analysiscuts->SetTitle("Cuts for Dplus Analysis and CF");
 189     analysiscuts->SetPtBins(nptbins+1,ptbins);
 190     analysiscuts->SetCuts(nvars,nptbins,anacutsval);
 191     analysiscuts->AddTrackCuts(esdTrackCuts);
 192     analysiscuts->SetUsePID(kTRUE);
 193     analysiscuts->SetUseImpParProdCorrCut(kFALSE);
 194     cout<<"This is the odject I'm going to save:"<<nptbins<<endl;
 195
 196     analysiscuts->PrintAll();
 197     TFile* fout=new TFile("DplustoKpipiCuts.root","recreate");
 198     fout->cd();
 199     analysiscuts->Write();
 200     fout->Close();
 201
 202
 203 }
 204
 205
 206
 207
 208
 209
 210 void makeInputAliAnalysisTaskSEDplusPbPb(){
 211
 212  //  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");
 213
 214  AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
 215   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 216   //default
 217     esdTrackCuts->SetRequireTPCRefit(kTRUE);
 218     esdTrackCuts->SetRequireITSRefit(kTRUE);
 219     //esdTrackCuts->SetMinNClustersITS(4); // default is 5
 220     esdTrackCuts->SetMinNClustersTPC(70);
 221     esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 222                                            AliESDtrackCuts::kAny);
 223     // default is kBoth, otherwise kAny
 224     esdTrackCuts->SetMinDCAToVertexXY(0.);
 225     esdTrackCuts->SetPtRange(0.8,1.e10);
 226
 227     TString cent="";
 228     //centrality selection (Pb-Pb)
 229     Float_t minc=0.,maxc=20.;
 230   const Int_t nptbins=12;
 231     Float_t* ptbins;
 232     ptbins=new Float_t[nptbins+1];
 233     ptbins[0]=0.;
 234     ptbins[1]=1.;
 235     ptbins[2]=2.;
 236     ptbins[3]=3.;
 237     ptbins[4]=4.;
 238     ptbins[5]=5.;
 239     ptbins[6]=6.;
 240     ptbins[7]=8.;
 241     ptbins[8]=9.;
 242     ptbins[9]=10.;
 243     ptbins[10]=12.;
 244     ptbins[11]=16.;
 245     ptbins[12]=24.;
 246     const Int_t nvars=14;
 247
 248
 249
 250
 251   Float_t** anacutsval;
 252     anacutsval=new Float_t*[nvars];
 253     for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
 254
 255
 256      Int_t ic=0;
 257     for(Int_t ipt=0;ipt<nptbins;ipt++){
 258       anacutsval[ic][ipt]=0.2;
 259     }
 260
 261
 262     ic=1;
 263
 264     for(Int_t ipt=7;ipt<nptbins;ipt++){
 265       anacutsval[1][ipt]=1.0;
 266     }
 267
 268     anacutsval[1][0]=0.8;
 269     anacutsval[1][1]=0.8;
 270     anacutsval[1][2]=0.8;
 271     anacutsval[1][3]=0.8;
 272     anacutsval[1][4]=0.9;
 273     anacutsval[1][5]=0.9;
 274     anacutsval[1][6]=0.8;
 275
 276
 277
 278
 279     ic=2;
 280     for(Int_t ipt=0;ipt<nptbins;ipt++){
 281       anacutsval[ic][ipt]=0.8;
 282     }
 283
 284
 285
 286
 287     ic=3;
 288     for(Int_t ipt=0;ipt<nptbins;ipt++){
 289       anacutsval[ic][ipt]=0.;
 290     }
 291     ic=4;
 292     for(Int_t ipt=0;ipt<nptbins;ipt++){
 293       anacutsval[ic][ipt]=0.;
 294     }
 295     ic=5;
 296     for(Int_t ipt=0;ipt<nptbins;ipt++){
 297       anacutsval[ic][ipt]=0.01;
 298     }
 299
 300
 301     ic=6;
 302     for(Int_t ipt=5;ipt<nptbins;ipt++){
 303      anacutsval[ic][ipt]=0.023333;
 304     }
 305      anacutsval[6][0]=0.022100;
 306      anacutsval[6][1]=0.022100;
 307      anacutsval[6][2]=0.034;
 308      anacutsval[6][3]=0.020667;
 309      anacutsval[6][4]=0.020667;
 310
 311
 312      ic=7;
 313     for(Int_t ipt=7;ipt<nptbins;ipt++){
 314      anacutsval[ic][ipt]=0.19;
 315     }
 316
 317     anacutsval[7][0]=0.08;
 318     anacutsval[7][1]=0.08;
 319     anacutsval[7][2]=0.17;
 320     anacutsval[7][3]=0.14;
 321     anacutsval[7][4]=0.14;
 322     anacutsval[7][5]=0.19;
 323     anacutsval[7][6]=0.14;
 324
 325
 326
 327
 328     ic=8;
 329     for(Int_t ipt=5;ipt<nptbins;ipt++){
 330       anacutsval[ic][ipt]=2.0;
 331     }
 332
 333     anacutsval[8][0]=0.8;
 334     anacutsval[8][1]=0.8;
 335     anacutsval[8][2]=1.1;
 336     anacutsval[8][3]=0.5;
 337     anacutsval[8][4]=0.5;
 338
 339
 340
 341     ic=9;
 342     for(Int_t ipt=7;ipt<nptbins;ipt++){
 343       anacutsval[ic][ipt]=0.997;
 344     }
 345
 346     anacutsval[9][0]=0.995;
 347     anacutsval[9][1]=0.995;
 348     anacutsval[9][2]=0.997;
 349     anacutsval[9][3]=0.998;
 350     anacutsval[9][4]= 0.998;
 351     anacutsval[9][5]=0.995;
 352     anacutsval[9][6]=0.995;
 353
 354
 355
 356
 357
 358
 359     ic=10;
 360     for(Int_t ipt=3;ipt<nptbins;ipt++){
 361       anacutsval[ic][ipt]=0.000883;
 362     }
 363
 364     anacutsval[10][0]=0.0055;
 365     anacutsval[10][1]=0.0055;
 366     anacutsval[10][2]= 0.0028;
 367
 368
 369
 370     ic=11;
 371     for(Int_t ipt=0;ipt<nptbins;ipt++){
 372       anacutsval[ic][ipt]=10000000000.;
 373     }
 374
 375
 376      ic=12;
 377     for(Int_t ipt=7;ipt<nptbins;ipt++){
 378       anacutsval[ic][ipt]=0.;
 379     }
 380
 381     anacutsval[12][0]=0.;
 382     anacutsval[12][1]=0.;
 383
 384     anacutsval[12][2]=0.;
 385     anacutsval[12][3]=12.;
 386     anacutsval[12][4]=12.;
 387     anacutsval[12][5]=12.0;
 388     anacutsval[12][6]=10.;
 389
 390
 391
 392   ic=13;
 393     for(Int_t ipt=7;ipt<nptbins;ipt++){
 394       anacutsval[ic][ipt]=0.;
 395     }
 396
 397     anacutsval[13][0]=0.;
 398     anacutsval[13][1]=0.;
 399
 400     anacutsval[13][2]=0.;
 401     anacutsval[13][3]=0.998571;
 402     anacutsval[13][4]=0.998571;
 403     anacutsval[13][5]=0.998571;
 404     anacutsval[13][6]=0.997143;
 405
 406
 407
 408
 409     AliRDHFCutsDplustoKpipi* analysiscuts=new AliRDHFCutsDplustoKpipi();
 410     analysiscuts->SetName("AnalysisCuts");
 411     analysiscuts->SetTitle("Cuts for Dplus Analysis and CF");
 412     analysiscuts->SetPtBins(nptbins+1,ptbins);
 413     analysiscuts->SetCuts(nvars,nptbins,anacutsval);
 414     analysiscuts->AddTrackCuts(esdTrackCuts);
 415
 416     analysiscuts->SetUsePID(kTRUE);
 417     analysiscuts->SetUseImpParProdCorrCut(kFALSE);
 418     analysiscuts->SetOptPileup(kFALSE);
 419     analysiscuts->SetUseAOD049(kTRUE);
 420     analysiscuts->SetMinCentrality(minc);
 421     analysiscuts->SetMaxCentrality(maxc);
 422     cent=Form("%.0f%.0f",minc,maxc);
 423     analysiscuts->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
 424     analysiscuts->SetMinPtCandidate(3.);
 425     analysiscuts->SetMaxPtCandidate(10000.);
 426     cout<<"This is the odject I'm going to save:"<<nptbins<<endl;
 427
 428     analysiscuts->PrintAll();
 429     TFile* fout=new TFile("DplustoKpipiCuts.root","recreate");
 430     fout->cd();
 431     analysiscuts->Write();
 432     fout->Close();
 433
 434 }
 435
 436
 437
 438 void makeInputAliAnalysisTaskSESignificanceMaximization(){
 439
 440   AliRDHFCutsDplustoKpipi* RDHFDplustoKpipi=new AliRDHFCutsDplustoKpipi();
 441   RDHFDplustoKpipi->SetName("loosercuts");
 442   RDHFDplustoKpipi->SetTitle("Cuts for significance maximization");
 443
 444   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
 445   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 446   //default
 447   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 448   esdTrackCuts->SetRequireITSRefit(kTRUE);
 449   //esdTrackCuts->SetMinNClustersITS(4); // default is 5
 450   esdTrackCuts->SetMinNClustersTPC(70);
 451   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 452                                          AliESDtrackCuts::kAny);
 453   // default is kBoth, otherwise kAny
 454   esdTrackCuts->SetMinDCAToVertexXY(0.);
 455   esdTrackCuts->SetPtRange(0.3,1.e10);
 456
 457   RDHFDplustoKpipi->AddTrackCuts(esdTrackCuts);
 458
 459   const Int_t nvars=12;
 460
 461   const Int_t nptbins=4;
 462   Float_t* ptbins;
 463   ptbins=new Float_t[nptbins+1];
 464   ptbins[0]=0.;
 465   ptbins[1]=2.;
 466   ptbins[2]=3.;
 467   ptbins[3]=5.;
 468   ptbins[4]=99999.;
 469
 470   RDHFDplustoKpipi->SetPtBins(nptbins+1,ptbins);
 471
 472   //setting my cut values
 473
 474   Float_t** prodcutsval;
 475   prodcutsval=new Float_t*[nvars];
 476   for(Int_t ic=0;ic<nvars;ic++){prodcutsval[ic]=new Float_t[nptbins];}
 477   for(Int_t ipt=0;ipt<nptbins;ipt++){
 478     prodcutsval[0][ipt]=0.2;
 479     prodcutsval[1][ipt]=0.4;
 480     prodcutsval[2][ipt]=0.4;
 481     prodcutsval[3][ipt]=0.;
 482     prodcutsval[4][ipt]=0.;
 483     prodcutsval[5][ipt]=0.01;
 484     prodcutsval[6][ipt]=0.06;
 485     prodcutsval[7][ipt]=0.02;
 486     prodcutsval[8][ipt]=0.;
 487     prodcutsval[9][ipt]=0.85;
 488     prodcutsval[10][ipt]=0.;
 489     prodcutsval[11][ipt]=10000000.0;
 490
 491   }
 492
 493   RDHFDplustoKpipi->SetCuts(nvars,nptbins,prodcutsval);
 494
 495   Int_t nvarsforopt=RDHFDplustoKpipi->GetNVarsForOpt();
 496   const Int_t dim=5; //set this!!
 497   Bool_t *boolforopt;
 498   boolforopt=new Bool_t[nvars];
 499   if(dim>nvarsforopt){
 500     cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
 501     return;
 502   } else {
 503     if(dim==nvarsforopt){
 504       boolforopt=RDHFDplustoKpipi->GetVarsForOpt();
 505     }else{
 506       TString *names;
 507       names=new TString[nvars];
 508       TString answer="";
 509       Int_t checktrue=0;
 510       names=RDHFDplustoKpipi->GetVarNames();
 511       for(Int_t i=0;i<nvars;i++){
 512         cout<<names[i]<<" for opt? (y/n)"<<endl;
 513         cin>>answer;
 514         if(answer=="y") {
 515           boolforopt[i]=kTRUE;
 516           checktrue++;
 517         }
 518         else boolforopt[i]=kFALSE;
 519       }
 520       if (checktrue!=dim) {
 521         cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
 522         return;
 523       }
 524       RDHFDplustoKpipi->SetVarsForOpt(dim,boolforopt);
 525     }
 526   }
 527
 528
 529   Float_t tighterval[dim][nptbins];
 530   //sigmavert
 531   //declength
 532   //Pmax
 533   //costhetapoint
 534   //sumd02
 535
 536   //number of steps for each variable is 4 now
 537   //set this!!
 538
 539   // tighterval[var][ptbin]
 540   tighterval[0][0]=0.022100;
 541   tighterval[0][1]=0.034;
 542   tighterval[0][2]=0.020667;
 543   tighterval[0][3]=0.023333;
 544
 545   tighterval[1][0]=0.08;
 546   tighterval[1][1]=0.09;
 547   tighterval[1][2]=0.095;
 548   tighterval[1][3]=0.115;
 549
 550   tighterval[2][0]=1.;
 551   tighterval[2][1]=1.;
 552   tighterval[2][2]=1.;
 553   tighterval[2][3]=1.;
 554
 555   tighterval[3][0]=0.979;
 556   tighterval[3][1]=0.9975;
 557   tighterval[3][2]=0.995;
 558   tighterval[3][3]=0.9975;
 559
 560   tighterval[4][0]=0.0055;
 561   tighterval[4][1]=0.0028;
 562   tighterval[4][2]=0.000883;
 563   tighterval[4][3]=0.000883;
 564
 565   TString name="";
 566   Int_t arrdim=dim*nptbins;
 567   TClonesArray max("TParameter<float>",arrdim);
 568   for (Int_t ipt=0;ipt<nptbins;ipt++){
 569     for(Int_t ival=0;ival<dim;ival++){
 570       name=Form("par%dptbin%d",ival,ipt);
 571       new(max[ipt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][ipt]);
 572     }
 573   }
 574
 575   TFile* fout=new TFile("cuts4SignifMaximDplus.root","recreate");   //set this!!
 576   fout->cd();
 577   RDHFDplustoKpipi->Write();
 578   max.Write();
 579   fout->Close();
 580
 581 }