PWG3/vertexingHF/macros/makeTFile4CutsLctopKpi.C

   1 #include <Riostream.h>
   2 #include <TFile.h>
   3 #include <AliRDHFCutsLctopKpi.h>
   4 #include <AliAODPidHF.h>
   5 #include <TClonesArray.h>
   6 #include <TParameter.h>
   7
   8
   9 //Use:
  10 //Set hard coded commentet with //set this!!
  11 // root[] .L makeInput...C++
  12 // root[] makeInputAliAnalysisTaskSE...()
  13 //similar macros for the other D mesons
  14
  15 //Author: Rosa Romita, r.romita@gsi.de
  16
  17
  18 //macro to make a .root file which contains an AliRDHFCutsLctopKpi for AliAnalysisTaskSELambdac task
  19
  20 void makeInputAliAnalysisTaskSELctopKpi(){
  21
  22   AliRDHFCutsLctopKpi* RDHFLctopKpiProd=new AliRDHFCutsLctopKpi();
  23   RDHFLctopKpiProd->SetName("LctopKpiProdCuts");
  24   RDHFLctopKpiProd->SetTitle("Production cuts for Lc analysis");
  25
  26   AliRDHFCutsLctopKpi* RDHFLctopKpiAn=new AliRDHFCutsLctopKpi();
  27   RDHFLctopKpiAn->SetName("LctopKpiAnalysisCuts");
  28   RDHFLctopKpiAn->SetTitle("Analysis cuts for Lc analysis");
  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   RDHFLctopKpiProd->AddTrackCuts(esdTrackCuts);
  45   RDHFLctopKpiAn->AddTrackCuts(esdTrackCuts);
  46
  47   const Int_t nvars=13;
  48
  49   const Int_t nptbins=4;
  50   Float_t* ptbins;
  51   ptbins=new Float_t[nptbins+1];
  52   ptbins[0]=0.;
  53   ptbins[1]=2.;
  54   ptbins[2]=3.;
  55   ptbins[3]=4.;
  56   ptbins[4]=12.;
  57
  58
  59   Float_t** prodcutsval;
  60   prodcutsval=new Float_t*[nvars];
  61   for(Int_t iv=0;iv<nvars;iv++){
  62     prodcutsval[iv]=new Float_t[nptbins];
  63   }
  64
  65   //  inv. mass [GeV]
  66   // pTK [GeV/c]
  67   // pTPi [GeV/c]
  68   // d0K [cm]   lower limit!
  69   // d0Pi [cm]  lower limit!
  70   // dist12 (cm)
  71   // sigmavert (cm)
  72   // dist prim-sec (cm)
  73   // pM=Max{pT1,pT2,pT3} (GeV/c)
  74   // cosThetaPoint
  75   // Sum d0^2 (cm^2)
  76   // dca cut (cm)
  77   for(Int_t ipt=0;ipt<nptbins;ipt++){
  78     prodcutsval[0][ipt]=0.18;
  79     prodcutsval[1][ipt]=0.4;
  80     prodcutsval[2][ipt]=0.5;
  81     prodcutsval[3][ipt]=0.;
  82     prodcutsval[4][ipt]=0.;
  83     prodcutsval[5][ipt]=0.01;
  84     prodcutsval[6][ipt]=0.06;
  85     prodcutsval[7][ipt]=0.005;
  86     prodcutsval[8][ipt]=0.7;
  87     prodcutsval[9][ipt]=0.;
  88     prodcutsval[10][ipt]=0.;
  89     prodcutsval[11][ipt]=0.05;
  90     prodcutsval[12][ipt]=0.4;
  91   }
  92
  93   RDHFLctopKpiProd->SetPtBins(nptbins+1,ptbins);
  94   RDHFLctopKpiProd->SetCuts(nvars,nptbins,prodcutsval);
  95
  96   Float_t** anacutsval;
  97   anacutsval=new Float_t*[nvars];
  98   for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
  99   for(Int_t ipt2=0;ipt2<nptbins;ipt2++){
 100    anacutsval[0][ipt2]=0.18;
 101    anacutsval[1][ipt2]=0.4;
 102    anacutsval[2][ipt2]=0.5;
 103    anacutsval[3][ipt2]=0.;
 104    anacutsval[4][ipt2]=0.;
 105    anacutsval[5][ipt2]=0.01;
 106    anacutsval[6][ipt2]=0.06;
 107    anacutsval[7][ipt2]=0.005;
 108    anacutsval[8][ipt2]=0.7;
 109    anacutsval[9][ipt2]=0.;
 110    anacutsval[10][ipt2]=0.;
 111    anacutsval[11][ipt2]=0.05;
 112    anacutsval[12][ipt2]=0.4;
 113   }
 114
 115   // pt kaon
 116   anacutsval[1][0]=0.5;
 117   anacutsval[1][1]=0.85;
 118   anacutsval[1][2]=0.9;
 119   anacutsval[1][3]=0.4;
 120   //pt proton
 121   anacutsval[2][0]=0.5;
 122   anacutsval[2][1]=0.6;
 123   anacutsval[2][2]=0.9;
 124   anacutsval[2][3]=0.9;
 125
 126   //pt pion
 127   anacutsval[12][0]=0.475;
 128   anacutsval[12][1]=0.75;
 129   anacutsval[12][2]=0.75;
 130   anacutsval[12][3]=0.7;
 131
 132   anacutsval[5][0]=0.02;
 133   anacutsval[5][1]=0.025;
 134   anacutsval[5][2]=0.02;
 135   anacutsval[5][3]=0.01;
 136
 137   anacutsval[7][0]=0.00625;
 138   anacutsval[7][1]=0.0125;
 139   anacutsval[7][2]=0.005;
 140   anacutsval[7][3]=0.007;
 141
 142   anacutsval[9][0]=0.5;
 143   anacutsval[9][1]=0.2;
 144   anacutsval[9][2]=0.6;
 145   anacutsval[9][3]=0.;
 146
 147   anacutsval[10][0]=0.00125;
 148
 149   RDHFLctopKpiAn->SetPtBins(nptbins+1,ptbins);
 150   RDHFLctopKpiAn->SetCuts(nvars,nptbins,anacutsval);
 151
 152  //  RDHFLc->SetRecoKF(); //set this if you want to recompute the secondary vertex with the KF package
 153     //pid settings
 154     //1. kaon: default one
 155   AliAODPidHF* pidObjK=new AliAODPidHF();
 156   Double_t sigmasK[5]={3.,1.,1.,3.,2.};
 157   pidObjK->SetSigma(sigmasK);
 158   pidObjK->SetAsym(kTRUE);
 159   pidObjK->SetMatch(1);
 160   pidObjK->SetTPC(kTRUE);
 161   pidObjK->SetTOF(kTRUE);
 162   pidObjK->SetITS(kTRUE);
 163   Double_t plimK[2]={0.5,0.8};
 164   pidObjK->SetPLimit(plimK,2);
 165   pidObjK->SetTOFdecide(kTRUE);
 166
 167   RDHFLctopKpiProd->SetPidHF(pidObjK);
 168   RDHFLctopKpiAn->SetPidHF(pidObjK);
 169
 170     //2. pion
 171   AliAODPidHF* pidObjpi=new AliAODPidHF();
 172   pidObjpi->SetTPC(kTRUE);
 173   Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
 174   pidObjpi->SetSigma(sigmaspi);
 175 //  pidObjpi->SetTOFdecide(kTRUE);
 176
 177   RDHFLctopKpiProd->SetPidpion(pidObjpi);
 178   RDHFLctopKpiAn->SetPidpion(pidObjpi);
 179
 180   // 3. proton
 181   AliAODPidHF* pidObjp=new AliAODPidHF();
 182   Double_t sigmasp[5]={3.,1.,1.,3.,2.};
 183   pidObjp->SetSigma(sigmasp);
 184   pidObjp->SetAsym(kTRUE);
 185   pidObjp->SetMatch(1);
 186   pidObjp->SetTPC(kTRUE);
 187   pidObjp->SetTOF(kTRUE);
 188   pidObjp->SetITS(kTRUE);
 189   Double_t plimp[2]={1.,2.};
 190   pidObjp->SetPLimit(plimp,2);
 191   pidObjp->SetTOFdecide(kTRUE);
 192
 193   RDHFLctopKpiProd->SetPidprot(pidObjp);
 194   RDHFLctopKpiAn->SetPidprot(pidObjp);
 195
 196   Bool_t pidflag=kTRUE;
 197   RDHFLctopKpiAn->SetUsePID(pidflag);
 198   RDHFLctopKpiProd->SetUsePID(pidflag);
 199   if(pidflag) cout<<"PID is used"<<endl;
 200   else cout<<"PID is not used"<<endl;
 201
 202   cout<<"This is the object I'm going to save:"<<endl;
 203   RDHFLctopKpiProd->PrintAll();
 204   RDHFLctopKpiAn->PrintAll();
 205   TFile* fout=new TFile("cuts4LctopKpi.root","RECREATE");
 206   fout->cd();
 207   RDHFLctopKpiProd->Write();
 208   RDHFLctopKpiAn->Write();
 209   fout->Close();
 210   delete fout;
 211   delete pidObjp;
 212   delete pidObjpi;
 213   delete pidObjK;
 214   delete RDHFLctopKpiProd;
 215   delete RDHFLctopKpiAn;
 216
 217 }
 218
 219
 220 //macro to make a .root file (for significance maximization) which contains an AliRDHFCutsLctopKpi with loose set of cuts  and TParameter with the tighest value of these cuts
 221
 222 void makeInputAliAnalysisTaskSESignificanceMaximization(){
 223
 224   AliRDHFCutsLctopKpi* RDHFLctopKpi=new AliRDHFCutsLctopKpi();
 225   RDHFLctopKpi->SetName("loosercuts");
 226   RDHFLctopKpi->SetTitle("Cuts for significance maximization");
 227
 228   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
 229   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 230   //default
 231   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 232   esdTrackCuts->SetMinNClustersTPC(70);
 233   esdTrackCuts->SetRequireITSRefit(kTRUE);
 234   esdTrackCuts->SetMinNClustersITS(4);
 235
 236   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 237   esdTrackCuts->SetMinDCAToVertexXY(0.);
 238   esdTrackCuts->SetEtaRange(-0.8,0.8);
 239   esdTrackCuts->SetPtRange(0.3,1.e10);
 240
 241   RDHFLctopKpi->AddTrackCuts(esdTrackCuts);
 242
 243   const Int_t nvars=13;
 244
 245   const Int_t nptbins=4; //change this when adding pt bins!
 246   Float_t ptbins[nptbins+1];
 247   ptbins[0]=0.;
 248   ptbins[1]=2.;
 249   ptbins[2]=3.;
 250   ptbins[3]=4.;
 251   ptbins[4]=12.;
 252
 253   RDHFLctopKpi->SetPtBins(nptbins+1,ptbins);
 254
 255   Float_t** rdcutsvalmine;
 256   rdcutsvalmine=new Float_t*[nvars];
 257   for(Int_t iv=0;iv<nvars;iv++){
 258     rdcutsvalmine[iv]=new Float_t[nptbins];
 259   }
 260
 261   //setting my cut values
 262   //  inv. mass [GeV]
 263   // pTK [GeV/c]
 264   // pTP [GeV/c]
 265   // d0K [cm]   lower limit!
 266   // d0Pi [cm]  lower limit!
 267   // dist12 (cm)
 268   // sigmavert (cm)
 269   // dist prim-sec (cm)
 270   // pM=Max{pT1,pT2,pT3} (GeV/c)
 271   // cosThetaPoint
 272   // Sum d0^2 (cm^2)
 273   // dca cut (cm)
 274   // pt pion
 275   Float_t cutsMatrixLctopKpiStand[nptbins][nvars]=
 276   {{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4},
 277    {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4},
 278    {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4},
 279    {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05,0.4}};
 280
 281   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
 282   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
 283   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
 284   for (Int_t ibin=0;ibin<nptbins;ibin++){
 285     for (Int_t ivar = 0; ivar<nvars; ivar++){
 286       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixLctopKpiStand[ibin][ivar];
 287     }
 288   }
 289   RDHFLctopKpi->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
 290
 291
 292   Int_t nvarsforopt=RDHFLctopKpi->GetNVarsForOpt();
 293   Int_t dim=4; //set this!!
 294   Bool_t *boolforopt;
 295   boolforopt=new Bool_t[nvars];
 296   if(dim>nvarsforopt){
 297     cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
 298     return;
 299   } else {
 300     if(dim==nvarsforopt){
 301       boolforopt=RDHFLctopKpi->GetVarsForOpt();
 302     }else{
 303       TString *names;
 304       names=new TString[nvars];
 305       TString answer="";
 306       Int_t checktrue=0;
 307       names=RDHFLctopKpi->GetVarNames();
 308       for(Int_t i=0;i<nvars;i++){
 309         cout<<names[i]<<" for opt? (y/n)"<<endl;
 310         cin>>answer;
 311         if(answer=="y") {
 312           boolforopt[i]=kTRUE;
 313           checktrue++;
 314         }
 315         else boolforopt[i]=kFALSE;
 316       }
 317       if (checktrue!=dim) {
 318         cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
 319         return;
 320       }
 321       RDHFLctopKpi->SetVarsForOpt(dim,boolforopt);
 322     }
 323   }
 324
 325
 326   Float_t tighterval[dim][nptbins];
 327   // 5 dist12 (cm)        <---
 328   // 7 dist prim-sec (cm) <---
 329   // 9 cosThetaPoint      <---
 330   // 10 Sum d0^2 (cm^2)   <---
 331   // 11 dca cut (cm)
 332   // {0.18,0.4,0.5,0.,0.,0.01 <--- ,0.06,0.005 <--- ,0.7,0. <--- ,0. <--- ,0.05 }
 333
 334   //number of steps for each variable is set in the AddTask arguments (default=8)
 335   //set this!!
 336   for(Int_t ipt=0;ipt<nptbins;ipt++){
 337     tighterval[0][ipt]=0.03;
 338     tighterval[1][ipt]=0.02;
 339     tighterval[2][ipt]=1.;
 340     tighterval[3][ipt]=0.01;
 341   }
 342
 343
 344
 345   TString name="";
 346   Int_t arrdim=dim*nptbins;
 347   cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
 348   TClonesArray max("TParameter<float>",arrdim);
 349   for(Int_t ival=0;ival<dim;ival++){
 350     for(Int_t jpt=0;jpt<nptbins;jpt++){
 351       name=Form("par%dptbin%d",ival,jpt);
 352       cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
 353       new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
 354     }
 355   }
 356
 357   Bool_t flagPID=kFALSE;
 358   RDHFLctopKpi->SetUsePID(flagPID);
 359
 360   RDHFLctopKpi->PrintAll();
 361   printf("Use PID? %s\n",flagPID ? "yes" : "no");
 362
 363     //pid settings
 364     //1. kaon: default one
 365   AliAODPidHF* pidObjK=new AliAODPidHF();
 366   Double_t sigmasK[5]={3.,1.,1.,3.,2.};
 367   pidObjK->SetSigma(sigmasK);
 368   pidObjK->SetAsym(kTRUE);
 369   pidObjK->SetMatch(1);
 370   pidObjK->SetTPC(kTRUE);
 371   pidObjK->SetTOF(kTRUE);
 372   pidObjK->SetITS(kTRUE);
 373   Double_t plimK[2]={0.5,0.8};
 374   pidObjK->SetPLimit(plimK,2);
 375   pidObjK->SetTOFdecide(kTRUE);
 376
 377   RDHFLctopKpi->SetPidHF(pidObjK);
 378
 379     //2. pion
 380   AliAODPidHF* pidObjpi=new AliAODPidHF();
 381   pidObjpi->SetTPC(kTRUE);
 382   Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
 383   pidObjpi->SetSigma(sigmaspi);
 384   pidObjpi->SetTOFdecide(kTRUE);
 385
 386   RDHFLctopKpi->SetPidpion(pidObjpi);
 387
 388   // 3. proton
 389   AliAODPidHF* pidObjp=new AliAODPidHF();
 390   Double_t sigmasp[5]={3.,1.,1.,3.,2.};
 391   pidObjp->SetSigma(sigmasp);
 392   pidObjp->SetAsym(kTRUE);
 393   pidObjp->SetMatch(1);
 394   pidObjp->SetTPC(kTRUE);
 395   pidObjp->SetTOF(kTRUE);
 396   pidObjp->SetITS(kTRUE);
 397   Double_t plimp[2]={1.,2.};
 398   pidObjp->SetPLimit(plimp,2);
 399   pidObjp->SetTOFdecide(kTRUE);
 400
 401   RDHFLctopKpi->SetPidprot(pidObjp);
 402
 403   //activate pileup rejection (for pp)
 404   //RDHFLctopKpi->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
 405
 406   //Do not recalculate the vertex
 407   RDHFLctopKpi->SetRemoveDaughtersFromPrim(kFALSE); //activate for pp
 408
 409   TString cent="";
 410   //centrality selection (Pb-Pb)
 411   Float_t minc=20,maxc=80;
 412   RDHFLctopKpi->SetMinCentrality(minc);
 413   RDHFLctopKpi->SetMaxCentrality(maxc);
 414   cent=Form("%.0f%.0f",minc,maxc);
 415   //  RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
 416   RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
 417
 418   //temporary
 419   //  RDHFLctopKpi->SetFixRefs();
 420
 421   TFile* fout=new TFile(Form("cuts4SignifMaxim%s%s%sRecVtx%sPileupRej.root", RDHFLctopKpi->GetUseCentrality()==0 ? "pp" : "PbPb",cent.Data(),RDHFLctopKpi->GetIsPrimaryWithoutDaughters() ? "" : "No",RDHFLctopKpi->GetOptPileUp() ? "" : "No"),"recreate");   //set this!!
 422
 423   fout->cd();
 424   RDHFLctopKpi->Write();
 425   max.Write();
 426   fout->Close();
 427
 428 }
 429