PWGHF/vertexingHF/macros/AddTaskCFVertexingHFLctoV0bachelorTMVA.C

   1 const Double_t ymin  = -1.2 ;
   2 const Double_t ymax  =  1.2 ;
   3 const Double_t cosminTS = -1.05;
   4 const Double_t cosmaxTS =  1.05;
   5 const Double_t cosmin = 0.7;
   6 const Double_t cosmax =  1.02;
   7 const Double_t cTmin = 0;  // micron
   8 const Double_t cTmax = 300;  // micron
   9 const Double_t dcamin = 0;  // micron
  10 const Double_t dcamax = 600;  // micron
  11 const Double_t d0xd0min = -80000;  // micron
  12 const Double_t d0xd0max = 100000;  // micron
  13 const Double_t phimin = 0.0;
  14 const Int_t    mintrackrefsTPC = 2 ;
  15 const Int_t    mintrackrefsITS = 3 ;
  16 const Int_t    charge  = 1 ;
  17 const Int_t    minclustersTPC = 50 ;
  18 // cuts
  19 const Double_t ptmin = 0.1;
  20 const Double_t ptmax = 9999.;
  21 const Double_t etamin = -0.9;
  22 const Double_t etamax = 0.9;
  23 const Double_t zmin = -15;
  24 const Double_t zmax = 15;
  25 const Int_t    minITSClusters = 5;
  26
  27 const Float_t centmin_0_10 = 0.;
  28 const Float_t centmax_0_10 = 10.;
  29 const Float_t centmin_10_60 = 10.;
  30 const Float_t centmax_10_60 = 60.;
  31 const Float_t centmin_60_100 = 60.;
  32 const Float_t centmax_60_100 = 100.;
  33 const Float_t centmax = 100.;
  34 const Float_t fakemin = -0.5;
  35 const Float_t fakemax = 2.5.;
  36 const Float_t cosminXY = 0.95;
  37 const Float_t cosmaxXY = 1.0;
  38 const Float_t normDecLXYmin = 0;
  39 const Float_t normDecLXYmax = 20;
  40 const Float_t multmin_0_20 = 0;
  41 const Float_t multmax_0_20 = 20;
  42 const Float_t multmin_20_50 = 20;
  43 const Float_t multmax_20_50 = 50;
  44 const Float_t multmin_50_80 = 50;
  45 const Float_t multmax_50_80 = 80;
  46 const Float_t multmin_80_100 = 80;
  47 const Float_t multmax_80_100 = 100;
  48 const Float_t multmin_100_400 = 100;
  49 const Float_t multmax_100_400 = 400;
  50
  51
  52 //----------------------------------------------------
  53
  54 AliCFTaskVertexingHF *AddTaskCFVertexingHFLctoV0bachelorTMVA(const char* cutFile = "CutsLc2pK0STMVA_20140311.root",
  55                                                              TString cutObjectName = "LctoV0AnalysisCuts",
  56                                                              TString suffix = "TMVA",
  57                                                              Int_t configuration = AliCFTaskVertexingHF::kCheetah,
  58                                                              Bool_t isKeepDfromB = kTRUE,
  59                                                              Bool_t isKeepDfromBOnly = kFALSE,
  60                                                              Int_t pdgCode = 4122,
  61                                                              Char_t isSign = 2,
  62                                                              Bool_t useWeight = kFALSE,
  63                                                              Bool_t useFlatPtWeight = kFALSE,
  64                                                              Bool_t useZWeight = kFALSE,
  65                                                              Bool_t useNchWeight = kFALSE,
  66                                                              Bool_t useNtrkWeight = kFALSE,
  67                                                              TString estimatorFilename="",
  68                                                              Int_t multiplicityEstimator = AliCFTaskVertexingHF::kNtrk10,
  69                                                              Bool_t isPPData = kFALSE,
  70                                                              Bool_t isPPbData = kFALSE,
  71                                                              Double_t refMult = 9.26,
  72                                                              Bool_t isFineNtrkBin = kFALSE)
  73 {
  74   printf("Adding CF task using cuts from file %s\n",cutFile);
  75   if (configuration == AliCFTaskVertexingHF::kSnail){
  76     printf("The configuration is set to be SLOW --> all the variables will be used to fill the CF\n");
  77   }
  78   else if (configuration == AliCFTaskVertexingHF::kCheetah){
  79     printf("The configuration is set to be FAST --> using only pt, y, ct, phi, zvtx, centrality, fake, multiplicity to fill the CF\n");
  80   }
  81   else{
  82     printf("The configuration is not defined! returning\n");
  83     return;
  84   }
  85
  86   gSystem->Sleep(2000);
  87
  88   // isSign = 0 --> D* only
  89   // isSign = 1 --> D*bar only
  90   // isSign = 2 --> D* + D*bar
  91
  92   TString expected;
  93   if (isSign == 0 && pdgCode < 0){
  94     AliError(Form("Error setting PDG code (%d) and sign (0 --> D* only): they are not compatible, returning"));
  95     return 0x0;
  96   }
  97   else if (isSign == 1 && pdgCode > 0){
  98     AliError(Form("Error setting PDG code (%d) and sign (1 --> D*bar only): they are not compatible, returning"));
  99     return 0x0;
 100   }
 101   else if (isSign > 2 || isSign < 0){
 102     AliError(Form("Sign not valid (%d, possible values are 0, 1, 2), returning"));
 103     return 0x0;
 104   }
 105
 106   TFile* fileCuts = TFile::Open(cutFile);
 107   if(!fileCuts || (fileCuts && !fileCuts->IsOpen())){
 108     AliError("Wrong cut file");
 109     return 0x0;
 110   }
 111
 112   AliRDHFCutsLctoV0 *cutsLcTMVA = (AliRDHFCutsLctoV0*)fileCuts->Get(cutObjectName.Data());
 113
 114   // check that the fKeepD0fromB flag is set to true when the fKeepD0fromBOnly flag is true
 115   //  for now the binning is the same than for all D's
 116   if (isKeepDfromBOnly) isKeepDfromB = true;
 117
 118   Double_t ptmin_0_6;
 119   Double_t ptmax_0_6;
 120   Double_t ptmin_6_8;
 121   Double_t ptmax_6_8;
 122   Double_t ptmin_8_16;
 123   Double_t ptmax_8_16;
 124   Double_t ptmin_16_24;
 125   Double_t ptmax_16_24;
 126
 127   ptmin_0_6 =  0.0 ;
 128   ptmax_0_6 =  6.0 ;
 129   ptmin_6_8 =  6.0 ;
 130   ptmax_6_8 =  8.0 ;
 131   ptmin_8_16 =  8.0 ;
 132   ptmax_8_16 =  16.0 ;
 133   ptmin_16_24 =  16.0 ;
 134   ptmax_16_24 =  24.0 ;
 135
 136   //CONTAINER DEFINITION
 137   Info("AliCFTaskVertexingHF","SETUP CONTAINER");
 138   const Double_t phimax = 2*TMath::Pi();
 139   UInt_t nstep = 10; //number of selection steps: MC with limited acceptance, MC, Acceptance, Vertex, Refit, Reco (no cuts), RecoAcceptance, RecoITSClusters (RecoAcceptance included), RecoPPR (RecoAcceptance+RecoITSCluster included), RecoPID
 140
 141   //const UInt_t ipT, iy, icosThetaStar, ipTpi, ipTk, icT, idca, id0xd0, ipointing, iphi, izvtx, icent, ifake, ipointingXY, iNormDecayLXY, imult;
 142   const Int_t nbiny  = 24 ; //bins in y
 143   const Int_t nbincosThetaStar  = 42 ; //bins in cosThetaStar
 144   const Int_t nbincT  = 15 ; //bins in cT
 145   const Int_t nbindca  = 20 ; //bins in dca
 146   const Int_t nbind0xd0  = 90 ; //bins in d0xd0
 147   const Int_t nbinpointing  = 50 ; //bins in cosPointingAngle
 148   const Int_t nbinphi  = 18 ; //bins in Phi
 149   const Int_t nbinzvtx  = 30 ; //bins in z vertex
 150   const Int_t nbincent = 28;  //bins in centrality
 151   const Int_t nbincent_0_10 = 4;  //bins in centrality between 0 and 10
 152   const Int_t nbincent_10_60 = 20;  //bins in centrality between 10 and 60
 153   const Int_t nbincent_60_100 = 4;  //bins in centrality between 60 and 100
 154   const Int_t nbinfake = 3;  //bins in fake
 155   const Int_t nbinpointingXY = 50;  //bins in cosPointingAngleXY
 156   const Int_t nbinnormDecayLXY = 20;  //bins in NormDecayLengthXY
 157   const Int_t nbinmult = 49;  //bins in multiplicity (total number)
 158   const Int_t nbinmult_0_20 = 20; //bins in multiplicity between 0 and 20
 159   const Int_t nbinmult_20_50 = 15; //bins in multiplicity between 20 and 50
 160   const Int_t nbinmult_50_80 = 10; //bins in multiplicity between 50 and 80
 161   const Int_t nbinmult_80_100 = 4; //bins in multiplicity between 80 and 100
 162   const Int_t nbinmult_100_400 = 3; //bins in multiplicity between 100 and 400
 163   if(isPPbData) nbinmult += nbinmult_100_400;
 164
 165
 166   //the sensitive variables, their indices
 167
 168   const Int_t nvarTot   = 16 ; //number of variables on the grid:pt, y, cosThetaStar, pTpi, pTk, cT, dca, d0pi, d0K, d0xd0, cosPointingAngle, phi, z, centrality, fake, cosPointingAngleXY, normDecayLengthXY, multiplicity
 169
 170   // variables' indices
 171   const UInt_t ipT = 0;
 172   const UInt_t iy  = 1;
 173   const UInt_t icosThetaStar  = 2;
 174   const UInt_t ipTpi  = 3;
 175   const UInt_t ipTk  = 4;
 176   const UInt_t icT  = 5;
 177   const UInt_t idca  = 6;
 178   const UInt_t id0xd0  = 7;
 179   const UInt_t ipointing  = 8;
 180   const UInt_t iphi  = 9;
 181   const UInt_t izvtx  = 10;
 182   const UInt_t icent = 11;
 183   const UInt_t ifake = 12;
 184   const UInt_t ipointingXY = 13;
 185   const UInt_t inormDecayLXY = 14;
 186   const UInt_t imult = 15;
 187
 188   //Setting the bins: pt, ptPi, and ptK are considered seprately because for them you can either define the binning by hand, or using the cuts file
 189
 190   //arrays for the number of bins in each dimension
 191   Int_t iBin[nvarTot];
 192
 193   const Int_t nbinpt = cutsLcTMVA->GetNPtBins(); // bins in pT
 194   iBin[ipT] = nbinpt;
 195   iBin[ipTpi] = nbinpt;
 196   iBin[ipTk] = nbinpt;
 197   Double_t *binLimpT = new Double_t[iBin[ipT]+1];
 198   Double_t *binLimpTpi = new Double_t[iBin[ipTpi]+1];
 199   Double_t *binLimpTk = new Double_t[iBin[ipTk]+1];
 200   // values for bin lower bounds
 201   Float_t* floatbinLimpT = cutsLcTMVA->GetPtBinLimits();
 202   for (Int_t ibin0 = 0 ; ibin0<iBin[ipT]+1; ibin0++){
 203     binLimpT[ibin0] = (Double_t)floatbinLimpT[ibin0];
 204     binLimpTpi[ibin0] = (Double_t)floatbinLimpT[ibin0];
 205     binLimpTk[ibin0] = (Double_t)floatbinLimpT[ibin0];
 206   }
 207   for(Int_t i=0; i<=nbinpt; i++) printf("binLimpT[%d]=%f\n",i,binLimpT[i]);
 208
 209   printf("pT: nbin (from cuts file) = %d\n",nbinpt);
 210
 211   // Fine Ntrk binning setting
 212   Double_t *binLimmultFine;
 213   Int_t nbinmultTmp = nbinmult;
 214   if (isFineNtrkBin) {
 215     Int_t nbinLimmultFine = 100;
 216     if (isPPbData) nbinLimmultFine = 200;
 217     const UInt_t nbinMultFine = nbinLimmultFine;
 218     binLimmultFine = new Double_t[nbinMultFine+1];
 219     for (Int_t ibin0 = 0; ibin0 < nbinMultFine+1; ibin0++) {
 220       binLimmultFine[ibin0] = ibin0;
 221     }
 222     nbinmultTmp = nbinLimmultFine;
 223   }
 224   const Int_t nbinmultTot = nbinmultTmp;
 225
 226   // defining now the binning for the other variables:
 227
 228   AliLog::SetClassDebugLevel("AliCFManager",AliLog::kInfo);
 229
 230   iBin[iy] = nbiny;
 231   iBin[icosThetaStar] = nbincosThetaStar;
 232   iBin[icT] = nbincT;
 233   iBin[idca] = nbindca;
 234   iBin[id0xd0] = nbind0xd0;
 235   iBin[ipointing] = nbinpointing;
 236   iBin[iphi] = nbinphi;
 237   iBin[izvtx] = nbinzvtx;
 238   iBin[icent] = nbincent;
 239   iBin[ifake] = nbinfake;
 240   iBin[ipointingXY] = nbinpointingXY;
 241   iBin[inormDecayLXY] = nbinnormDecayLXY;
 242   iBin[imult] = nbinmultTot;
 243
 244   //arrays for lower bounds :
 245   Double_t *binLimy = new Double_t[iBin[iy]+1];
 246   Double_t *binLimcosThetaStar = new Double_t[iBin[icosThetaStar]+1];
 247   Double_t *binLimcT = new Double_t[iBin[icT]+1];
 248   Double_t *binLimdca = new Double_t[iBin[idca]+1];
 249   Double_t *binLimd0xd0 = new Double_t[iBin[id0xd0]+1];
 250   Double_t *binLimpointing = new Double_t[iBin[ipointing]+1];
 251   Double_t *binLimphi = new Double_t[iBin[iphi]+1];
 252   Double_t *binLimzvtx = new Double_t[iBin[izvtx]+1];
 253   Double_t *binLimcent = new Double_t[iBin[icent]+1];
 254   Double_t *binLimfake = new Double_t[iBin[ifake]+1];
 255   Double_t *binLimpointingXY = new Double_t[iBin[ipointingXY]+1];
 256   Double_t *binLimnormDecayLXY = new Double_t[iBin[inormDecayLXY]+1];
 257   Double_t *binLimmult = new Double_t[iBin[imult]+1];
 258
 259
 260   // y
 261   for(Int_t i=0; i<=nbiny; i++) binLimy[i] = (Double_t)ymin  + (ymax-ymin)  /nbiny*(Double_t)i ;
 262
 263   // cosThetaStar
 264   for(Int_t i=0; i<=nbincosThetaStar; i++) binLimcosThetaStar[i] = (Double_t)cosminTS  + (cosmaxTS-cosminTS)  /nbincosThetaStar*(Double_t)i ;
 265
 266   // cT
 267   for(Int_t i=0; i<=nbincT; i++) binLimcT[i] = (Double_t)cTmin  + (cTmax-cTmin)  /nbincT*(Double_t)i ;
 268
 269   // dca
 270   for(Int_t i=0; i<=nbindca; i++) binLimdca[i] = (Double_t)dcamin  + (dcamax-dcamin)  /nbindca*(Double_t)i ;
 271
 272   // d0xd0
 273   for(Int_t i=0; i<=nbind0xd0; i++) binLimd0xd0[i] = (Double_t)d0xd0min  + (d0xd0max-d0xd0min)  /nbind0xd0*(Double_t)i ;
 274
 275   // cosPointingAngle
 276   for(Int_t i=0; i<=nbinpointing; i++) binLimpointing[i] = (Double_t)cosmin  + (cosmax-cosmin)  /nbinpointing*(Double_t)i ;
 277
 278   // Phi
 279   for(Int_t i=0; i<=nbinphi; i++) binLimphi[i] = (Double_t)phimin  + (phimax-phimin)  /nbinphi*(Double_t)i ;
 280
 281   // z Primary Vertex
 282   for(Int_t i=0; i<=nbinzvtx; i++) {
 283     binLimzvtx[i] = (Double_t)zmin  + (zmax-zmin)  /nbinzvtx*(Double_t)i ;
 284   }
 285
 286   // centrality
 287   for(Int_t i=0; i<=nbincent_0_10; i++) binLimcent[i] = (Double_t)centmin_0_10 + (centmax_0_10-centmin_0_10)/nbincent_0_10*(Double_t)i ;
 288   if (binLimcent[nbincent_0_10] != centmin_10_60)  {
 289     Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for cent - 1st range - differs from expected!\n");
 290   }
 291   for(Int_t i=0; i<=nbincent_10_60; i++) binLimcent[i+nbincent_0_10] = (Double_t)centmin_10_60 + (centmax_10_60-centmin_10_60)/nbincent_10_60*(Double_t)i ;
 292   if (binLimcent[nbincent_0_10+nbincent_10_60] != centmin_60_100)  {
 293     Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for cent - 2st range - differs from expected!\n");
 294   }
 295   for(Int_t i=0; i<=nbincent_60_100; i++) binLimcent[i+nbincent_0_10+nbincent_10_60] = (Double_t)centmin_60_100 + (centmax_60_100-centmin_60_100)/nbincent_60_100*(Double_t)i ;
 296
 297   // fake
 298   for(Int_t i=0; i<=nbinfake; i++) {
 299     binLimfake[i] = (Double_t)fakemin  + (fakemax-fakemin)/nbinfake * (Double_t)i;
 300   }
 301
 302   // cosPointingAngleXY
 303   for(Int_t i=0; i<=nbinpointingXY; i++) binLimpointingXY[i] = (Double_t)cosminXY  + (cosmaxXY-cosminXY)  /nbinpointingXY*(Double_t)i ;
 304
 305   // normDecayLXY
 306   for(Int_t i=0; i<=nbinnormDecayLXY; i++) binLimnormDecayLXY[i] = (Double_t)normDecLXYmin  + (normDecLXYmax-normDecLXYmin)  /nbinnormDecayLXY*(Double_t)i ;
 307
 308
 309   // multiplicity
 310   for(Int_t i=0; i<=nbinmult_0_20; i++) binLimmult[i] = (Double_t)multmin_0_20 + (multmax_0_20-multmin_0_20)/nbinmult_0_20*(Double_t)i ;
 311   if (binLimmult[nbinmult_0_20] != multmin_20_50)  {
 312     Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for mult - 1st range - differs from expected!\n");
 313   }
 314   for(Int_t i=0; i<=nbinmult_20_50; i++) binLimmult[i+nbinmult_0_20] = (Double_t)multmin_20_50 + (multmax_20_50-multmin_20_50)/nbinmult_20_50*(Double_t)i ;
 315   if (binLimmult[nbinmult_0_20+nbinmult_20_50] != multmin_50_80)  {
 316     Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for mult - 2nd range - differs from expected!\n");
 317   }
 318   for(Int_t i=0; i<=nbinmult_50_80; i++) binLimmult[i+nbinmult_0_20+nbinmult_20_50] = (Double_t)multmin_50_80 + (multmax_50_80-multmin_50_80)/nbinmult_50_80*(Double_t)i ;
 319   if (binLimmult[nbinmult_0_20+nbinmult_20_50+nbinmult_50_80] != multmin_80_100)  {
 320     Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for mult - 3rd range - differs from expected!\n");
 321   }
 322   for(Int_t i=0; i<=nbinmult_80_100; i++) binLimmult[i+nbinmult_0_20+nbinmult_20_50+nbinmult_50_80] = (Double_t)multmin_80_100 + (multmax_80_100-multmin_80_100)/nbinmult_80_100*(Double_t)i ;
 323   if (binLimmult[nbinmult_0_20+nbinmult_20_50+nbinmult_50_80+nbinmult_80_100] != multmin_100_400) {
 324     Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for mult - 4th range - differs from expected!\n");
 325   }
 326   if(isPPbData){
 327     for (Int_t i = 0; i<=nbinmult_100_400; i++) binLimmult[i+nbinmult_0_20+nbinmult_20_50+nbinmult_50_80+nbinmult_80_100] = (Double_t)multmin_100_400 + (multmax_100_400-multmin_100_400)/nbinmult_100_400*(Double_t)i ;
 328   }
 329
 330   //one "container" for MC
 331   TString nameContainer="";
 332   if(!isKeepDfromB) {
 333     nameContainer = "CFHFccontainer0";
 334   }
 335   else  if(isKeepDfromBOnly){
 336     nameContainer = "CFHFccontainer0DfromB";
 337   }
 338   else  {
 339     nameContainer = "CFHFccontainer0allD";
 340   }
 341   nameContainer += suffix;
 342   //Setting up the container grid...
 343
 344   AliCFContainer* container;
 345
 346   if (configuration == AliCFTaskVertexingHF::kSnail){
 347     container = new AliCFContainer(nameContainer,"container for tracks",nstep,nvarTot,iBin);
 348     //setting the bin limits
 349     printf("pt\n");
 350     container -> SetBinLimits(ipT, binLimpT);
 351     printf("y\n");
 352     container -> SetBinLimits(iy, binLimy);
 353     printf("cts\n");
 354     container -> SetBinLimits(icosThetaStar, binLimcosThetaStar);
 355     printf("ptPi\n");
 356     container -> SetBinLimits(ipTpi, binLimpTpi);
 357     printf("ptK\n");
 358     container -> SetBinLimits(ipTk, binLimpTk);
 359     printf("cT\n");
 360     container -> SetBinLimits(icT, binLimcT);
 361     printf("dca\n");
 362     container -> SetBinLimits(idca, binLimdca);
 363     printf("d0xd0\n");
 364     container -> SetBinLimits(id0xd0, binLimd0xd0);
 365     printf("pointing\n");
 366     container -> SetBinLimits(ipointing, binLimpointing);
 367     printf("phi\n");
 368     container -> SetBinLimits(iphi, binLimphi);
 369     printf("z\n");
 370     container -> SetBinLimits(izvtx, binLimzvtx);
 371     printf("cent\n");
 372     container -> SetBinLimits(icent, binLimcent);
 373     printf("fake\n");
 374     container -> SetBinLimits(ifake, binLimfake);
 375     printf("pointingXY\n");
 376     container -> SetBinLimits(ipointingXY, binLimpointingXY);
 377     printf("normDecayLXY\n");
 378     container -> SetBinLimits(inormDecayLXY, binLimnormDecayLXY);
 379     printf("multiplicity\n");
 380
 381     if (isFineNtrkBin) container->SetBinLimits(imult, binLimmultFine);
 382     else               container->SetBinLimits(imult, binLimmult);
 383
 384     container -> SetVarTitle(ipT,"pt");
 385     container -> SetVarTitle(iy,"y");
 386     container -> SetVarTitle(icosThetaStar, "cosThetaStar");
 387     container -> SetVarTitle(ipTpi, "ptpi");
 388     container -> SetVarTitle(ipTk, "ptK");
 389     container -> SetVarTitle(icT, "ct");
 390     container -> SetVarTitle(idca, "dca");
 391     container -> SetVarTitle(id0xd0, "d0xd0");
 392     container -> SetVarTitle(ipointing, "pointing");
 393     container -> SetVarTitle(iphi, "phi");
 394     container -> SetVarTitle(izvtx, "zvtx");
 395     container -> SetVarTitle(icent, "centrality");
 396     container -> SetVarTitle(ifake, "fake");
 397     container -> SetVarTitle(ipointingXY, "piointingXY");
 398     container -> SetVarTitle(inormDecayLXY, "normDecayLXY");
 399     container -> SetVarTitle(imult, "multiplicity");
 400   }
 401   else if (configuration == AliCFTaskVertexingHF::kCheetah){
 402     //arrays for the number of bins in each dimension
 403     const Int_t nvar = 8;
 404
 405     const UInt_t ipTFast = 0;
 406     const UInt_t iyFast = 1;
 407     const UInt_t icTFast = 2;
 408     const UInt_t iphiFast = 3;
 409     const UInt_t izvtxFast = 4;
 410     const UInt_t icentFast = 5;
 411     const UInt_t ifakeFast = 6;
 412     const UInt_t imultFast = 7;
 413
 414     Int_t iBinFast[nvar];
 415     iBinFast[ipTFast] = iBin[ipT];
 416     iBinFast[iyFast] = iBin[iy];
 417     iBinFast[icTFast] = iBin[icT];
 418     iBinFast[iphiFast] = iBin[iphi];
 419     iBinFast[izvtxFast] = iBin[izvtx];
 420     iBinFast[icentFast] = iBin[icent];
 421     iBinFast[ifakeFast] = iBin[ifake];
 422     iBinFast[imultFast] = iBin[imult];
 423
 424     container = new AliCFContainer(nameContainer,"container for tracks", nstep, nvar, iBinFast);
 425     printf("pt\n");
 426     container -> SetBinLimits(ipTFast, binLimpT);
 427     printf("y\n");
 428     container -> SetBinLimits(iyFast, binLimy);
 429     printf("ct\n");
 430     container -> SetBinLimits(icTFast, binLimcT);
 431     printf("phi\n");
 432     container -> SetBinLimits(iphiFast, binLimphi);
 433     printf("zvtx\n");
 434     container -> SetBinLimits(izvtxFast, binLimzvtx);
 435     printf("centrality\n");
 436     container -> SetBinLimits(icentFast, binLimcent);
 437     printf("fake\n");
 438     container -> SetBinLimits(ifakeFast, binLimfake);
 439     printf("multiplicity\n");
 440     if(isFineNtrkBin) container -> SetBinLimits(imultFast, binLimmultFine);
 441     else              container -> SetBinLimits(imultFast, binLimmult);
 442
 443     container -> SetVarTitle(ipTFast, "pt");
 444     container -> SetVarTitle(iyFast, "y");
 445     container -> SetVarTitle(icTFast, "ct");
 446     container -> SetVarTitle(iphiFast, "phi");
 447     container -> SetVarTitle(izvtxFast, "zvtx");
 448     container -> SetVarTitle(icentFast, "centrality");
 449     container -> SetVarTitle(ifakeFast, "fake");
 450     container -> SetVarTitle(imultFast, "multiplicity");
 451   }
 452
 453   container -> SetStepTitle(0, "MCLimAcc");
 454   container -> SetStepTitle(1, "MC");
 455   container -> SetStepTitle(2, "MCAcc");
 456   container -> SetStepTitle(3, "RecoVertex");
 457   container -> SetStepTitle(4, "RecoRefit");
 458   container -> SetStepTitle(5, "Reco");
 459   container -> SetStepTitle(6, "RecoAcc");
 460   container -> SetStepTitle(7, "RecoITSCluster");
 461   container -> SetStepTitle(8, "RecoCuts");
 462   container -> SetStepTitle(9, "RecoPID");
 463
 464   //return container;
 465
 466   //CREATE THE  CUTS -----------------------------------------------
 467
 468   // Gen-Level kinematic cuts
 469   AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
 470
 471   //Particle-Level cuts:
 472   AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
 473   Bool_t useAbsolute = kTRUE;
 474   if (isSign != 2){
 475     useAbsolute = kFALSE;
 476   }
 477   mcGenCuts->SetRequirePdgCode(pdgCode, useAbsolute);  // kTRUE set in order to include D0_bar
 478   mcGenCuts->SetAODMC(1); //special flag for reading MC in AOD tree (important)
 479
 480   // Acceptance cuts:
 481   AliCFAcceptanceCuts* accCuts = new AliCFAcceptanceCuts("accCuts", "Acceptance cuts");
 482   AliCFTrackKineCuts *kineAccCuts = new AliCFTrackKineCuts("kineAccCuts","Kine-Acceptance cuts");
 483   kineAccCuts->SetPtRange(ptmin, ptmax);
 484   kineAccCuts->SetEtaRange(etamin, etamax);
 485
 486   // Rec-Level kinematic cuts
 487   AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts", "rec-level kine cuts");
 488
 489   AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts", "rec-level quality cuts");
 490
 491   AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts", "rec-level isPrimary cuts");
 492
 493   printf("CREATE MC KINE CUTS\n");
 494   TObjArray* mcList = new TObjArray(0) ;
 495   mcList->AddLast(mcKineCuts);
 496   mcList->AddLast(mcGenCuts);
 497
 498   printf("CREATE ACCEPTANCE CUTS\n");
 499   TObjArray* accList = new TObjArray(0) ;
 500   accList->AddLast(kineAccCuts);
 501
 502   printf("CREATE RECONSTRUCTION CUTS\n");
 503   TObjArray* recList = new TObjArray(0) ;   // not used!!
 504   recList->AddLast(recKineCuts);
 505   recList->AddLast(recQualityCuts);
 506   recList->AddLast(recIsPrimaryCuts);
 507
 508   TObjArray* emptyList = new TObjArray(0);
 509
 510   //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
 511   printf("CREATE INTERFACE AND CUTS\n");
 512   AliCFManager* man = new AliCFManager() ;
 513   man->SetParticleContainer(container);
 514   man->SetParticleCutsList(0 , mcList); // MC, Limited Acceptance
 515   man->SetParticleCutsList(1 , mcList); // MC
 516   man->SetParticleCutsList(2 , accList); // Acceptance
 517   man->SetParticleCutsList(3 , emptyList); // Vertex
 518   man->SetParticleCutsList(4 , emptyList); // Refit
 519   man->SetParticleCutsList(5 , emptyList); // AOD
 520   man->SetParticleCutsList(6 , emptyList); // AOD in Acceptance
 521   man->SetParticleCutsList(7 , emptyList); // AOD with required n. of ITS clusters
 522   man->SetParticleCutsList(8 , emptyList); // AOD Reco (PPR cuts implemented in Task)
 523   man->SetParticleCutsList(9 , emptyList); // AOD Reco PID
 524
 525   // Get the pointer to the existing analysis manager via the static access method.
 526   //==============================================================================
 527   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 528   if (!mgr) {
 529     ::Error("AddTaskCompareHF", "No analysis manager to connect to.");
 530     return NULL;
 531   }
 532   //CREATE THE TASK
 533   printf("CREATE TASK\n");
 534
 535   // create the task
 536   AliCFTaskVertexingHF *task = new AliCFTaskVertexingHF("AliCFTaskVertexingHF", cutsLcTMVA);
 537   task->SetConfiguration(configuration);
 538   task->SetFillFromGenerated(kFALSE);
 539   task->SetCFManager(man); //here is set the CF manager
 540   task->SetDecayChannel(22);
 541   task->SetUseCascadeTaskForLctoV0bachelor(kTRUE);
 542   task->SetUseAdditionalCuts(kTRUE);
 543   task->SetUseCutsForTMVA(kFALSE);
 544   task->SetUseFlatPtWeight(useFlatPtWeight);
 545   task->SetUseWeight(useWeight);
 546   task->SetUseZWeight(useZWeight);
 547   task->SetSign(isSign);
 548   task->SetCentralitySelection(kFALSE);
 549   task->SetFakeSelection(0);
 550   //  task->SetRejectCandidateIfNotFromQuark(kTRUE); // put to false if you want to keep HIJING D0!!
 551   task->SetRejectCandidateIfNotFromQuark(kFALSE); // put to false if you want to keep HIJING D0!!
 552   task->SetUseMCVertex(kFALSE); // put to true if you want to do studies on pp
 553   //task->SetPtWeightsFromDataPbPb276overLHC12a17a();
 554   task->SetCutOnMomConservation(0.0000005);
 555
 556   if (isKeepDfromB && !isKeepDfromBOnly) task->SetDselection(2);
 557   if (isKeepDfromB && isKeepDfromBOnly) task->SetDselection(1);
 558
 559   TF1* funcWeight = 0x0;
 560   if (task->GetUseWeight()) {
 561     funcWeight = (TF1*)fileCuts->Get("funcWeight");
 562     if (funcWeight == 0x0){
 563       Printf("FONLL Weights will be used");
 564     }
 565     else {
 566       task->SetWeightFunction(funcWeight);
 567       Printf("User-defined Weights will be used. The function being:");
 568       task->GetWeightFunction()->Print();
 569     }
 570   }
 571
 572   if (useNchWeight || useNtrkWeight){
 573     TH1F *hNchPrimaries;
 574     TH1F *hNchMeasured;
 575     if (isPPbData) hNchPrimaries = (TH1F*)fileCuts->Get("hNtrUnCorrEvWithCandWeight");
 576     else hNchPrimaries = (TH1F*)fileCuts->Get("hGenPrimaryParticlesInelGt0");
 577     hNchMeasured = (TH1F*)fileCuts->Get("hNchMeasured");
 578     if(hNchPrimaries) {
 579       task->SetUseNchWeight(kTRUE);
 580       task->SetMCNchHisto(hNchPrimaries);
 581       if(isPPbData) task->SetUseNchTrackletsWeight();
 582     } else {
 583       AliFatal("Histogram for multiplicity weights not found");
 584       return 0x0;
 585     }
 586     if(hNchMeasured) task->SetMeasuredNchHisto(hNchMeasured);
 587     if(useNtrkWeight) task->SetUseNchTrackletsWeight();
 588   }
 589
 590   task->SetMultiplicityEstimator(multiplicityEstimator);
 591   if(estimatorFilename.EqualTo("") ) {
 592     printf("Estimator file not provided, multiplicity corrected histograms will not be filled\n");
 593     task->SetUseZvtxCorrectedNtrkEstimator(kFALSE);
 594   } else{
 595     TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
 596     if(!fileEstimator)  {
 597       AliFatal("File with multiplicity estimator not found");
 598       return;
 599     }
 600     task->SetUseZvtxCorrectedNtrkEstimator(kTRUE);
 601     task->SetReferenceMultiplcity(refMult);
 602
 603     if (isPPbData) {  //load multiplicity estimators for pPb
 604       task->SetIsPPbData(kTRUE);
 605       const Char_t* periodNames[2] = {"LHC13b", "LHC13c"};
 606       TProfile *multEstimatorAvg[2];
 607       for (Int_t ip=0; ip < 2; ip++) {
 608         multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));
 609         if (!multEstimatorAvg[ip]) {
 610           AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
 611           return;
 612         }
 613       }
 614       task->SetMultiplVsZProfileLHC13b(multEstimatorAvg[0]);
 615       task->SetMultiplVsZProfileLHC13c(multEstimatorAvg[1]);
 616
 617     } else {    //load multiplicity estimators for pp
 618       const Char_t* periodNames[4] = {"LHC10b", "LHC10c", "LHC10d", "LHC10e"};
 619       TProfile* multEstimatorAvg[4];
 620
 621       for(Int_t ip=0; ip<4; ip++) {
 622         multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));
 623         if (!multEstimatorAvg[ip]) {
 624           AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
 625           return;
 626         }
 627       }
 628       task->SetMultiplVsZProfileLHC10b(multEstimatorAvg[0]);
 629       task->SetMultiplVsZProfileLHC10c(multEstimatorAvg[1]);
 630       task->SetMultiplVsZProfileLHC10d(multEstimatorAvg[2]);
 631       task->SetMultiplVsZProfileLHC10e(multEstimatorAvg[3]);
 632     }
 633   }
 634   Printf("***************** CONTAINER SETTINGS *****************");
 635   Printf("decay channel = %d", (Int_t)task->GetDecayChannel());
 636   Printf("FillFromGenerated = %d", (Int_t)task->GetFillFromGenerated());
 637   Printf("Dselection = %d", (Int_t)task->GetDselection());
 638   Printf("UseWeight = %d", (Int_t)task->GetUseWeight());
 639   if (task->GetUseWeight()) {
 640     if(funcWeight) Printf("User-defined Weight function");
 641     else Printf("FONLL will be used for the weights");
 642   }
 643
 644   Printf("Use Nch weight = %d", (Int_t)task->GetUseNchWeight());
 645   Printf("Sign = %d", (Int_t)task->GetSign());
 646   Printf("Centrality selection = %d", (Int_t)task->GetCentralitySelection());
 647   Printf("Fake selection = %d", (Int_t)task->GetFakeSelection());
 648   Printf("RejectCandidateIfNotFromQuark selection = %d", (Int_t)task->GetRejectCandidateIfNotFromQuark());
 649   Printf("UseMCVertex selection = %d", (Int_t)task->GetUseMCVertex());
 650   Printf("UseCascadeTask = %d", (Int_t)task->GetUseCascadeTaskForLctoV0bachelor());
 651   Printf("***************END CONTAINER SETTINGS *****************\n");
 652
 653
 654   //-----------------------------------------------------------//
 655   //   create correlation matrix for unfolding - only eta-pt   //
 656   //-----------------------------------------------------------//
 657
 658   Bool_t AcceptanceUnf = kTRUE; // unfold at acceptance level, otherwise PPR
 659
 660   Int_t thnDim[4];
 661
 662   //first half  : reconstructed
 663   //second half : MC
 664
 665   thnDim[0] = iBin[0];
 666   thnDim[2] = iBin[0];
 667   thnDim[1] = iBin[1];
 668   thnDim[3] = iBin[1];
 669
 670   TString nameCorr="";
 671   if(!isKeepDfromB) {
 672     nameCorr="CFHFcorr0";
 673   }
 674   else  if(isKeepDfromBOnly){
 675     nameCorr= "CFHFcorr0KeepDfromBOnly";
 676   }
 677   else  {
 678     nameCorr="CFHFcorr0allD";
 679
 680   }
 681   nameCorr += suffix;
 682
 683   THnSparseD* correlation = new THnSparseD(nameCorr,"THnSparse with correlations",4,thnDim);
 684   Double_t** binEdges = new Double_t[2];
 685
 686   // set bin limits
 687
 688   binEdges[0]= binLimpT;
 689   binEdges[1]= binLimy;
 690
 691   correlation->SetBinEdges(0,binEdges[0]);
 692   correlation->SetBinEdges(2,binEdges[0]);
 693
 694   correlation->SetBinEdges(1,binEdges[1]);
 695   correlation->SetBinEdges(3,binEdges[1]);
 696
 697   correlation->Sumw2();
 698
 699   // correlation matrix ready
 700   //------------------------------------------------//
 701
 702   task->SetCorrelationMatrix(correlation); // correlation matrix for unfolding
 703
 704   // Create and connect containers for input/output
 705
 706   // ------ input data ------
 707   AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();
 708
 709   // ----- output data -----
 710
 711   TString outputfile = AliAnalysisManager::GetCommonFileName();
 712   TString output1name="", output2name="", output3name="",output4name="", output5name="";
 713   output2name=nameContainer;
 714   output3name=nameCorr;
 715   output4name= "Cuts";
 716   output5name= "coutProfDst";
 717   if(!isKeepDfromB) {
 718     outputfile += ":PWG3_D2H_CFtaskLctoK0SpCascade";
 719     output1name="CFHFchist0";
 720     output3name+="_cOnly";
 721     output4name+="_cOnly";
 722     output5name+="_cOnly";
 723   }
 724   else  if(isKeepDfromBOnly){
 725     outputfile += ":PWG3_D2H_CFtaskLctoK0SpCascadeKeepDfromBOnly";
 726     output1name="CFHFchist0DfromB";
 727     output3name+="_bOnly";
 728     output4name+="_bOnly";
 729     output5name+="_bOnly";
 730   }
 731   else{
 732     outputfile += ":PWG3_D2H_CFtaskLctoK0SpCascadeKeepDfromB";
 733     output1name="CFHFchist0allD";
 734     output3name+="_all";
 735     output4name+="_all";
 736     output5name+="_all";
 737   }
 738
 739   outputfile += suffix;
 740   output1name += suffix;
 741   output4name += suffix;
 742   output5name += suffix;
 743
 744   //now comes user's output objects :
 745   // output TH1I for event counting
 746   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(output1name, TH1I::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
 747   // output Correction Framework Container (for acceptance & efficiency calculations)
 748   AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(output2name, AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
 749   // Unfolding - correlation matrix
 750   AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(output3name, THnSparseD::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
 751   // cuts
 752   AliAnalysisDataContainer *coutput4 = mgr->CreateContainer(output4name, AliRDHFCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());
 753   // estimators list
 754   AliAnalysisDataContainer *coutput5 = mgr->CreateContainer(output5name, TList::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());
 755
 756   mgr->AddTask(task);
 757
 758   mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
 759   mgr->ConnectOutput(task,1,coutput1);
 760   mgr->ConnectOutput(task,2,coutput2);
 761   mgr->ConnectOutput(task,3,coutput3);
 762   mgr->ConnectOutput(task,4,coutput4);
 763   mgr->ConnectOutput(task,5,coutput5);
 764   return task;
 765
 766 }