PWG3/vertexingHF/macros/MakeCuts4Charm4Prong.C

   1 //gSystem->AddIncludePath("-I/$ALICE_ROOT/PWG3/vertexingHF");
   2
   3 #include <Riostream.h>
   4 #include <TFile.h>
   5 //#include <AliRDHFCutsD0toKpipipi.h>
   6 #include <TClonesArray.h>
   7 #include <TParameter.h>
   8
   9 //Author: Chiara Bianchin, cbianchi@pd.infn.it, Fabio Colamaria, fabio.colamaria@ba.infn.it
  10
  11 //macro to make a .root file which contains an AliRDHFCutsD0toKpipipi for AliAnalysisTaskSESelectHF4Prong task
  12
  13 void MakeCuts4Charm4Prong(){
  14
  15 gSystem->Load("libTree.so");
  16 gSystem->Load("libGeom.so");
  17 gSystem->Load("libPhysics.so");
  18 gSystem->Load("libVMC.so");
  19 gSystem->Load("libMinuit.so");
  20 gSystem->Load("libSTEERBase.so");
  21 gSystem->Load("libESD.so");
  22 gSystem->Load("libAOD.so");
  23 gSystem->Load("libANALYSIS.so");
  24 gSystem->Load("libANALYSISalice.so");
  25 gSystem->Load("libCORRFW.so");
  26 gSystem->Load("libPWG3base.so");
  27 gSystem->Load("libPWG3vertexingHF.so");
  28 gSystem->Load("libPWG3muon.so");
  29
  30   AliRDHFCutsD0toKpipipi* RDHFCharm4Prong=new AliRDHFCutsD0toKpipipi();
  31   RDHFCharm4Prong->SetName("Charm4ProngCuts");
  32   RDHFCharm4Prong->SetTitle("Cuts for D0 analysis");
  33
  34   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
  35   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
  36   //default
  37   esdTrackCuts->SetRequireTPCRefit(kTRUE);
  38   esdTrackCuts->SetRequireITSRefit(kTRUE);
  39   esdTrackCuts->SetMinNClustersITS(4); // default is 5
  40   //esdTrackCuts->SetMinNClustersTPC(70);
  41   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
  42                                          AliESDtrackCuts::kAny);
  43  // default is kBoth, otherwise kAny
  44   esdTrackCuts->SetMinDCAToVertexXY(0.);
  45   esdTrackCuts->SetPtRange(0.3,1.e10);
  46
  47   RDHFCharm4Prong->AddTrackCuts(esdTrackCuts);
  48
  49
  50
  51   const Int_t nvars=9;
  52
  53   const Int_t nptbins=5;
  54   Float_t* ptbins;
  55   ptbins=new Float_t[nptbins+1];
  56   ptbins[0]=0.;
  57   ptbins[1]=4.5;
  58   ptbins[2]=6.;
  59   ptbins[3]=8.;
  60   ptbins[4]=12.;
  61   ptbins[5]=16.;
  62
  63   RDHFCharm4Prong->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]");
  75     //     printf("    DCA [cm]");
  76     //     printf("    Dist 2-trk Vtx to PrimVtx [cm]");
  77     //     printf("    Dist 3-trk Vtx to PrimVtx [cm]");
  78     //     printf("    Dist 4-trk Vtx to PrimVtx [cm]");
  79     //     printf("    cosThetaPoint");
  80     //     printf("    pt [GeV/c]");
  81     //     printf("    rho mass");
  82     //     printf("    PID cut");
  83
  84   //setting my cut values
  85   rdcutsvalmine[0][0]=0.2;
  86   rdcutsvalmine[1][0]=0.025;
  87   rdcutsvalmine[2][0]=0.0900;
  88   rdcutsvalmine[3][0]=0.0900;
  89   rdcutsvalmine[4][0]=0.0600;
  90   rdcutsvalmine[5][0]=0.995;
  91   rdcutsvalmine[6][0]=0.;
  92   rdcutsvalmine[7][0]=0.1;
  93   rdcutsvalmine[8][0]=0.;
  94
  95   rdcutsvalmine[0][1]=0.2;
  96   rdcutsvalmine[1][1]=0.015;
  97   rdcutsvalmine[2][1]=0.1200;
  98   rdcutsvalmine[3][1]=0.1000;
  99   rdcutsvalmine[4][1]=0.0850;
 100   rdcutsvalmine[5][1]=0.99;
 101   rdcutsvalmine[6][1]=0.;
 102   rdcutsvalmine[7][1]=0.1;
 103   rdcutsvalmine[8][1]=0.;
 104
 105   rdcutsvalmine[0][2]=0.2;
 106   rdcutsvalmine[1][2]=0.025;
 107   rdcutsvalmine[2][2]=0.0975;
 108   rdcutsvalmine[3][2]=0.0900;
 109   rdcutsvalmine[4][2]=0.0800;
 110   rdcutsvalmine[5][2]=0.99;
 111   rdcutsvalmine[6][2]=0.;
 112   rdcutsvalmine[7][2]=0.1;
 113   rdcutsvalmine[8][2]=0.;
 114
 115   rdcutsvalmine[0][3]=0.2;
 116   rdcutsvalmine[1][3]=0.015;
 117   rdcutsvalmine[2][3]=0.0975;
 118   rdcutsvalmine[3][3]=0.0700;
 119   rdcutsvalmine[4][3]=0.0750;
 120   rdcutsvalmine[5][3]=0.9825;
 121   rdcutsvalmine[6][3]=0.;
 122   rdcutsvalmine[7][3]=0.1;
 123   rdcutsvalmine[8][3]=0.;
 124
 125   rdcutsvalmine[0][4]=rdcutsvalmine[0][3];
 126   rdcutsvalmine[1][4]=rdcutsvalmine[1][3];
 127   rdcutsvalmine[2][4]=rdcutsvalmine[2][3];
 128   rdcutsvalmine[3][4]=rdcutsvalmine[3][3];
 129   rdcutsvalmine[4][4]=rdcutsvalmine[4][3];
 130   rdcutsvalmine[5][4]=rdcutsvalmine[5][3];
 131   rdcutsvalmine[6][4]=rdcutsvalmine[6][3];
 132   rdcutsvalmine[7][4]=rdcutsvalmine[7][3];
 133   rdcutsvalmine[8][4]=rdcutsvalmine[8][3];
 134
 135   RDHFCharm4Prong->SetCuts(nvars,nptbins,rdcutsvalmine);
 136   cout<<"This is the odject I'm going to save:"<<endl;
 137   RDHFCharm4Prong->PrintAll();
 138   TFile* fout=new TFile("Charm4ProngCutsDef_16GeV.root","recreate");   //set this!!
 139   fout->cd();
 140   RDHFCharm4Prong->Write();
 141   fout->Close();
 142
 143 }
 144
 145 void MakeCuts4Charm4ProngForMaxim(){
 146
 147 gSystem->Clear();
 148 gSystem->Load("libTree.so");
 149 gSystem->Load("libGeom.so");
 150 gSystem->Load("libPhysics.so");
 151 gSystem->Load("libVMC.so");
 152 gSystem->Load("libMinuit.so");
 153 gSystem->Load("libSTEERBase.so");
 154 gSystem->Load("libESD.so");
 155 gSystem->Load("libAOD.so");
 156 gSystem->Load("libANALYSIS.so");
 157 gSystem->Load("libANALYSISalice.so");
 158 gSystem->Load("libCORRFW.so");
 159 gSystem->Load("libPWG3base.so");
 160 gSystem->Load("libPWG3vertexingHF.so");
 161 gSystem->Load("libPWG3muon.so");
 162
 163   AliRDHFCutsD0toKpipipi* RDHFCharm4Prong=new AliRDHFCutsD0toKpipipi();
 164   RDHFCharm4Prong->SetName("loosercuts");
 165   RDHFCharm4Prong->SetTitle("Cuts for significance maximization");
 166
 167   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
 168   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 169   //default
 170   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 171   esdTrackCuts->SetRequireITSRefit(kTRUE);
 172   //esdTrackCuts->SetMinNClustersITS(4);
 173   //esdTrackCuts->SetMinNClustersTPC(120);
 174
 175   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 176   esdTrackCuts->SetMinDCAToVertexXY(0.);
 177   esdTrackCuts->SetEtaRange(-0.8,0.8);
 178   esdTrackCuts->SetPtRange(0.3,1.e10);
 179
 180   RDHFCharm4Prong->AddTrackCuts(esdTrackCuts);
 181
 182   const Int_t nvars=9;
 183
 184   const Int_t nptbins=5; //change this when adding pt bins!
 185   Float_t ptbins[nptbins+1];
 186   ptbins[0]=0.;
 187   ptbins[1]=4.5;
 188   ptbins[2]=6.;
 189   ptbins[3]=8.;
 190   ptbins[4]=12.;
 191   ptbins[5]=16.; //o 25!
 192
 193   RDHFCharm4Prong->SetPtBins(nptbins+1,ptbins);
 194
 195   Float_t** rdcutsvalmine;
 196   rdcutsvalmine=new Float_t*[nvars];
 197   for(Int_t iv=0;iv<nvars;iv++){
 198     rdcutsvalmine[iv]=new Float_t[nptbins];
 199   }
 200
 201   //setting my cut values
 202     //cuts order
 203     //     printf("    |M-MD0| [GeV]");
 204     //     printf("    DCA [cm]");
 205     //     printf("    Dist 2-trk Vtx to PrimVtx [cm]");
 206     //     printf("    Dist 3-trk Vtx to PrimVtx [cm]");
 207     //     printf("    Dist 4-trk Vtx to PrimVtx [cm]");
 208     //     printf("    cosThetaPoint");
 209     //     printf("    pt [GeV/c]");
 210     //     printf("    rho mass");
 211     //     printf("    PID cut");
 212
 213   Float_t cutsMatrixD0toKpipipiStand[nptbins][nvars]={{0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
 214
 215 {0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
 216
 217 {0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
 218
 219 {0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
 220
 221 {0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.}};
 222
 223   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
 224   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
 225   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
 226   for (Int_t ibin=0;ibin<nptbins;ibin++){
 227     for (Int_t ivar = 0; ivar<nvars; ivar++){
 228       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpipipiStand[ibin][ivar];
 229     }
 230   }
 231   RDHFCharm4Prong->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
 232
 233   Int_t nvarsforopt=RDHFCharm4Prong->GetNVarsForOpt();
 234   Int_t dim=5; //set this!!
 235   Bool_t *boolforopt;
 236   boolforopt=new Bool_t[nvars];
 237   if(dim>nvarsforopt){
 238     cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
 239     return;
 240   } else {
 241     if(dim==nvarsforopt){
 242       boolforopt=RDHFCharm4Prong->GetVarsForOpt();
 243     }else{
 244       TString *names;
 245       names=new TString[nvars];
 246       TString answer="";
 247       Int_t checktrue=0;
 248       names=RDHFCharm4Prong->GetVarNames();
 249       for(Int_t i=0;i<nvars;i++){
 250         cout<<names[i]<<" for opt? (y/n)"<<endl;
 251         cin>>answer;
 252         if(answer=="y") {
 253           boolforopt[i]=kTRUE;
 254           checktrue++;
 255         }
 256         else boolforopt[i]=kFALSE;
 257       }
 258       if (checktrue!=dim) {
 259         cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
 260         return;
 261       }
 262       RDHFCharm4Prong->SetVarsForOpt(dim,boolforopt);
 263     }
 264   }
 265
 266 // Tight values for variables to be maximized
 267 // (DCA, Dist12, Dist3, Dist4, CosThetaPoint)
 268 // Indexes: variable, ptbin
 269
 270   Float_t tighterval[5][5];
 271
 272   //number of steps for each variable is set in the AddTask arguments (default=8)
 273
 274   tighterval[0][0]=0.025;
 275   tighterval[1][0]=0.09;
 276   tighterval[2][0]=0.09;
 277   tighterval[3][0]=0.06;
 278   tighterval[4][0]=0.955;
 279
 280   tighterval[0][1]=0.015;
 281   tighterval[1][1]=0.12;
 282   tighterval[2][1]=0.10;
 283   tighterval[3][1]=0.085;
 284   tighterval[4][1]=0.99;
 285
 286   tighterval[0][2]=0.025;
 287   tighterval[1][2]=0.0975;
 288   tighterval[2][2]=0.09;
 289   tighterval[3][2]=0.08;
 290   tighterval[4][2]=0.99;
 291
 292   tighterval[0][3]=0.015;
 293   tighterval[1][3]=0.0975;
 294   tighterval[2][3]=0.07;
 295   tighterval[3][3]=0.075;
 296   tighterval[4][3]=0.9825;
 297
 298   tighterval[0][4]=0.015;
 299   tighterval[1][4]=0.0975;
 300   tighterval[2][4]=0.07;
 301   tighterval[3][4]=0.075;
 302   tighterval[4][4]=0.9825;
 303
 304   TString name="";
 305   Int_t arrdim=dim*nptbins;
 306   cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
 307   TClonesArray max("TParameter<float>",arrdim);
 308   for(Int_t ival=0;ival<dim;ival++){
 309     for(Int_t jpt=0;jpt<nptbins;jpt++){
 310       name=Form("par%dptbin%d",ival,jpt);
 311       cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
 312       new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
 313     }
 314   }
 315
 316   Bool_t flagPID=kFALSE;
 317   RDHFCharm4Prong->SetUsePID(flagPID);
 318
 319   RDHFCharm4Prong->PrintAll();
 320   printf("Use PID? %s\n",flagPID ? "yes" : "no");
 321
 322 /*
 323   //pid settings
 324   AliAODPidHF* pidObj=new AliAODPidHF();
 325   //pidObj->SetName("pid4D0");
 326   Int_t mode=1;
 327   const Int_t nlims=2;
 328   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
 329   Bool_t compat=kTRUE; //effective only for this mode
 330   Bool_t asym=kTRUE;
 331   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
 332   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
 333   pidObj->SetMatch(mode);
 334   pidObj->SetPLimit(plims,nlims);
 335   pidObj->SetSigma(sigmas);
 336   pidObj->SetCompat(compat);
 337   pidObj->SetTPC(kTRUE);
 338   pidObj->SetTOF(kTRUE);
 339   RDHFCharm4Prong->SetPidHF(pidObj);
 340
 341   RDHFCharm4Prong->SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
 342 */
 343
 344   //activate pileup rejection (for pp)
 345   RDHFCharm4Prong->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
 346
 347   //Do not recalculate the vertex
 348   RDHFCharm4Prong->SetRemoveDaughtersFromPrim(kTRUE); //activate for pp
 349
 350 /* Per il pp levo la centralità!
 351 //Metto kCentOff! Per il Pb c'era kCentV0M
 352   TString cent="";
 353   //centrality selection (Pb-Pb)
 354   Float_t minc=20,maxc=80;
 355   RDHFCharm4Prong->SetMinCentrality(minc);
 356   RDHFCharm4Prong->SetMaxCentrality(maxc);
 357   cent=Form("%.0f%.0f",minc,maxc);
 358   RDHFCharm4Prong->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
 359 */
 360
 361   //temporary
 362 //  RDHFCharm4Prong->SetFixRefs();
 363
 364   TFile* fout=new TFile("Charm4ProngCutsForMaxim.root","recreate");
 365   fout->cd();
 366   RDHFCharm4Prong->Write();
 367   max.Write();
 368   fout->Close();
 369 }