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

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