PWGHF/hfe/macros/configs/PbPb/ConfigHFE_FLOW_TOFTPC.C

   1 TH2D *GetCorrectionsJpsiMean(TString map){
   2   if (gSystem->AccessPathName(gSystem->ExpandPathName(map.Data()))){
   3     Warning("ConfigHFE","Centrality map not found: %s",map.Data());
   4     printf("ConfigHFE: Centrality map not found: %s\n",map.Data());
   5     return kFALSE;
   6   }
   7   printf("File %s\n",map.Data());
   8   TFile *f = TFile::Open(map.Data());
   9   if(!f->IsOpen()) return kFALSE;
  10   f->ls();
  11   return (TH2D*)f->Get("centroidCentEta");
  12 }
  13 TH2D *GetCorrectionsJpsiWidth(TString map){
  14   if (gSystem->AccessPathName(gSystem->ExpandPathName(map.Data()))){
  15     Warning("ConfigHFE","Centrality map not found: %s",map.Data());
  16     printf("ConfigHFE: Centrality map not found: %s\n",map.Data());
  17     return kFALSE;
  18   }
  19   printf("File %s\n",map.Data());
  20   TFile *f = TFile::Open(map.Data());
  21   if(!f->IsOpen()) return kFALSE;
  22   f->ls();
  23   return (TH2D*)f->Get("widthCentEta");
  24 }
  25 TF1* GetCentralityCorrection(TString listname="LHC11h"){
  26
  27   TString etaMap="$ALICE_ROOT/PWGHF/hfe/macros/configs/PbPb/CentCorrMapsTPC.root";
  28
  29   if (gSystem->AccessPathName(gSystem->ExpandPathName(etaMap.Data()))){
  30     Error("ConfigHFEpbpb","Eta map not found: %s",etaMap.Data());
  31     return 0;
  32   }
  33
  34   TFile f(etaMap.Data());
  35   if (!f.IsOpen()) return 0;
  36   gROOT->cd();
  37   TList *keys=f.GetListOfKeys();
  38
  39   for (Int_t i=0; i<keys->GetEntries(); ++i){
  40     TString kName=keys->At(i)->GetName();
  41     TPRegexp reg(kName);
  42     if (reg.MatchB(listname)){
  43       printf("Using Eta Correction Function: %s\n",kName.Data());
  44       return (TF1*)f.Get(kName.Data());
  45     }
  46   }
  47   return 0;
  48 }
  49 TF1* GetEtaCorrection(TString listname="LHC11h"){
  50
  51   TString etaMap="$ALICE_ROOT/PWGHF/hfe/macros/configs/PbPb/EtaCorrMapsTPC.root";
  52
  53   if (gSystem->AccessPathName(gSystem->ExpandPathName(etaMap.Data()))){
  54     Error("ConfigHFEpbpb","Eta map not found: %s",etaMap.Data());
  55     return 0;
  56   }
  57
  58   TFile f(etaMap.Data());
  59   if (!f.IsOpen()) return 0;
  60   gROOT->cd();
  61   TList *keys=f.GetListOfKeys();
  62
  63   for (Int_t i=0; i<keys->GetEntries(); ++i){
  64     TString kName=keys->At(i)->GetName();
  65     TPRegexp reg(kName);
  66     if (reg.MatchB(listname)){
  67       printf("Using Eta Correction Function: %s\n",kName.Data());
  68       return (TF1*)f.Get(kName.Data());
  69     }
  70   }
  71   return 0;
  72 }
  73 TF1* Getv2Contamination_30_40(){
  74
  75   TString v2Map="$ALICE_ROOT/PWGHF/hfe/macros/configs/PbPb/raw_TF1_pi_v2_3040.root";
  76   //TString v2Map="$TRAIN_ROOT/bailhach_backgroundhfe/raw_TF1_pi_v2_3040.root";
  77   TString trainRoot=gSystem->Getenv("TRAIN_ROOT");
  78
  79   if (gSystem->AccessPathName(gSystem->ExpandPathName(v2Map.Data()))){
  80     //Error("ConfigPbPb2010_Cent","v2 map not found: %s",v2Map.Data());
  81     return 0;
  82   }
  83
  84   TFile f(v2Map.Data());
  85   if (!f.IsOpen()) return 0;
  86   gROOT->cd();
  87   TF1 *fg = (TF1 *) f.Get("flogisticTadd");
  88   return fg;
  89
  90 }
  91
  92 TF1* Getv2Contamination_40_50(){
  93
  94   TString v2Map="$ALICE_ROOT/PWGHF/hfe/macros/configs/PbPb/raw_TF1_pi_v2_4050.root";
  95   //TString v2Map="$TRAIN_ROOT/bailhach_backgroundhfe/raw_TF1_pi_v2_4050.root";
  96   TString trainRoot=gSystem->Getenv("TRAIN_ROOT");
  97
  98   if (gSystem->AccessPathName(gSystem->ExpandPathName(v2Map.Data()))){
  99     //Error("ConfigPbPb2010_Cent","v2 map not found: %s",v2Map.Data());
 100     return 0;
 101   }
 102
 103   TFile f(v2Map.Data());
 104   if (!f.IsOpen()) return 0;
 105   gROOT->cd();
 106   TF1 *fg = (TF1 *) f.Get("flogisticTadd");
 107   return fg;
 108
 109 }
 110
 111
 112 TF1* Getv2Contamination_20_30(){
 113
 114   TString v2Map="$ALICE_ROOT/PWGHF/hfe/macros/configs/PbPb/raw_TF1_pi_v2_2030.root";
 115   //TString v2Map="$TRAIN_ROOT/bailhach_backgroundhfe/raw_TF1_pi_v2_2030.root";
 116   TString trainRoot=gSystem->Getenv("TRAIN_ROOT");
 117
 118   if (gSystem->AccessPathName(gSystem->ExpandPathName(v2Map.Data()))){
 119     //Error("ConfigPbPb2010_Cent","v2 map not found: %s",v2Map.Data());
 120     return 0;
 121   }
 122
 123   TFile f(v2Map.Data());
 124   if (!f.IsOpen()) return 0;
 125   gROOT->cd();
 126   TF1 *fg = (TF1 *) f.Get("flogisticTadd");
 127   return fg;
 128
 129 }
 130
 131 Double_t Contamination_40_50(const Double_t *x, const Double_t *par)
 132 {
 133
 134   if(x[0] < 2.5) return 0.0326072;
 135   Double_t value = -0.00908646-0.00201815*x[0]+0.0026871*x[0]*x[0];
 136   if(x[0] >= 2.5) return value;
 137
 138
 139 }
 140 AliAnalysisTaskFlowTPCTOFEPSP* ConfigHFE_FLOW_TOFTPC(Bool_t useMC, TString appendix,Int_t trigger,Int_t aodfilter=-1,Bool_t scalarProduct=kFALSE,Bool_t cutPileup=kTRUE,Int_t variableMr = 1,Int_t tpcCls=110, Double_t tpcClsr=60.,Int_t tpcClspid=80, Int_t itsCls=4, Int_t pixellayer=2, Double_t dcaxy=100., Double_t dcaz=200.,  Double_t tofsig=30., Double_t *tpcdedx=NULL, Int_t vzero=1, Int_t debuglevel=0, Double_t etarange=80, Bool_t withetacorrection=kFALSE, Bool_t withmultcorrection=kFALSE, Double_t ITSclustersback=0,Double_t minTPCback=-2.0,Double_t maxTPCback=5.0)
 141 {
 142   //
 143   // HFE flow task
 144   //
 145   Double_t tpcsharedfraction=11;
 146   Double_t chi2peritscl=36.;
 147
 148   // multEst = 0 or 10 (VZERO centrality, good, bad runs for Jpsi), = 1 (GetNumberOfESDtracks, Theo)
 149   Int_t variableM = variableMr;
 150   if(variableMr==10) variableM = 0;
 151
 152   printf("Summary settings flow task\n");
 153   printf("filter %d\n",aodfilter);
 154   printf("TPC number of tracking clusters %d\n",tpcCls);
 155   printf("TPC ratio clusters %f\n",tpcClsr*0.01);
 156   printf("TPC number of pid clusters %d\n",tpcClspid);
 157   printf("Maximal fraction of TPC shared cluster %f\n",tpcsharedfraction*0.01);
 158   printf("ITS number of clusters %d\n",itsCls);
 159   printf("Maximal chi2 per ITS cluster %f\n",chi2peritscl);
 160   printf("Requirement on the pixel layer %d\n",pixellayer);
 161   printf("dcaxy %f\n",dcaxy*0.01);
 162   printf("dcaz %f\n",dcaz*0.01);
 163   printf("TOF sigma %f\n",tofsig*0.1);
 164   printf("TPC min sigma cut 0: %f\n",tpcdedx[0]);
 165   printf("TPC min sigma cut 1: %f\n",tpcdedx[1]);
 166   printf("TPC min sigma cut 2: %f\n",tpcdedx[2]);
 167   printf("TPC min sigma cut 3: %f\n",tpcdedx[3]);
 168   printf("TPC min sigma cut 4: %f\n",tpcdedx[4]);
 169   printf("TPC min sigma cut 5: %f\n",tpcdedx[5]);
 170   printf("TPC min sigma cut 6: %f\n",tpcdedx[6]);
 171   printf("TPC min sigma cut 7: %f\n",tpcdedx[7]);
 172   printf("VZERO event plane %d\n",vzero);
 173   printf("Debug level %d\n",debuglevel);
 174   printf("Etarange %f\n",etarange*0.01);
 175   printf("TPC dE/dx Eta correction %d\n",withetacorrection);
 176   printf("TPC dE/dx multiplicity correction %d\n",withmultcorrection);
 177   printf("Number of ITS back clusters %d\n",(Int_t)ITSclustersback);
 178   printf("Min TPC back %f\n",minTPCback);
 179   printf("Max TPC back %f\n",maxTPCback);
 180   printf("PileUp cut %d\n",cutPileup);
 181   printf("Scalar Product %d\n",scalarProduct);
 182   printf("Multiplicity Esimator %d\n",variableM);
 183
 184   // Cut HFE
 185   AliHFEcuts *hfecuts = new AliHFEcuts("hfeCuts","HFE Standard Cuts");
 186   hfecuts->CreateStandardCuts();
 187   hfecuts->SetEtaRange(etarange*0.01);
 188   hfecuts->SetMinNClustersTPC(tpcCls);
 189   hfecuts->SetMinNClustersTPCPID(tpcClspid);
 190   hfecuts->SetMinRatioTPCclusters(tpcClsr*0.01);
 191   hfecuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
 192   hfecuts->SetFractionOfSharedTPCClusters(tpcsharedfraction*0.01);
 193   hfecuts->SetMinNClustersITS(itsCls);
 194   hfecuts->SetCutITSpixel(pixellayer);
 195   hfecuts->SetCheckITSLayerStatus(kFALSE);
 196   hfecuts->SetMaxImpactParam(dcaxy*0.01,dcaz*0.01);
 197   hfecuts->SetPtRange(0.1,10.);
 198   hfecuts->SetAcceptKinkMothers();
 199
 200   hfecuts->SetVertexRange(10.);
 201   hfecuts->SetUseSPDVertex(kTRUE);
 202   hfecuts->SetUseCorrelationVertex();
 203
 204   hfecuts->SetTOFPIDStep(kTRUE);
 205
 206   // Cut HFE background
 207   AliESDtrackCuts *hfeBackgroundCuts = new AliESDtrackCuts();
 208   hfeBackgroundCuts->SetName("backgroundcuts");
 209   //hfeBackgroundCuts->SetAcceptKinkDaughters(kFALSE);
 210   hfeBackgroundCuts->SetRequireTPCRefit(kTRUE);
 211   hfeBackgroundCuts->SetRequireITSRefit(kTRUE);
 212   if(ITSclustersback>0) hfeBackgroundCuts->SetMinNClustersITS(ITSclustersback);
 213   hfeBackgroundCuts->SetEtaRange(-0.8,0.8);
 214   hfeBackgroundCuts->SetRequireSigmaToVertex(kTRUE);
 215   hfeBackgroundCuts->SetMaxChi2PerClusterTPC(3.5);
 216   hfeBackgroundCuts->SetMinNClustersTPC(100);
 217   hfeBackgroundCuts->SetPtRange(0.5,1e10);
 218
 219   // Name
 220   printf("appendix %s\n", appendix.Data());
 221
 222   // The task
 223   AliAnalysisTaskFlowTPCTOFEPSP *task = new AliAnalysisTaskFlowTPCTOFEPSP(Form("HFE_%s", appendix.Data()));
 224   //task->SelectCollisionCandidates(trigger);
 225   task->SetTriggerUsed(trigger);
 226   task->SetDebugLevel(1);
 227   task->SetVariableMultiplicity(variableM);
 228   //task->GetPIDQAManager()->SetHighResolutionEtaHistos();
 229   task->GetPIDQAManager()->SetMidResolutionHistos();
 230   //task->GetPIDQAManager()->SetFillMultiplicity();
 231   task->SetHFECuts(hfecuts);
 232   task->SetHFEBackgroundCuts(hfeBackgroundCuts);
 233   if(aodfilter > 0) {
 234     printf("ON AOD filter %d\n",aodfilter);
 235     task->SetFilter(aodfilter);
 236   }
 237   else {
 238     printf("Default AOD filter 1<<4\n");
 239   }
 240   if(useMC) {
 241     task->SetMCPID(kTRUE);
 242     //task->SetUseMCReactionPlane(kTRUE);
 243     task->SetAfterBurnerOn(kTRUE);
 244     task->SetV1V2V3V4V5(0.0,0.2,0.0,0.0,0.0);
 245   }
 246   if(vzero>=1) task->SetVZEROEventPlane(kTRUE);
 247   if(vzero==2) task->SetVZEROEventPlaneA(kTRUE);
 248   if(vzero==3) task->SetVZEROEventPlaneC(kTRUE);
 249   if((vzero==4) && useMC) task->SetUseMCReactionPlane(kTRUE);
 250
 251   // Debug level
 252   task->SetDebugLevel(debuglevel);
 253   if(debuglevel==1) task->SetMonitorEventPlane(kTRUE);
 254   if(debuglevel==2) task->SetMonitorContamination(kTRUE);
 255   if(debuglevel==3) task->SetMonitorPhotonic(kTRUE);
 256   if(debuglevel==4) task->SetMonitorWithoutPID(kTRUE);
 257   if(debuglevel==5) task->SetMonitorTrackCuts(kTRUE);
 258   if(debuglevel==6) task->SetMonitorQCumulant(kTRUE);
 259
 260
 261   // Define PID
 262   AliHFEpid *pid = task->GetPID();
 263   if(useMC) pid->SetHasMCData(kTRUE);
 264   pid->AddDetector("TOF", 0);
 265   pid->AddDetector("TPC", 1);
 266
 267   if(withetacorrection || withmultcorrection) {
 268     AliHFEpidTPC *tpcpid = pid->GetDetPID(AliHFEpid::kTPCpid);
 269     // Jpsi
 270     if(variableMr==0) {
 271       TH2D *meanc = GetCorrectionsJpsiMean("$ALICE_ROOT/PWGHF/hfe/macros/configs/PbPb/jpsietacentcorr2_2011.root");
 272       TH2D *widthc = GetCorrectionsJpsiWidth("$ALICE_ROOT/PWGHF/hfe/macros/configs/PbPb/jpsietacentcorr2_2011.root");
 273       //printf("Set\n");
 274       if(meanc && widthc) {
 275         tpcpid->SetJpsiCorrections(meanc,widthc);
 276         //printf("Set the histos\n");
 277       }
 278     }
 279     if(variableMr==10) {
 280       TH2D *meanc = GetCorrectionsJpsiMean("$ALICE_ROOT/PWGHF/hfe/macros/configs/PbPb/jpsietacentcorr3_2011.root");
 281       TH2D *widthc = GetCorrectionsJpsiWidth("$ALICE_ROOT/PWGHF/hfe/macros/configs/PbPb/jpsietacentcorr3_2011.root");
 282       if(meanc && widthc) tpcpid->SetJpsiCorrections(meanc,widthc);
 283     }
 284     // Theo
 285     if(variableMr==1 || variableMr==2) {
 286       task->GetPIDQAManager()->SetFillMultiplicity();
 287       TF1 *etaCorrMean = GetEtaCorrection("LHC11h_etaCorrMean");
 288       TF1 *etaCorrWdth = GetEtaCorrection("LHC11h_etaCorrWidth");
 289       if(etaCorrMean && etaCorrWdth && withetacorrection){
 290         tpcpid->SetEtaCorrections(etaCorrMean, etaCorrWdth);
 291         printf("TPC dE/dx Eta correction %p %p\n",etaCorrMean,etaCorrWdth);
 292       }
 293       TF1 *centCorrMean = GetCentralityCorrection("LHC11h_multCorrMean");
 294       TF1 *centCorrWdth = GetCentralityCorrection("LHC11h_multCorrWidth");
 295       if(centCorrMean && centCorrWdth && withmultcorrection){
 296         tpcpid->SetCentralityCorrections(centCorrMean, centCorrWdth);
 297         printf("TPC dE/dx multiplicity correction %p %p\n",centCorrMean,centCorrWdth);
 298       }
 299     }
 300   }
 301
 302   task->SetPileUpCut(cutPileup);
 303   task->SetUseSP(scalarProduct);
 304
 305   //pid->AddDetector("BAYES", 0);
 306   //pid->ConfigureBayesDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF+AliPIDResponse::kDetTRD);
 307   //pid->ConfigureBayesPIDThreshold(0.9);
 308
 309
 310   TString datatype=gSystem->Getenv("CONFIG_FILE");
 311
 312   if(!useMC) {
 313
 314     if(withetacorrection || withmultcorrection) {
 315
 316       pid->ConfigureTPCasymmetric(0.0,9999.,tpcdedx[0],3.);
 317
 318     } else {
 319
 320       Double_t params_centr_0_5[1];
 321       Double_t params_centr_5_10[1];
 322       Double_t params_centr_10_20[1];
 323       Double_t params_centr_20_30[1];
 324       Double_t params_centr_30_40[1];
 325       Double_t params_centr_40_50[1];
 326       Double_t params_centr_50_60[1];
 327       Double_t params_centr_per[1];
 328
 329       params_centr_0_5[0]=tpcdedx[0];  // cut tuned for 0-5%
 330       params_centr_5_10[0]=tpcdedx[1]; // cut tuned for 5-10%
 331       params_centr_10_20[0]=tpcdedx[2];
 332       params_centr_20_30[0]=tpcdedx[3];
 333       params_centr_30_40[0]=tpcdedx[4];
 334       params_centr_40_50[0]=tpcdedx[5];
 335       params_centr_50_60[0]=tpcdedx[6];
 336       params_centr_per[0]=tpcdedx[7];
 337
 338       char *cutmodel;
 339       cutmodel="pol0";
 340
 341       for(Int_t a=0;a<11;a++)
 342         {
 343           if(a>6)  pid->ConfigureTPCcentralityCut(a,cutmodel,params_centr_per,3.0);      //  60-80%
 344           if(a==0) pid->ConfigureTPCcentralityCut(a,cutmodel,params_centr_0_5,3.0);      //  0-5%
 345           if(a==1) pid->ConfigureTPCcentralityCut(a,cutmodel,params_centr_5_10,3.0);     //  5-10%
 346           if(a==2) pid->ConfigureTPCcentralityCut(a,cutmodel,params_centr_10_20,3.0);    //  10-20%
 347           if(a==3) pid->ConfigureTPCcentralityCut(a,cutmodel,params_centr_20_30,3.0);    //  20-30%
 348           if(a==4) pid->ConfigureTPCcentralityCut(a,cutmodel,params_centr_30_40,3.0);    //  30-40%
 349           if(a==5) pid->ConfigureTPCcentralityCut(a,cutmodel,params_centr_40_50,3.0);    //  40-50%
 350           if(a==6) pid->ConfigureTPCcentralityCut(a,cutmodel,params_centr_50_60,3.0);    //  50-60%
 351         }
 352     }
 353
 354   }
 355
 356   if(tofsig>0.) pid->ConfigureTOF(tofsig*0.1);
 357
 358   AliHFEpidTOF *tofpid = pid->GetDetPID(AliHFEpid::kTOFpid);
 359   if(tofsig<=0.) {
 360     Double_t uppercuttof = TMath::Abs(tofsig)*0.1;
 361     Double_t paramsTOFlowercut[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};
 362     Double_t paramsTOFuppercut[12] ={uppercuttof,uppercuttof,uppercuttof,uppercuttof,uppercuttof,uppercuttof,uppercuttof,uppercuttof,uppercuttof,uppercuttof,uppercuttof,uppercuttof};
 363     for(Int_t a=0;a<11;a++)
 364       {
 365         tofpid->SetTOFnSigmaBandCentrality(paramsTOFlowercut[a],paramsTOFuppercut[a],a);
 366       }
 367   }
 368
 369
 370   // Define hadron contamination (Only for 2011, TPC middle)
 371   TF1 *hBackground_20_30 = new TF1("hadronicBackgroundFunction_20_30","[0]+[1]*x+[2]*x*x+[3]*x*x*x", 0., 200.);
 372   hBackground_20_30->SetParameter(0, -0.165789);
 373   hBackground_20_30->SetParameter(1, 0.218694);
 374   hBackground_20_30->SetParameter(2, -0.076635);
 375   hBackground_20_30->SetParameter(3, 0.00947502);
 376   task->SetContamination(hBackground_20_30,3);
 377   TF1 *fv2_20_30 = Getv2Contamination_20_30();
 378   if(fv2_20_30) {
 379     fv2_20_30->SetName("fv2_20_30");
 380     task->SetV2Contamination(fv2_20_30,3);
 381   }
 382   // 30-40%
 383   TF1 *hBackground_30_40 = new TF1("hadronicBackgroundFunction_30_40","[0]+[1]*x+[2]*x*x+[3]*x*x*x", 0., 200.);
 384   hBackground_30_40->SetParameter(0, -0.072222);
 385   hBackground_30_40->SetParameter(1, 0.132098);
 386   hBackground_30_40->SetParameter(2, -0.0561759);
 387   hBackground_30_40->SetParameter(3, 0.00789356);
 388   task->SetContamination(hBackground_30_40,4);
 389   TF1 *fv2_30_40 = Getv2Contamination_30_40();
 390   if(fv2_30_40) {
 391     fv2_30_40->SetName("fv2_30_40");
 392     task->SetV2Contamination(fv2_30_40,4);
 393   }
 394   // 40-50%
 395   TF1 *hBackground_40_50 = new TF1("hadronicBackgroundFunction_40_50",Contamination_40_50,0.,200.,0);
 396   printf("contamination 5 in pt range 0.5 %f\n",hBackground_40_50->Eval(0.5));
 397   printf("contamination 5 in pt range 1.5 %f\n",hBackground_40_50->Eval(1.5));
 398   printf("contamination 5 in pt range 3.5 %f\n",hBackground_40_50->Eval(3.5));
 399   task->SetContamination(hBackground_40_50,5);
 400   TF1 *fv2_40_50 = Getv2Contamination_40_50();
 401   if(fv2_40_50) {
 402     fv2_40_50->SetName("fv2_40_50");
 403     task->SetV2Contamination(fv2_40_50,5);
 404   }
 405
 406   // Define PID TOF Only
 407   AliHFEpid *pidTOFOnly = task->GetPIDTOFOnly();
 408   if(useMC) pidTOFOnly->SetHasMCData(kTRUE);
 409   pidTOFOnly->AddDetector("TOF", 0);
 410   pidTOFOnly->ConfigureTOF(tofsig*0.1);
 411
 412   // Define PID background
 413   AliHFEpid *pidbackground = task->GetPIDBackground();
 414   if(useMC) pidbackground->SetHasMCData(kTRUE);
 415   //pidbackground->AddDetector("TOF", 0);
 416   pidbackground->AddDetector("TPC", 0);
 417   //pidbackground->ConfigureTOF(3.0);
 418   pidbackground->ConfigureTPCasymmetric(0.0,9999.,minTPCback,maxTPCback);
 419   task->SetMaxopeningtheta(9999.0);
 420   task->SetMaxopeningphi(9999.0);
 421   task->SetMaxopening3D(0.1);
 422   // Always AliKF and no mass constraint
 423   task->SetAlgorithmMA(kFALSE);
 424   task->SetMassConstraint(kFALSE);
 425
 426   printf("*************************************\n");
 427   printf("Configuring standard Task:\n");
 428   task->Print();
 429   pid->PrintStatus();
 430   printf("*************************************\n");
 431   return task;
 432
 433 }