PWGHF/hfe/macros/configs/pPb/ConfigHFEnpepPb.C

   1 Bool_t ReadContaminationFunctionsBeauty(TString filename, TF1 **functions, double sigma){
   2   TFile *in = TFile::Open(Form("$ALICE_ROOT/PWGHF/hfe/macros/configs/pPb/%s", filename.Data()));
   3   gROOT->cd();
   4   int isig = static_cast<int>(sigma * 100.);
   5   if (isig == -44) isig = -42;
   6   if (isig == 6) isig = 9;
   7   printf("Getting hadron background for the sigma cut: %d\n", isig);
   8   bool status = kTRUE;
   9   for(int icent = 0; icent < 12; icent++){
  10     functions[icent] = dynamic_cast<TF1 *>(in->Get(Form("hback_%d_%d", isig, icent)));
  11     if(functions[icent]) printf("Config for centrality class %d found\n", icent);
  12     else{
  13       printf("Config for the centrality class %d not found\n", icent);
  14       status = kFALSE;
  15     }
  16   }
  17   delete in;
  18   return status;
  19 }
  20
  21 Bool_t ReadContaminationFunctions(TString filename, TF1 **functions, double sigma){
  22   TFile *in = TFile::Open(Form("$ALICE_ROOT/PWGHF/hfe/macros/configs/pPb/%s", filename.Data()));
  23   gROOT->cd();
  24   int isig = static_cast<int>(sigma * 100.);
  25   printf("Getting hadron background for the sigma cut: %d\n", isig);
  26   bool status = kTRUE;
  27   for(int icent = 0; icent < 12; icent++){
  28     functions[icent] = dynamic_cast<TF1 *>(in->Get(Form("hback_%d_%d", isig, icent)));
  29     if(functions[icent]) printf("Config for centrality class %d found\n", icent);
  30     else{
  31       printf("Config for the centrality class %d not found\n", icent);
  32       status = kFALSE;
  33     }
  34   }
  35   delete in;
  36   return status;
  37 }
  38
  39
  40 AliAnalysisTaskHFE* ConfigHFEnpepPb(Bool_t useMC, Bool_t isAOD, Bool_t isBeauty, TString appendix,
  41                 UChar_t TPCcl=70, UChar_t TPCclPID = 80,
  42                 UChar_t ITScl=3, Double_t DCAxy=1000., Double_t DCAz=1000.,
  43                 Double_t* tpcdEdxcutlow=NULL, Double_t* tpcdEdxcuthigh=NULL,
  44                 Double_t TOFs=3., Int_t TOFmis=0,
  45                 Int_t itshitpixel = 0, Int_t icent,
  46                 Double_t etami=-0.8, Double_t etama=0.8,
  47                 Double_t phimi=-1., Double_t phima=-1.,
  48                 Double_t assETAm=-0.8, Double_t assETAp=0.8,
  49                 Double_t assMinPt=0.2, Int_t assITS=2,
  50                 Int_t assTPCcl=100, Int_t assTPCPIDcl=80,
  51                 Double_t assDCAr=1.0, Double_t assDCAz=2.0,
  52                 Double_t *assTPCSminus=NULL, Double_t *assTPCSplus=NULL,
  53                 Double_t assITSpid=-3.,
  54                 Double_t assTOFs=3.,
  55                 Bool_t useCat1Tracks = kTRUE, Bool_t useCat2Tracks = kTRUE, Int_t weightlevelback = -1,
  56                 Bool_t nonPhotonicElectronBeauty = kFALSE, Bool_t ipCharge = kFALSE, Bool_t ipOpp = kFALSE, Int_t ipsys = 0)
  57 {
  58   Bool_t kAnalyseTaggedTracks = kTRUE;
  59   Bool_t kApplyPreselection = kFALSE;
  60
  61   //***************************************//
  62   //        Setting up the HFE cuts        //
  63   //***************************************//
  64
  65   AliHFEcuts *hfecuts = new AliHFEcuts(appendix,"HFE cuts for pPb");
  66   //hfecuts->SetQAOn();
  67   hfecuts->CreateStandardCuts();
  68   hfecuts->SetMinNClustersTPC(TPCcl);
  69   hfecuts->SetMinNClustersTPCPID(TPCclPID);
  70   hfecuts->SetMinNClustersITS(ITScl);
  71   hfecuts->SetMinRatioTPCclusters(0.6);
  72   hfecuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
  73   hfecuts->SetCutITSpixel(itshitpixel);
  74   hfecuts->SetCheckITSLayerStatus(kFALSE);
  75   hfecuts->SetEtaRange(etami,etama);
  76   if(phimi >= 0. && phima >= 0) hfecuts->SetPhiRange(phimi,phima);
  77   hfecuts->SetRejectKinkDaughters();
  78   hfecuts->SetAcceptKinkMothers();
  79   if(isAOD) hfecuts->SetAODFilterBit(4);
  80
  81   //if((iPixelAny==AliHFEextraCuts::kAny) || (iPixelAny==AliHFEextraCuts::kSecond))
  82
  83   hfecuts->SetMaxImpactParam(DCAxy,DCAz);
  84   hfecuts->SetUseMixedVertex(kTRUE);
  85   hfecuts->SetVertexRange(10.);
  86   // New pPb cuts (February 2013)
  87   hfecuts->SetUseCorrelationVertex();
  88   hfecuts->SetSPDVtxResolutionCut();
  89   hfecuts->SetpApileupCut();
  90
  91   Bool_t ipSig = kFALSE;
  92
  93   hfecuts->SetIPcutParam(0.0054,0.078,-0.56,0,ipSig,ipCharge,ipOpp);
  94   if(ipsys==1) hfecuts->SetIPcutParam(0.0054,0.057,-0.66,0,ipSig,ipCharge,ipOpp);
  95   if(ipsys==2) hfecuts->SetIPcutParam(0.012,0.088,-0.65,0,ipSig,ipCharge,ipOpp);
  96
  97   if(isBeauty) hfecuts->SetProductionVertex(0,100,0,100);
  98
  99   // TOF settings:
 100   Int_t usetof=0;
 101   Bool_t kTOFmis=kFALSE;
 102   if (TOFs>0.){
 103     usetof = 1;
 104     printf("CONFIGURATION FILE: TOF is used \n");
 105     hfecuts->SetTOFPIDStep(kTRUE);
 106     printf("CONFIGURATION FILE: TOF PID step is requested !!!! \n");
 107     if (TOFmis>0){
 108       kTOFmis = kTRUE;
 109       printf("CONFIGURATION FILE: TOF mismatch rejection is set ON \n");
 110     }
 111   }
 112
 113   //***************************************//
 114   //        Setting up the task            //
 115   //***************************************//
 116
 117   AliAnalysisTaskHFE *task = new AliAnalysisTaskHFE(Form("HFEtask%s",appendix.Data()));
 118   printf("task %p\n", task);
 119   task->SetpPbAnalysis();
 120   if(!isAOD) task->SetRemoveFirstEventInChunk();
 121   task->SetRemovePileUp(kFALSE);
 122   task->SetHFECuts(hfecuts);
 123   task->GetPIDQAManager()->SetHighResolutionHistos();
 124   task->SetRejectKinkMother(kFALSE);
 125
 126   // Determine the centrality estimator
 127   task->SetCentralityEstimator("V0A");
 128   if (icent == 2) task->SetCentralityEstimator("V0M");
 129   else if (icent == 3) task->SetCentralityEstimator("CL1");
 130   else if (icent == 4) task->SetCentralityEstimator("ZNA");
 131
 132   //***************************************//
 133   //        Prepare preselection           //
 134   // This mimics the ESD->AOD filter in    //
 135   // case of the ESD analysis and selects  //
 136   // only tracks which will be selected in //
 137   // the AOD analysis with the given filter//
 138   // bit. Not to be applied for AODS.      //
 139   // For pPb the cuts used are (bit 4)     //
 140   // esdTrackCutsHG0 from file $ALICE_ROOT///
 141   // ANALYSIS/macros/AddTaskESDFilter.C    //
 142   //***************************************//
 143
 144   if(kApplyPreselection){
 145     AliESDtrackCuts* esdfilter = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 146     esdfilter->SetMaxDCAToVertexXY(2.4);
 147     esdfilter->SetMaxDCAToVertexZ(3.2);
 148     esdfilter->SetDCAToVertex2D(kTRUE);
 149
 150     task->SetHFECutsPreselect(esdfilter);
 151     printf("Put a preselection cut\n");
 152     task->SetFillNoCuts(kTRUE);
 153   }
 154
 155   //***************************************//
 156   //          Variable manager             //
 157   //***************************************//
 158   // Define Variables
 159   Double_t ptbinning[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
 160   Double_t etabinning[17] = {-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};
 161
 162   Int_t sizept=(sizeof(ptbinning)/sizeof(double))-1;
 163   Int_t sizeeta=(sizeof(etabinning)/sizeof(double))-1;
 164
 165   AliHFEvarManager *vm = task->GetVarManager();
 166   vm->AddVariable("pt", sizept, ptbinning);
 167   vm->AddVariable("eta", sizeeta, -0.8,0.8);
 168   vm->AddVariable("phi",18, -0, 2*TMath::Pi());
 169   vm->AddVariable("charge");
 170   vm->AddVariable("source");
 171   vm->AddVariable("centrality");
 172
 173   // For the moment, remove the part dedicated to the background subtraction.
 174   // It will be implemented in a different way, reading it from a root file.
 175
 176   //***************************************//
 177   //          Configure the PID            //
 178   //***************************************//
 179
 180   // Define PID
 181   AliHFEpid *pid = task->GetPID();
 182   if(useMC) pid->SetHasMCData(kTRUE);
 183
 184   if (usetof){
 185     pid->AddDetector("TOF", 0);
 186     pid->AddDetector("TPC", 1);
 187   } else {
 188     pid->AddDetector("TPC", 0);
 189   }
 190
 191   // Configure TPC PID
 192   // do the identical thing in data and MC
 193   Double_t paramsTPCdEdxcutlow[12] ={0.0, 0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
 194   if(tpcdEdxcutlow) memcpy(paramsTPCdEdxcutlow,tpcdEdxcutlow,sizeof(paramsTPCdEdxcutlow));
 195
 196   Double_t paramsTPCdEdxcuthigh[12] ={3.0, 3.0, 3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
 197   if(tpcdEdxcuthigh) memcpy(paramsTPCdEdxcuthigh,tpcdEdxcuthigh,sizeof(paramsTPCdEdxcuthigh));
 198
 199   char *cutmodel;
 200
 201   if(useMC){ // constant (default) cut for MC
 202       cutmodel="pol0(0)";
 203       Double_t params[1];
 204       params[0]=paramsTPCdEdxcutlow[0];
 205       pid->ConfigureTPCdefaultCut(cutmodel, params,tpcdEdxcuthigh[0]);
 206   } else { // correct for mean shift in data
 207       cutmodel="min(pol1(0),pol0(2))";
 208       Double_t params[3];
 209       //params[0]=-0.12; params[1]=0.14; params[2]=0.09;
 210       params[0]=-0.21 + paramsTPCdEdxcutlow[0];
 211       params[1]=0.14;
 212       params[2]=paramsTPCdEdxcutlow[0];
 213       pid->ConfigureTPCdefaultCut(cutmodel, params,tpcdEdxcuthigh[0]);
 214   }
 215   /*
 216   char *cutmodel;
 217   cutmodel="pol0";
 218
 219   for(Int_t a=0;a<11;a++){
 220     // Not necessary anymore, since the pPb case is handled similarly to the pp case
 221     //   cout << a << " " << paramsTPCdEdxcut[a] << endl;
 222     Double_t tpcparamlow[1]={paramsTPCdEdxcutlow[a]};
 223     Float_t tpcparamhigh=paramsTPCdEdxcuthigh[a];
 224     pid->ConfigureTPCcentralityCut(a,cutmodel,tpcparamlow,tpcparamhigh);
 225   }
 226   pid->ConfigureTPCdefaultCut(cutmodel,paramsTPCdEdxcutlow,paramsTPCdEdxcuthigh[0]); // After introducing the pPb flag, pPb is merged with pp and this line defines the cut
 227   */
 228
 229   // Configure TOF PID
 230   if (usetof){
 231     pid->ConfigureTOF(TOFs);
 232     AliHFEpidTOF *tofpid = pid->GetDetPID(AliHFEpid::kTOFpid);
 233     if (kTOFmis){
 234       tofpid->SetRejectTOFmismatch();
 235     }
 236   }
 237
 238   // To make different upper TOF cut to see contamination effect
 239   // The below two lines should be removed after this check
 240   //AliHFEpidTOF *tofpid = pid->GetDetPID(AliHFEpid::kTOFpid);
 241   //if(TOFs<3.) tofpid->SetTOFnSigmaBand(-3,TOFs); //only to check the assymmetric tof cut
 242
 243   // Load hadron background
 244   if(!useMC){
 245     Bool_t status = kTRUE;
 246     TF1 *hBackground[12];
 247     if(isAOD==1) {
 248         if (usetof)  status = ReadContaminationFunctions("hadroncontamination_AOD139_TOFPID_pPb_eta06.root", hBackground, tpcdEdxcutlow[0]);
 249         else status = ReadContaminationFunctions("hadroncontamination_AOD139_noTOFPID_pPb_eta06.root", hBackground, tpcdEdxcutlow[0]);
 250     }
 251     else if (isBeauty==1) {
 252         status = ReadContaminationFunctionsBeauty("hadroncontamination_ESD_Beauty_TOFPID_pPb_eta06.root", hBackground, tpcdEdxcutlow[0]);
 253     }
 254   else  status = ReadContaminationFunctions("hadroncontamination_TOFTPC_pPb_eta06_newsplines_try3.root", hBackground, tpcdEdxcutlow[0]);
 255     for(Int_t a=0;a<12;a++) {
 256       //printf("back %f \n",hBackground[a]);
 257       if(status) task->SetBackGroundFactorsFunction(hBackground[a],a);
 258       else printf("not all background functions found\n");
 259     }
 260   }
 261
 262   //***************************************//
 263   //       Configure NPE plugin            //
 264   //***************************************//
 265
 266   AliHFENonPhotonicElectron *backe = new AliHFENonPhotonicElectron(Form("HFEBackGroundSubtractionPID2%s",appendix.Data()),"Background subtraction");  //appendix
 267     //Setting the Cuts for the Associated electron-pool
 268   AliHFEcuts *hfeBackgroundCuts = new AliHFEcuts(Form("HFEBackSub%s",appendix.Data()),"Background sub Cuts");
 269   //  hfeBackgroundCuts->SetEtaRange(assETA);
 270   hfeBackgroundCuts->SetEtaRange(assETAm,assETAp);
 271   hfeBackgroundCuts->SetPtRange(assMinPt,20.);
 272
 273   hfeBackgroundCuts->SetMaxChi2perClusterTPC(4);
 274   hfeBackgroundCuts->SetMinNClustersITS(assITS);
 275   hfeBackgroundCuts->SetMinNClustersTPC(assTPCcl);
 276   hfeBackgroundCuts->SetMinNClustersTPCPID(assTPCPIDcl);
 277   hfeBackgroundCuts->SetMaxImpactParam(assDCAr,assDCAz);
 278   if(isAOD) hfeBackgroundCuts->SetAODFilterBit(4);
 279   hfeBackgroundCuts->SetQAOn();                         // QA
 280
 281   AliHFEpid *pidbackground = backe->GetPIDBackground();
 282   if(useMC) pidbackground->SetHasMCData(kTRUE);
 283
 284   if (assTOFs>0.){
 285     pidbackground->AddDetector("TOF", 0);
 286     pidbackground->AddDetector("TPC", 1);
 287   } else {
 288     pidbackground->AddDetector("TPC", 0);
 289   }
 290
 291   Double_t paramsTPCdEdxcutlowAssoc[12] ={-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0,-3.0};
 292   if(assTPCSminus) memcpy(paramsTPCdEdxcutlowAssoc,assTPCSminus,sizeof(paramsTPCdEdxcutlowAssoc));
 293
 294   Double_t paramsTPCdEdxcuthighAssoc[12] ={3.0, 3.0, 3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
 295   if(assTPCSplus) memcpy(paramsTPCdEdxcuthighAssoc,assTPCSplus,sizeof(paramsTPCdEdxcuthighAssoc));
 296
 297   char *cutmodelAssoc;
 298   cutmodelAssoc="pol0";
 299   for(Int_t a=0;a<11;a++){
 300     // Not necessary anymore, since the pPb case is handled similarly to the pp case
 301     //   cout << a << " " << paramsTPCdEdxcut[a] << endl;
 302     Double_t tpcparamlow[1]={paramsTPCdEdxcutlowAssoc[a]};
 303     Float_t tpcparamhigh=paramsTPCdEdxcuthighAssoc[a];
 304     pidbackground->ConfigureTPCcentralityCut(a,cutmodelAssoc,tpcparamlow,tpcparamhigh);
 305   }
 306   pidbackground->ConfigureTPCdefaultCut(cutmodelAssoc,paramsTPCdEdxcutlowAssoc,paramsTPCdEdxcuthighAssoc[0]); // After introducing the pPb flag, pPb is merged with pp and this line defines the cut
 307   //backe->GetPIDBackgroundQAManager()->SetHighResolutionHistos();
 308
 309   if (assTOFs>0.){
 310     pidbackground->ConfigureTOF(TOFs);
 311   }
 312
 313   backe->SetHFEBackgroundCuts(hfeBackgroundCuts);
 314
 315   // Selection of associated tracks for the pool
 316   if(useCat1Tracks) backe->SelectCategory1Tracks(kTRUE);
 317   if(useCat2Tracks){
 318     backe->SelectCategory2Tracks(kTRUE);
 319     backe->SetITSMeanShift(-0.5);
 320     backe->SetITSnSigmaHigh(assITSpid);
 321     Double_t assITSminus = -1.0 * assITSpid;
 322     backe->SetITSnSigmaLow(assITSminus);
 323     //backe->SetminPt(assMinPt);
 324   }
 325
 326   // apply opening angle cut to reduce file size
 327   backe->SetMaxInvMass(0.3);
 328   backe->SetPtBinning(sizept, ptbinning);
 329   backe->SetEtaBinning(sizeeta, etabinning);
 330   // MC weight
 331   if(useMC) {
 332     //printf("test put weight %d\n",weightlevelback);
 333     if((weightlevelback >=0) && (weightlevelback < 3)) backe->SetWithWeights(weightlevelback);
 334   }
 335   task->SetHFEBackgroundSubtraction(backe);
 336
 337   //task->SetWeightHist();
 338
 339   //***************************************//
 340   //          V0 tagged tracks             //
 341   //***************************************//
 342
 343   if(kAnalyseTaggedTracks){
 344     AliHFEcuts *v0trackCuts = new AliHFEcuts("V0trackCuts", "Track Cuts for tagged track Analysis");
 345     v0trackCuts->CreateStandardCuts();
 346     v0trackCuts->SetMinNClustersTPC(TPCcl);
 347     v0trackCuts->SetMinNClustersTPCPID(TPCclPID);
 348     v0trackCuts->SetMinRatioTPCclusters(0.6);
 349     v0trackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
 350     v0trackCuts->SetMinNClustersITS(1);
 351     v0trackCuts->SetCutITSpixel(AliHFEextraCuts::kFirst);
 352     v0trackCuts->SetCheckITSLayerStatus(kFALSE);
 353     v0trackCuts->UnsetVertexRequirement();
 354     //hfecuts->SetSigmaToVertex(10);
 355     if(usetof) v0trackCuts->SetTOFPIDStep(kTRUE);
 356     v0trackCuts->SetQAOn();
 357
 358     task->SwitchOnPlugin(AliAnalysisTaskHFE::kTaggedTrackAnalysis);
 359     task->SetTaggedTrackCuts(v0trackCuts);
 360     task->SetCleanTaggedTrack(kTRUE);
 361   }
 362
 363   // QA
 364   printf("task %p\n", task);
 365   task->SetQAOn(AliAnalysisTaskHFE::kPIDqa);
 366   task->SetQAOn(AliAnalysisTaskHFE::kMCqa);
 367   task->SwitchOnPlugin(AliAnalysisTaskHFE::kDEstep);
 368   if(nonPhotonicElectronBeauty) task->SwitchOnPlugin(AliAnalysisTaskHFE::kNonPhotonicElectronBeauty);
 369   else task->SwitchOnPlugin(AliAnalysisTaskHFE::kNonPhotonicElectron);
 370
 371   printf("*************************************\n");
 372   printf("Configuring standard Task:\n");
 373   task->PrintStatus();
 374   pid->PrintStatus();
 375   printf("*************************************\n");
 376   return task;
 377 }