PWGLF/SPECTRA/PiKaPr/TOF/pPb502/macros/TOFpid.C

   1 #include "CommonDefs.C"
   2
   3 #define DISPLAY 0
   4
   5 /* SIGNAL SHAPE */
   6 Bool_t GAUSSIAN_SIGNAL = kFALSE;
   7 Bool_t GAUSSIANTAIL_SIGNAL = kFALSE;
   8 Bool_t GAUSSIANTAIL2_SIGNAL = kFALSE;
   9 Bool_t GAUSSIANPLUSGAUSSIANTAIL_SIGNAL = kFALSE;
  10 Bool_t GAUSSIANPLUSEXPONENTIAL_SIGNAL = kFALSE;
  11 Bool_t EXPECTED_SIGNAL_TAIL = kTRUE;
  12 Bool_t EXPECTED_SIGNAL_TEMPLATE = kTRUE;
  13 Bool_t EXPECTED_SIGNAL_FIT = kFALSE;
  14 /* SIGNAL PARAMETERS */
  15 Bool_t FIX_SIGNAL_MEAN = kFALSE;
  16 Bool_t FIX_SIGNAL_SIGMA = kFALSE;
  17 Bool_t FIX_SIGNAL_TAIL = kFALSE;
  18 Float_t SCALE_SIGNAL_SIGMA = 1.;
  19 Float_t SCALE_SIGNAL_TAIL = 1.;
  20 /* OTHER STUFF */
  21 Char_t *SIGNAL_PARAM_FILE = NULL;//"signalParamFile.root";
  22 Float_t DEFAULT_SIGNAL_MEAN = 0.;
  23 Float_t MIN_SIGNAL_MEAN = -0.2;
  24 Float_t MAX_SIGNAL_MEAN = 0.2.;
  25 Float_t DEFAULT_SIGNAL_SIGMA = 1.;
  26 Float_t MIN_SIGNAL_SIGMA = 0.8;
  27 Float_t MAX_SIGNAL_SIGMA = 1.2;
  28 Float_t DEFAULT_SIGNAL_TAIL = 1.;
  29 Float_t MIN_SIGNAL_TAIL = 0.5;
  30 Float_t MAX_SIGNAL_TAIL = 1.5;
  31 /* BACKGROUND */
  32 Bool_t EXPECTED_BACKGROUND_TAIL = kTRUE;
  33 Bool_t EXPECTED_BACKGROUND_TEMPLATE = kTRUE;
  34 Bool_t EXPECTED_BACKGROUND_FIT = kFALSE;
  35 Bool_t GAUSSIAN_BACKGROUND = kFALSE;
  36 Bool_t USE_ELECTRON_BACKGROUND = kTRUE;
  37 /* BACKGROUND PARAMETERS */
  38 Bool_t FIX_BACKGROUND_MEAN = kTRUE;
  39 Bool_t FIX_BACKGROUND_SIGMA = kTRUE;
  40 Bool_t FIX_BACKGROUND_TAIL = kTRUE;
  41 Float_t SCALE_BACKGROUND_SIGMA = 1.;
  42 Float_t SCALE_BACKGROUND_TAIL = 1.;
  43 /* MISMATCH */
  44 Bool_t NO_MISMATCH = kFALSE;
  45 Bool_t EXPECTED_MISMATCH = kTRUE;
  46 Bool_t EXPONENTIAL_MISMATCH = kFALSE;
  47 Bool_t UNIFORM_MISMATCH = kFALSE;
  48 Bool_t DOUBLEEXPONENTIAL_MISMATCH = kFALSE;
  49 Bool_t UNIFORMPLUSEXPONENTIAL_MISMATCH = kFALSE;
  50 /* AUTOMATIC CONSTRAINS */
  51 Float_t FIT_ELECTRONS_PT_MIN = 0.;
  52 Float_t FIT_ELECTRONS_PT_MAX = 1.0;
  53 Float_t FIT_MUONS_PT_MIN = 0.7;
  54 Float_t FIT_MUONS_PT_MAX = 0.;
  55 Float_t FIT_PIONS_PT_MIN = 0.;
  56 Float_t FIT_PIONS_PT_MAX = 5.;
  57
  58 Float_t CONSTRAINSIGNAL_LIMIT = 0.;
  59 Float_t CONSTRAINTAIL_LIMIT = 0.;
  60
  61 enum EFitParams_t {
  62   kMean,
  63   kSigma,
  64   kTail,
  65   kTotalCounts,
  66   kIntegralCounts,
  67   kSignalCounts,
  68   kBkg1Counts,
  69   kBkg2Counts,
  70   kBkg3Counts,
  71   kBkg4Counts,
  72   kMismatchCounts,
  73   kNFitParams
  74 };
  75
  76 /* fit output params name */
  77 const Char_t *fitParamName[kNFitParams] = {
  78   "Mean",
  79   "Sigma",
  80   "Tail",
  81   "TotalCounts",
  82   "IntegralCounts",
  83   "SignalCounts",
  84   "Bkg1Counts",
  85   "Bkg2Counts",
  86   "Bkg3Counts",
  87   "Bkg4Counts",
  88   "MismatchCounts"
  89 };
  90
  91 /* fit output params title */
  92 const Char_t *fitParamTitle[kNFitParams] = {
  93   "Signal mean;p_{T} (GeV/c);#mu (ps)",
  94   "Signal sigma;p_{T} (GeV/c);#sigma (ps)",
  95   "Signal tail;p_{T} (GeV/c);#sigma_{tail} (ps)",
  96   "Total counts;p_{T} (GeV/c);counts",
  97   "Total counts within fit range;p_{T} (GeV/c);counts",
  98   "Signal counts;p_{T} (GeV/c);counts",
  99   "Background-1 counts;p_{T} (GeV/c);counts",
 100   "Background-2 counts;p_{T} (GeV/c);counts",
 101   "Background-3 counts;p_{T} (GeV/c);counts",
 102   "Background-4 counts;p_{T} (GeV/c);counts",
 103   "Mismatch counts within fit range;p_{T} (GeV/c);counts"
 104 };
 105
 106 /**************************************************************/
 107 /*** GENERATION OF TEMPLATE HISTOS ****************************/
 108 /**************************************************************/
 109
 110 const Int_t NmismatchTrials = 1;
 111 const Int_t NexpectedTrials = 1;
 112
 113 /**************************************************************/
 114 /*** HISTOS AND BINNING ***************************************/
 115 /**************************************************************/
 116
 117 /**************************************************************/
 118 /**************************************************************/
 119 enum EHistoParam_t {
 120   kCentrality,
 121   kTOFsigma,
 122   kPt,
 123   kTPCsigma,
 124   kNHistoParams
 125 };
 126 /**************************************************************/
 127 /**************************************************************/
 128 const Int_t NtofsigmaBins = 1750;
 129 Double_t tofsigmaBin[NtofsigmaBins + 1];
 130 Double_t tofsigmaMin = -100., tofsigmaMax = 250., tofsigmaStep = (tofsigmaMax - tofsigmaMin) / NtofsigmaBins;
 131 /**************************************************************/
 132 const Int_t NtofsignalBins = 2000;
 133 Double_t tofsignalBin[NtofsignalBins + 1];
 134 Double_t tofsignalMin = -24400., tofsignalMax = 24400., tofsignalStep = (tofsignalMax - tofsignalMin) / NtofsignalBins;
 135 /**************************************************************/
 136 const Int_t NtofresoBins = 500;
 137 Double_t tofresoBin[NtofsignalBins + 1];
 138 Double_t tofresoMin = 0., tofresoMax = 500., tofresoStep = (tofresoMax - tofresoMin) / NtofresoBins;
 139 /**************************************************************/
 140 /**************************************************************/
 141 const Int_t NtpcsigmaBins = 10;
 142 Double_t tpcsigmaBin[NtpcsigmaBins + 1];
 143 Double_t tpcsigmaMin = -5., tpcsigmaMax = 5., tpcsigmaStep = (tpcsigmaMax - tpcsigmaMin) / NtpcsigmaBins;
 144 /**************************************************************/
 145 Int_t NparamsBins[kNHistoParams] = {NcentralityBins, NtofsigmaBins, NptBins, NtpcsigmaBins};
 146 Double_t *paramBin[kNHistoParams] = {centralityBin, tofsigmaBin, ptBin, tpcsigmaBin};
 147 /**************************************************************/
 148
 149 /**************************************************************/
 150 /**************************************************************/
 151 /**************************************************************/
 152
 153 /**************************************************************/
 154
 155 AliTOFGeometry tofGeo;
 156 Float_t c = TMath::C() * 1.e2 / 1.e12; /* cm/ps */
 157 Float_t c_1 = 1. / c;
 158
 159 Double_t
 160 GenerateRandomHit(TH1F *hT0Fill, Double_t t0fill, Int_t index)
 161 {
 162
 163   Int_t det[5];
 164   Float_t length, timeexp, pos[3];
 165
 166   /* compute length and expected time */
 167   tofGeo.GetVolumeIndices(index, det);
 168   tofGeo.GetPosPar(det, pos);
 169   length = 0.;
 170   for (Int_t i = 0; i < 3; i++) length += pos[i] * pos[i];
 171   length = TMath::Sqrt(length);
 172   timeexp = length * c_1;
 173
 174   Double_t hittime = hT0Fill->GetRandom() - t0fill + timeexp;
 175   return hittime;
 176
 177 }
 178
 179 /**************************************************************/
 180
 181 TOFpid(const Char_t *filename, Int_t ipart, Int_t icharge, Int_t iipart, Bool_t dorapidityCut = kTRUE, Bool_t electronCut = kFALSE, Bool_t cutOnTPC = kFALSE, Float_t tpcsignalMin = -2., Float_t tpcsignalMax = 2., Int_t evMax = kMaxInt, Int_t startEv = 0, Bool_t mcFlag = kFALSE)
 182 {
 183
 184   printf("****************************************\n");
 185   printf("RUNNING TOF PID:\n");
 186   if (dorapidityCut)
 187     printf("RAPIDITY-CUT:   %s\n", AliPID::ParticleName(ipart));
 188   else
 189     printf("NO RAPIDITY-CUT\n");
 190   printf("CHARGE:         %s\n", chargeName[icharge]);
 191   printf("PARTICLE-ID:    %s\n", AliPID::ParticleName(iipart));
 192   if (electronCut)
 193     printf("-> ELECTRON CUT REQUESTED\n");
 194   if (cutOnTPC)
 195     printf(" -> CUT-ON-TPC [%3.1f,%3.1f]\n", tpcsignalMin, tpcsignalMax);
 196   printf("****************************************\n");
 197
 198   /* include path for ACLic */
 199   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
 200   gSystem->AddIncludePath("-I$ALICE_ROOT/TOF");
 201   /* load libraries */
 202   gSystem->Load("libANALYSIS");
 203   gSystem->Load("libANALYSISalice");
 204   /* build analysis task class */
 205   gROOT->LoadMacro("AliAnalysisParticle.cxx+g");
 206   gROOT->LoadMacro("AliAnalysisEvent.cxx+g");
 207   gROOT->LoadMacro("AliAnalysisTrack.cxx+g");
 208
 209   /* create TOF response with tail */
 210   //  gROOT->LoadMacro("~/ALICE.2011/ANALYSIS/TOFSpectraPbPb/macros/TOFsignal.C");
 211   gROOT->LoadMacro("TOFsignal.C");
 212   TF1 *fTOFsignal = new TF1("fTOFsignal", TOFsignal, -2440., 2440., 4);
 213   fTOFsignal->SetParameter(0, 1.);
 214   fTOFsignal->SetParameter(1, 0.);
 215   fTOFsignal->SetParameter(2, tofReso);
 216   fTOFsignal->SetParameter(3, tofTail);
 217
 218   /* open file, get tree and connect */
 219   TFile *filein = TFile::Open(filename);
 220   TTree *treein = (TTree *)filein->Get("aodTree");
 221   printf("got \"aodTree\": %d entries\n", treein->GetEntries());
 222   AliAnalysisEvent *analysisEvent = new AliAnalysisEvent();
 223   TClonesArray *analysisTrackArray = new TClonesArray("AliAnalysisTrack");
 224   AliAnalysisTrack *analysisTrack = NULL;
 225   treein->SetBranchAddress("AnalysisEvent", &analysisEvent);
 226   treein->SetBranchAddress("AnalysisTrack", &analysisTrackArray);
 227
 228   /* open hT0fill for mismatch */
 229   TFile *filein_T0Fill = TFile::Open(t0FillOnlineFileName);
 230   TH1F *hT0Fill = (TH1F *)filein_T0Fill->Get("hT0Fill");
 231   Double_t t0fill = t0Fill_offset;
 232
 233   /* open enabled flag map */
 234   TH1F *hEnabledFlag = NULL;
 235   if (enabledChannelsFileName) {
 236     TFile *enabledfile = TFile::Open(enabledChannelsFileName);
 237     hEnabledFlag = (TH1F *)enabledfile->Get("hEnabledFlag");
 238   }
 239
 240   /**************************************************************/
 241   /*** HISTOS ***************************************************/
 242   /**************************************************************/
 243
 244   /* run-time binning */
 245   for (Int_t ibin = 0; ibin < NtofsigmaBins + 1; ibin++)
 246     tofsigmaBin[ibin] = tofsigmaMin + ibin * tofsigmaStep;
 247   for (Int_t ibin = 0; ibin < NtofsignalBins + 1; ibin++)
 248     tofsignalBin[ibin] = tofsignalMin + ibin * tofsignalStep;
 249   for (Int_t ibin = 0; ibin < NtpcsigmaBins + 1; ibin++)
 250     tpcsigmaBin[ibin] = tpcsigmaMin + ibin * tpcsigmaStep;
 251   for (Int_t ibin = 0; ibin < NtofresoBins + 1; ibin++)
 252     tofresoBin[ibin] = tofresoMin + ibin * tofresoStep;
 253
 254   /* histos */
 255   TH1F *hAcceptedEvents = new TH1F(Form("hAcceptedEvents_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin);
 256   TH2I *hAcceptedTracks = new TH2I(Form("hAcceptedTracks_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin);
 257
 258   TH3I *hTOFreso = new TH3I(Form("hTOFreso_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofresoBins, tofresoBin);
 259
 260   TH3I *hTOFpid = new TH3I(Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin);
 261   TH3I *hTOFpid_withTZERO = new TH3I(Form("hTOFpid_withTZERO_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin);
 262   TH3I *hTOFmismatch = new TH3I(Form("hTOFmismatch_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin);
 263   TH3I *hTOFexpected[AliPID::kSPECIES];
 264   for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
 265     hTOFexpected[iiipart] = new TH3I(Form("hTOFexpected_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin);
 266   }
 267
 268   TH3I *hTOFpid_delta = new TH3I(Form("hTOFpid_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsignalBins, tofsignalBin);
 269   TH3I *hTOFmismatch_delta = new TH3I(Form("hTOFmismatch_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsignalBins, tofsignalBin);
 270   TH3I *hTOFexpected_delta[AliPID::kSPECIES];
 271   for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
 272     hTOFexpected_delta[iiipart] = new TH3I(Form("hTOFexpected_delta_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsignalBins, tofsignalBin);
 273   }
 274
 275   /**************************************************************/
 276   /**************************************************************/
 277   /**************************************************************/
 278
 279   /* TOF PID response */
 280   AliTOFPIDResponse tofResponse;
 281   tofResponse.SetTimeResolution(tofReso);
 282   /* TPC PID response */
 283   AliTPCPIDResponse *tpcResponse = AliAnalysisTrack::GetTPCResponse();
 284
 285   /* start stopwatch */
 286   TStopwatch timer;
 287   timer.Start();
 288
 289   /* verbose */
 290   printf("***** RUNNING for %s %s *****\n", chargeName[icharge], AliPID::ParticleName(ipart));
 291   if (!dorapidityCut)
 292     printf("***** NO RAPIDITY CUT REQUESTED *****\n");
 293   if (cutOnTPC) {
 294     printf("***** CUT-ON-TPC requested %3.1f-%3.1f *****\n", tpcsignalMin, tpcsignalMax);
 295   }
 296
 297   /* loop over events */
 298   Bool_t hastofpid;
 299   Int_t charge, index;
 300   UShort_t dedxN;
 301   Double_t cent, p, pt, mt, tofsignal, tpcsignal;
 302   Double_t dedx, bethe, deltadedx, dedx_sigma, ptpc;
 303   Double_t time, time_sigma, timezero, timezeroTOF, timezeroA, timezeroC, timezero_sigmaTOF, timezero_sigma, tof, tof_sigma, texp, texp_sigma, deltat, deltat_sigma, tof_rnd, tof_th, signal_smear, timezero_smear, texp_smear;
 304
 305   /* open TZERO calibfile */
 306   TH1 *hCentrality_TZEROA_mean = NULL;
 307   TH1 *hCentrality_TZEROA_sigma = NULL;
 308   TH1 *hCentrality_TZEROC_mean = NULL;
 309   TH1 *hCentrality_TZEROC_sigma = NULL;
 310   TH1 *hCentrality_TOF_mean = NULL;
 311   TH1 *hCentrality_TOF_TZEROA_mean = NULL;
 312   TH1 *hCentrality_TOF_TZEROC_mean = NULL;
 313   TH1 *hCentrality_TOF_TZEROTOF_mean = NULL;
 314   TH1 *hCentrality_TZEROA_reso = NULL;
 315   TH1 *hCentrality_TZEROC_reso = NULL;
 316   TH1 *hCentrality_TZEROAND_reso = NULL;
 317   if (1) {
 318     TFile *calibfile = TFile::Open("TZEROcalibration.root");
 319     hCentrality_TZEROA_mean = (TH1 *)calibfile->Get("hCentrality_TZEROA_mean");
 320     hCentrality_TZEROA_sigma = (TH1 *)calibfile->Get("hCentrality_TZEROA_sigma");
 321     hCentrality_TZEROC_mean = (TH1 *)calibfile->Get("hCentrality_TZEROC_calib");
 322     hCentrality_TZEROC_sigma = (TH1 *)calibfile->Get("hCentrality_TZEROC_sigma");
 323     hCentrality_TOF_mean = (TH1 *)calibfile->Get("hCentrality_TOF_mean");
 324     hCentrality_TOF_TZEROA_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROA_mean");
 325     hCentrality_TOF_TZEROC_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROC_mean");
 326     hCentrality_TOF_TZEROTOF_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROTOF_mean");
 327
 328     TFile *resofile = TFile::Open("TZEROresolution.root");
 329     hCentrality_TZEROA_reso = (TH1 *)resofile->Get("hTZEROA_reso");
 330     hCentrality_TZEROC_reso = (TH1 *)resofile->Get("hTZEROC_reso");
 331     hCentrality_TZEROAND_reso = (TH1 *)resofile->Get("hTZEROAND_reso");
 332   }
 333   Double_t TZEROA_mean;
 334   Double_t TZEROA_sigma;
 335   Double_t TZEROC_mean;
 336   Double_t TZEROC_sigma;
 337   Double_t TOF_mean;
 338   Double_t TOF_TZEROA_mean;
 339   Double_t TOF_TZEROC_mean;
 340   Double_t TOF_TZEROTOF_mean;
 341   Double_t TZEROA;
 342   Double_t TZEROA_reso;
 343   Bool_t hasTZEROA;
 344   Double_t TZEROC;
 345   Double_t TZEROC_reso;
 346   Bool_t hasTZEROC;
 347   Double_t TZEROAND;
 348   Double_t TZEROAND_reso;
 349   Bool_t hasTZEROAND;
 350   Double_t TZEROTOF;
 351   Double_t TZEROTOF_reso;
 352   Bool_t hasTZEROTOF;
 353   Double_t TZEROMEAN;
 354   Double_t TZEROMEAN_weight;
 355   Double_t TZEROBEST;
 356   Double_t TZEROBEST_reso;
 357
 358   Double_t param[kNHistoParams];
 359   for (Int_t iev = startEv; iev < treein->GetEntries() && iev < evMax; iev++) {
 360     /* get event */
 361     treein->GetEvent(iev);
 362     if (iev % 10000 == 0) printf("iev = %d\n", iev);
 363     /* check event */
 364     if (!analysisEvent->AcceptEvent(acceptEventType)) continue;
 365
 366     /*** ACCEPTED EVENT ***/
 367
 368     /* apply time-zero TOF correction */
 369     //    analysisEvent->ApplyTimeZeroTOFCorrection();
 370
 371     /* get centrality */
 372     cent = analysisEvent->GetCentralityPercentile(centralityEstimator);
 373
 374     /* vertex */
 375     Double_t vertexz = analysisEvent->GetVertexZ();
 376
 377     /* fill histos */
 378     hAcceptedEvents->Fill(cent);
 379
 380     /* TZERO corrections */
 381     Int_t icent;
 382     for (icent = 0; icent < NcentralityBins; icent++)
 383       if (cent < centralityBin[icent + 1])
 384         break;
 385     TZEROA_mean = hCentrality_TZEROA_mean ? hCentrality_TZEROA_mean->GetBinContent(icent + 1) : 0.;
 386     TZEROA_sigma = hCentrality_TZEROA_sigma ? hCentrality_TZEROA_sigma->GetBinContent(icent + 1) : 1000.;
 387     TZEROC_mean  = hCentrality_TZEROC_mean ? hCentrality_TZEROC_mean->GetBinContent(icent + 1) : 0.;
 388     TZEROC_sigma = hCentrality_TZEROC_sigma ? hCentrality_TZEROC_sigma->GetBinContent(icent + 1) : 1000.;
 389
 390     TOF_mean = hCentrality_TOF_mean ? hCentrality_TOF_mean->GetBinContent(icent + 1) : 0.;
 391     TOF_TZEROA_mean = hCentrality_TOF_TZEROA_mean ? hCentrality_TOF_TZEROA_mean->GetBinContent(icent + 1) : 0.;
 392     TOF_TZEROC_mean = hCentrality_TOF_TZEROC_mean ? hCentrality_TOF_TZEROC_mean->GetBinContent(icent + 1) : 0.;
 393     TOF_TZEROTOF_mean = hCentrality_TOF_TZEROTOF_mean ? hCentrality_TOF_TZEROTOF_mean->GetBinContent(icent + 1) : 0.;
 394
 395     TZEROA_reso = hCentrality_TZEROA_reso ? hCentrality_TZEROA_reso->GetBinContent(icent + 1) : 70.;
 396     TZEROC_reso = hCentrality_TZEROC_reso ? hCentrality_TZEROC_reso->GetBinContent(icent + 1) : 70.;
 397     TZEROAND_reso = hCentrality_TZEROAND_reso ? hCentrality_TZEROAND_reso->GetBinContent(icent + 1) : 50.;
 398
 399     /* simulate inefficient TZERO */
 400     Bool_t resetTZERO = kFALSE;
 401     if (forcetimezeroineff > 0.)
 402       if (gRandom->Uniform() < forcetimezeroineff)
 403         resetTZERO = kTRUE;
 404
 405     /* loop over tracks */
 406     for (Int_t itrk = 0; itrk < analysisTrackArray->GetEntries(); itrk++) {
 407       /* get track */
 408       analysisTrack = (AliAnalysisTrack *)analysisTrackArray->At(itrk);
 409       if (!analysisTrack) continue;
 410       /* check accepted track */
 411       if (!analysisTrack->AcceptTrack()) continue;
 412       /* check rapidity */
 413       if (dorapidityCut)
 414         if (analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift > rapidityMaxCut ||
 415             analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift < rapidityMinCut) continue;
 416
 417       /* check charge */
 418       charge = analysisTrack->GetSign() > 0. ? kPositive : kNegative;
 419       if (charge != icharge) continue;
 420
 421       /*** ACCEPTED TRACK ***/
 422
 423       /* get track info */
 424       p = analysisTrack->GetP();
 425       pt = analysisTrack->GetPt();
 426
 427       /* compute track mt */
 428       mt = TMath::Sqrt(pt * pt + AliPID::ParticleMass(ipart) * AliPID::ParticleMass(ipart));
 429
 430       /* get TPC info */
 431       dedx = analysisTrack->GetTPCdEdx();
 432       dedxN = analysisTrack->GetTPCdEdxN();
 433       ptpc = analysisTrack->GetTPCmomentum();
 434
 435       /* TPC signal */
 436       bethe = tpcResponse->GetExpectedSignal(ptpc, iipart);
 437       /* fix electron signal */
 438       if (iipart == AliPID::kElectron)
 439         bethe += 23.;
 440       deltadedx = dedx - bethe;
 441       dedx_sigma = 0.07 * bethe;
 442       tpcsignal = deltadedx / dedx_sigma;
 443
 444       /* electronCut requested, remove electrons */
 445       if (electronCut) {
 446         /* TPC signal */
 447         bethe = tpcResponse->GetExpectedSignal(ptpc, AliPID::kElectron);
 448         /* fix electron signal */
 449         bethe += 23.;
 450         deltadedx = dedx - bethe;
 451         dedx_sigma = 0.07 * bethe;
 452         tpcsignal = deltadedx / dedx_sigma;
 453         if (TMath::Abs(tpcsignal) < 1.5) continue;
 454       }
 455
 456       /* cut on TPC signal if requested */
 457       if (cutOnTPC && (tpcsignal < tpcsignalMin || tpcsignal > tpcsignalMax))
 458         continue;
 459
 460       /* fill histos */
 461       hAcceptedTracks->Fill(cent, pt);
 462
 463       /* set TOF pid flag */
 464       hastofpid = analysisTrack->HasTOFPID(hEnabledFlag);
 465
 466       /* check TOF pid */
 467       if (!hastofpid)
 468         continue;
 469
 470       /*** ACCEPTED TRACK WITH TOF PID ***/
 471
 472       /* apply expected time correction */
 473       //      analysisTrack->ApplyTOFExpectedTimeCorrection();
 474
 475       /* get TOF info */
 476       index = analysisTrack->GetTOFIndex();
 477       time = analysisTrack->GetTOFTime() - TOF_mean;
 478       time_sigma = tofReso;
 479       /* TZEROTOF */
 480       TZEROTOF = analysisEvent->GetTimeZeroTOF(analysisTrack->GetP());
 481       TZEROTOF_reso = analysisEvent->GetTimeZeroTOFSigma(analysisTrack->GetP());
 482       hasTZEROTOF = TZEROTOF_reso < 150.;
 483       if (hasTZEROTOF) {
 484         TZEROTOF += TOF_TZEROTOF_mean - TOF_mean;
 485         TZEROTOF_reso *= TZEROTOF_resoScaleFactor;
 486       }
 487       /* TZEROA */
 488       TZEROA = analysisEvent->GetTimeZeroT0(1) - TZEROA_shift;
 489       //      TZEROA_reso = TZEROA_resolution[icent];
 490       hasTZEROA = TMath::Abs(TZEROA - TZEROA_mean) < 3. * TZEROA_sigma;
 491       TZEROA += TOF_TZEROA_mean - TOF_mean;
 492       /* TZEROC */
 493       TZEROC = analysisEvent->GetTimeZeroT0(2) - TZEROC_shift;
 494       //      TZEROC_reso = TZEROC_resolution[icent];
 495       hasTZEROC = TMath::Abs(TZEROC - TZEROC_mean) < 3. * TZEROC_sigma;
 496       TZEROC += TOF_TZEROC_mean - TOF_mean;
 497       /* TZEROAND */
 498       TZEROAND = (TZEROA + TZEROC) * 0.5;
 499       //      TZEROAND_reso = TZEROAND_resolution[icent];
 500       hasTZEROAND = hasTZEROA && hasTZEROC;
 501       /* TZEROMEAN */
 502       TZEROMEAN = TZEROTOF / TZEROTOF_reso / TZEROTOF_reso;
 503       TZEROMEAN_weight = 1. / TZEROTOF_reso / TZEROTOF_reso;
 504       if (hasTZEROAND) {
 505         //      printf("TZEROAND\n");
 506         TZEROMEAN += TZEROAND / TZEROAND_reso / TZEROAND_reso;
 507         TZEROMEAN_weight = 1. / TZEROAND_reso / TZEROAND_reso;
 508       }
 509       else if (hasTZEROA) {
 510         //      printf("TZEROA\n");
 511         TZEROMEAN += TZEROA / TZEROA_reso / TZEROA_reso;
 512         TZEROMEAN_weight = 1. / TZEROA_reso / TZEROA_reso;
 513       }
 514       else if (hasTZEROC) {
 515         //      printf("TZEROC\n");
 516         TZEROMEAN += TZEROC / TZEROC_reso / TZEROC_reso;
 517         TZEROMEAN_weight = 1. / TZEROC_reso / TZEROC_reso;
 518       }
 519       timezero = TZEROMEAN / TZEROMEAN_weight;
 520       timezero_sigma = TMath::Sqrt(1. / TZEROMEAN_weight);
 521       /* TZEROBEST */
 522       TZEROBEST = TZEROTOF;
 523       TZEROBEST_reso = TZEROTOF_reso;
 524       if (hasTZEROAND && TZEROAND_reso < TZEROBEST_reso) {
 525         TZEROBEST = TZEROAND;
 526         TZEROBEST_reso = TZEROAND_reso;
 527       }
 528       else if (hasTZEROA && TZEROA_reso < TZEROBEST_reso) {
 529         TZEROBEST = TZEROA;
 530         TZEROBEST_reso = TZEROA_reso;
 531       }
 532       if (hasTZEROC && TZEROC_reso < TZEROBEST_reso) {
 533         TZEROBEST = TZEROC;
 534         TZEROBEST_reso = TZEROC_reso;
 535       }
 536       timezero = TZEROBEST;
 537       timezero_sigma = TZEROBEST_reso;
 538
 539       /* DEBUG */
 540       //      timezero = 0.;//TZEROTOF;
 541       //      timezero_sigma = 203.854691;//TZEROTOF_reso;
 542
 543       //      if (timezero == 0.)
 544       //        printf("%f %f\n", timezero, timezero_sigma);
 545
 546       timezero_sigma *= scaletimezerosigma;
 547
 548       if (resetTZERO) {
 549         timezero = 0.;
 550         timezero_sigma = timezero_spread;
 551       }
 552
 553
 554       tof = time - timezero;
 555       tof_sigma = TMath::Sqrt(time_sigma * time_sigma + timezero_sigma * timezero_sigma);
 556
 557       /* TOF expected time */
 558       texp = analysisTrack->GetTOFExpTime(iipart);
 559       texp_sigma = analysisTrack->GetTOFExpTimeSigma(iipart) * scaletexpreso[iipart];
 560
 561       /* TOF signal */
 562       deltat = tof - texp;
 563       deltat_sigma = TMath::Sqrt(tof_sigma * tof_sigma + texp_sigma * texp_sigma);
 564       tofsignal = deltat / deltat_sigma;
 565
 566       /* fill histo */
 567       hTOFpid->Fill(cent, pt, tofsignal);
 568       hTOFpid_delta->Fill(cent, p, deltat);
 569       hTOFreso->Fill(cent, pt, deltat_sigma);
 570       if (hasTZEROTOF || hasTZEROA || hasTZEROC || hasTZEROAND)
 571         hTOFpid_withTZERO->Fill(cent, pt, tofsignal);
 572
 573
 574       /*** EXPECTED MISMATCH ***/
 575
 576       /* loop to generated random hits */
 577       for (Int_t irnd = 0; irnd < NmismatchTrials; irnd++) {
 578
 579         /* generate ramdom tof values */
 580         tof_rnd = GenerateRandomHit(hT0Fill, t0fill, index);
 581
 582         /* TOF signal */
 583         deltat = tof_rnd - texp;
 584         tofsignal = deltat / deltat_sigma;
 585
 586         /* fill histo */
 587         hTOFmismatch->Fill(cent, pt, tofsignal);
 588         hTOFmismatch_delta->Fill(cent, p, deltat);
 589
 590       } /* end of loop over generated random hits */
 591
 592       /*** EXPECTED SIGNALS ***/
 593
 594       /* loop over other particles */
 595       for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
 596
 597         /* generate expected tof value */
 598         tof_th = analysisTrack->GetTOFExpTime(iiipart);
 599         texp_sigma = analysisTrack->GetTOFExpTimeSigma(iiipart) * scaletexpreso[iiipart];
 600
 601         /* loop over many trials */
 602         for (Int_t irnd = 0; irnd < NexpectedTrials; irnd++) {
 603
 604           /* tof response smearing */
 605           signal_smear = fTOFsignal->GetRandom();
 606           /* timezero resolution smearing */
 607           timezero_smear = gRandom->Gaus(0., timezero_sigma);
 608           /* texp resolution smearing */
 609           texp_smear = gRandom->Gaus(0., texp_sigma);
 610
 611           /* deltat and sigma */
 612           deltat = tof_th - texp + signal_smear + timezero_smear + texp_smear;
 613           tofsignal = deltat / deltat_sigma;
 614
 615           /* fill histo */
 616           hTOFexpected[iiipart]->Fill(cent, pt, tofsignal);
 617           hTOFexpected_delta[iiipart]->Fill(cent, p, deltat);
 618
 619         } /* end of loop over many trials */
 620       } /* end of loop over other particle */
 621     } /* end of loop over tracks */
 622   } /* end of loop over events */
 623
 624   /* stop stopwatch */
 625   timer.Stop();
 626   timer.Print();
 627
 628   /* output */
 629   TString outputstring = "TOFpid";
 630   if (!dorapidityCut)
 631     outputstring += "_noRapidityCut";
 632   if (electronCut)
 633     outputstring += "_electronCut";
 634   if (cutOnTPC)
 635     outputstring += Form("_cutOnTPC[%3.1f,%3.1f]", , tpcsignalMin, tpcsignalMax);
 636   outputstring += Form("_%s_%s_%sID.%s", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), filename);
 637   TFile *fileout = TFile::Open(outputstring.Data(), "RECREATE");
 638   hAcceptedEvents->Write();
 639   hAcceptedTracks->Write();
 640   hTOFpid->Write();
 641   hTOFreso->Write();
 642   hTOFpid_withTZERO->Write();
 643   hTOFmismatch->Write();
 644   for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
 645     hTOFexpected[iiipart]->Write();
 646   hTOFpid_delta->Write();
 647   hTOFmismatch_delta->Write();
 648   for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
 649     hTOFexpected_delta[iiipart]->Write();
 650
 651   fileout->Close();
 652
 653 }
 654
 655 //___________________________________________________________________________________
 656
 657 TOFspectra_defaultFit(const Char_t *filename)
 658 {
 659
 660 Bool_t EXPECTED_SIGNAL_TEMPLATE = kTRUE;
 661 Bool_t EXPECTED_SIGNAL_FIT = kFALSE;
 662 Bool_t EXPECTED_BACKGROUND_TEMPLATE = kFALSE;
 663 Bool_t EXPECTED_BACKGROUND_FIT = kTRUE;
 664
 665 }
 666
 667 //___________________________________________________________________________________
 668
 669 TOFspectra_defaultFit_fitElectrons(const Char_t *filename, Float_t electronLimit = 5.)
 670 {
 671
 672
 673   TOFspectra(filename, electronLimit);
 674 }
 675
 676 //___________________________________________________________________________________
 677
 678 TOFspectra_signalFit(Bool_t fixParams = kTRUE, Float_t scaleSigma = 1., Float_t scaleTail = 1.)
 679 {
 680
 681   EXPECTED_SIGNAL_TEMPLATE = kFALSE;
 682   EXPECTED_SIGNAL_FIT = kTRUE;
 683   FIX_SIGNAL_MEAN = fixParams;
 684   FIX_SIGNAL_SIGMA = fixParams;
 685   FIX_SIGNAL_TAIL = fixParams;
 686   SCALE_SIGNAL_SIGMA = scaleSigma;
 687   SCALE_SIGNAL_TAIL = scaleTail;
 688
 689 }
 690
 691 //___________________________________________________________________________________
 692
 693 TOFspectra_bkgFit(Bool_t fixParams = kTRUE, Float_t scaleSigma = 1., Float_t scaleTail = 1.)
 694 {
 695
 696   EXPECTED_BACKGROUND_TEMPLATE = kFALSE;
 697   EXPECTED_BACKGROUND_FIT = kTRUE;
 698   FIX_BACKGROUND_MEAN = fixParams;
 699   FIX_BACKGROUND_SIGMA = fixParams;
 700   FIX_BACKGROUND_TAIL = fixParams;
 701   SCALE_BACKGROUND_SIGMA = scaleSigma;
 702   SCALE_BACKGROUND_TAIL = scaleTail;
 703
 704   TOFspectra(filename);
 705
 706 }
 707
 708 //___________________________________________________________________________________
 709
 710 void
 711 TOFspectra(const Char_t *filename, Float_t electronLimit = 0.)
 712 {
 713
 714   for (Int_t icent = 0; icent < NcentralityBins; icent++)
 715     for (Int_t icharge = 0; icharge < kNCharges; icharge++)
 716       for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
 717         TOFspectrum(filename, ipart, icharge, ipart, icent, -1., -1., electronLimit);
 718
 719 }
 720
 721 //___________________________________________________________________________________
 722
 723 #define USENSIGMA 1
 724
 725 /* fit ranges */
 726 Double_t fitPtMin[AliPID::kSPECIES] = {0.5, 0.5, 0.3, 0.4, 0.5};
 727 Double_t fitPtMax[AliPID::kSPECIES] = {3.0, 3.0, 3.0, 3.0, 5.0};
 728 #if USENSIGMA
 729 Double_t fitSigmaMin[AliPID::kSPECIES] = {tofsigmaMin, tofsigmaMin, -25., -75., -65.};
 730 Double_t fitSigmaMax[AliPID::kSPECIES] = {tofsigmaMax, tofsigmaMax, 225., 200., 100.};
 731 #else
 732 Double_t fitSigmaMin[AliPID::kSPECIES] = {-24400., -24400., -24400., -24400., -24400.};
 733 Double_t fitSigmaMax[AliPID::kSPECIES] = {24400., 24400., 24400., 24400., 24400.};
 734 #endif
 735
 736 void
 737 TOFspectrum(const Char_t *filename, Int_t ipart, Int_t icharge, Int_t iipart, Int_t icent, Float_t ptMin = -1., Float_t ptMax = -1., Bool_t checkHistoFlag = kFALSE)
 738 {
 739
 740   printf("****************************************\n");
 741   printf("RUNNING SPECTRA FIT:\n");
 742   printf("RAPIDITY-CUT:   %s\n", AliPID::ParticleName(ipart));
 743   printf("RAPIDITY-CUT:   %s\n", AliPID::ParticleName(ipart));
 744   printf("CHARGE:         %s\n", chargeName[icharge]);
 745   printf("PARTICLE:       %s\n", AliPID::ParticleName(iipart));
 746   printf("CENTRALITY BIN: %d\n", icent);
 747   printf("****************************************\n");
 748
 749   /* mismatch balance functions */
 750   TF1 *fMismatchBalanceFunction[5][2];
 751   Double_t mismatchBalanceParam0[5][2] = {
 752     0., 0.,
 753     0., 0.,
 754     0.669488, 0.599374,
 755     -1.2582, -2.53613,
 756     5.48850e-01, 10.1622
 757   };
 758   Double_t mismatchBalanceParam1[5][2] = {
 759     0., 0.,
 760     0., 0.,
 761     -1.64894, -0.825764,
 762     -1.63683, -2.55156,
 763     -1.64894, -7.09531
 764   };
 765   for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
 766     for (Int_t iiicharge = 0; iiicharge < 2; iiicharge++) {
 767       fMismatchBalanceFunction[iiipart][iiicharge] = new TF1(Form(""), "1. - [0]*TMath::Exp(x*[1])", 0., 5.);
 768       fMismatchBalanceFunction[iiipart][iiicharge]->SetParameter(0, mismatchBalanceParam0[iiipart][iiicharge]);
 769       fMismatchBalanceFunction[iiipart][iiicharge]->SetParameter(1, mismatchBalanceParam1[iiipart][iiicharge]);
 770     }
 771   }
 772
 773   /* open data */
 774   TFile *filein = TFile::Open(filename);
 775
 776   /* get number of events */
 777   TH1F *hAcceptedEvents = (TH1F *)filein->Get(Form("hAcceptedEvents_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 778   if (!hAcceptedEvents) {
 779     printf("cannot find %s\n", Form("hAcceptedEvents_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 780     return;
 781   }
 782   Double_t nevents;
 783   if (icent < 0 || icent >= NcentralityBins)
 784     nevents = hAcceptedEvents->GetEntries();
 785   else
 786     nevents = hAcceptedEvents->Integral(icent + 1, icent + 1);
 787   printf("N_EVENTS      : %d\n", nevents);
 788   printf("****************************************\n");
 789
 790   /* get histos */
 791   TH2I *hAcceptedTracks = (TH2I *)filein->Get(Form("hAcceptedTracks_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 792   if (!hAcceptedTracks) {
 793     printf("cannot find %s\n", Form("hAcceptedTracks_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 794     //    return;
 795   }
 796 #if USENSIGMA
 797   TH3I *hTOFpid = (TH3I *)filein->Get(Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 798 #else
 799   TH3I *hTOFpid = (TH3I *)filein->Get(Form("hTOFpid_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 800 #endif
 801   if (!hTOFpid) {
 802     printf("cannot find %s\n", Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 803     return;
 804   }
 805 #if USENSIGMA
 806   TH3I *hTOFmismatch = (TH3I *)filein->Get(Form("hTOFmismatch_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 807 #else
 808   TH3I *hTOFmismatch = (TH3I *)filein->Get(Form("hTOFmismatch_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 809 #endif
 810   if (!hTOFmismatch) {
 811     printf("cannot find %s\n", Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
 812     return;
 813   }
 814   TH3I *hTOFexpected[AliPID::kSPECIES];
 815   for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
 816 #if USENSIGMA
 817     hTOFexpected[iiipart] = (TH3I *)filein->Get(Form("hTOFexpected_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)));
 818 #else
 819     hTOFexpected[iiipart] = (TH3I *)filein->Get(Form("hTOFexpected_delta_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)));
 820 #endif
 821     if (!hTOFexpected[iiipart]) {
 822       printf("cannot find %s\n", Form("hTOFexpected_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)));
 823       return;
 824     }
 825   }
 826
 827   /* setup centrality range */
 828   if (icent < 0 || icent >= NcentralityBins) {
 829     printf("WARNING: undefined centrality -> using 00-90\% range\n");
 830     if (hAcceptedTracks) hAcceptedTracks->GetXaxis()->SetRange(1, NcentralityBins);
 831     hTOFpid->GetXaxis()->SetRange(1, NcentralityBins);
 832     hTOFmismatch->GetXaxis()->SetRange(1, NcentralityBins);
 833     for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
 834       hTOFexpected[iiipart]->GetXaxis()->SetRange(1, NcentralityBins);
 835   }
 836   else {
 837     printf("***** FITTING CENTRALITY-BIN [%02d, %02d] %% *****\n", centralityBin[icent], centralityBin[icent + 1]);
 838     if (hAcceptedTracks) hAcceptedTracks->GetXaxis()->SetRange(icent + 1, icent + 1);
 839     hTOFpid->GetXaxis()->SetRange(icent + 1, icent + 1);
 840     hTOFmismatch->GetXaxis()->SetRange(icent + 1, icent + 1);
 841     for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
 842       hTOFexpected[iiipart]->GetXaxis()->SetRange(icent + 1, icent + 1);
 843   }
 844
 845   /* init flags */
 846   Bool_t requestedRange = kFALSE;
 847   Bool_t fitElectrons = kTRUE;
 848   Bool_t fitMuons = kTRUE;
 849   Bool_t fitPions = kTRUE;
 850
 851   /* setup pt range if requested */
 852   if (ptMin > -0.001 && ptMax > -0.001 && ptMax > ptMin) {
 853     printf("***** FITTING PT-BIN [%f, %f] GeV/c *****\n", ptMin, ptMax);
 854     requestedRange = kTRUE;
 855
 856     /* check electron-fit is allowed */
 857     fitElectrons = kTRUE;
 858     if ((ptMin + 0.001) < FIT_ELECTRONS_PT_MIN || (ptMax - 0.001) > FIT_ELECTRONS_PT_MAX)
 859       fitElectrons = kFALSE;
 860     /* check muon-fit is allowed */
 861     fitMuons = kTRUE;
 862     if ((ptMin + 0.001) < FIT_MUONS_PT_MIN || (ptMax - 0.001) > FIT_MUONS_PT_MAX)
 863       fitMuons = kFALSE;
 864     /* check pion-fit is allowed */
 865     fitPions = kTRUE;
 866     if ((ptMin + 0.001) < FIT_PIONS_PT_MIN || (ptMax - 0.001) > FIT_PIONS_PT_MAX)
 867       fitPions = kFALSE;
 868
 869
 870     hTOFpid->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
 871     hTOFmismatch->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
 872     for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
 873       hTOFexpected[iiipart]->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
 874   }
 875
 876   /* output */
 877   Char_t outfilename[1024];
 878   if (icent < 0 || icent >= NcentralityBins)
 879     sprintf(outfilename, "TOFspectrum_cent0090_%s_%s_%sID.root", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart));
 880   else {
 881     sprintf(outfilename, "TOFspectrum_cent%02d%02d_%s_%s_%sID.root", centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart));
 882   }
 883   TFile *fileout = TFile::Open(outfilename, "RECREATE");
 884   TDirectory *fitDir = fileout->mkdir("FitParams");
 885   /* canvas */
 886   TCanvas *canvas = new TCanvas("canvas");
 887   canvas->SetLogy();
 888   /* histo */
 889   TH1D *hFitParamHisto[kNFitParams];
 890   for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
 891     hFitParamHisto[iparam] = new TH1D(Form("h%s", fitParamName[iparam]), fitParamTitle[iparam], NptBins, ptBin);
 892
 893   /* loop over ptBins */
 894   for (Int_t ipt = 0; ipt < NptBins; ipt++) {
 895
 896     if (!requestedRange) {
 897       if ((ptBin[ipt] + 0.001) < fitPtMin[ipart] || (ptBin[ipt + 1] - 0.001) > fitPtMax[ipart]) continue;
 898       printf("***** FITTING PT-BIN [%f, %f] GeV/c *****\n", ptBin[ipt], ptBin[ipt + 1]);
 899
 900       /* check electron-fit is allowed */
 901       fitElectrons = kTRUE;
 902       if ((ptBin[ipt] + 0.001) < FIT_ELECTRONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_ELECTRONS_PT_MAX)
 903         fitElectrons = kFALSE;
 904       /* check muon-fit is allowed */
 905       fitMuons = kTRUE;
 906       if ((ptBin[ipt] + 0.001) < FIT_MUONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_MUONS_PT_MAX)
 907         fitMuons = kFALSE;
 908       /* check pion-fit is allowed */
 909       fitPions = kTRUE;
 910       if ((ptBin[ipt] + 0.001) < FIT_PIONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_PIONS_PT_MAX)
 911         fitPions = kFALSE;
 912
 913       hTOFpid->GetYaxis()->SetRange(ipt + 1, ipt + 1);
 914       hTOFmismatch->GetYaxis()->SetRange(ipt + 1, ipt + 1);
 915       for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
 916         hTOFexpected[iiipart]->GetYaxis()->SetRange(ipt + 1, ipt + 1);
 917     }
 918
 919     /* nsigma projections */
 920     TH1 *hSignal_py = hTOFpid->Project3D("z");
 921     TH1 *hMismatch_py = hTOFmismatch->Project3D("z");
 922     TH1 *hSignalExp_py[AliPID::kSPECIES];
 923     TH1 *hSignalExpTail_py[AliPID::kSPECIES];
 924     for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
 925       hSignalExp_py[iiipart] = hTOFexpected[iiipart]->Project3D("z");
 926       hSignalExpTail_py[iiipart] = hTOFexpected[iiipart]->Project3D("z");
 927     }
 928
 929     /* prepare histos for the fitter */
 930     Int_t partbkg1[AliPID::kSPECIES] = {AliPID::kKaon, AliPID::kKaon, AliPID::kKaon, AliPID::kPion, AliPID::kPion};
 931     Int_t partbkg2[AliPID::kSPECIES] = {AliPID::kProton, AliPID::kProton, AliPID::kProton, AliPID::kProton, AliPID::kKaon};
 932     Int_t partbkg3[AliPID::kSPECIES] = {AliPID::kPion, AliPID::kPion, AliPID::kElectron, AliPID::kElectron, AliPID::kElectron};
 933     Int_t partbkg4[AliPID::kSPECIES] = {AliPID::kMuon, AliPID::kElectron, AliPID::kMuon, AliPID::kMuon, AliPID::kMuon};
 934     TH1 *hSigExp_py, *hBkgExp1_py, *hBkgExp2_py, *hBkgExp3_py, *hBkgExp4_py;
 935     hSigExp_py = EXPECTED_SIGNAL_TAIL ? hSignalExpTail_py[iipart] : hSignalExp_py[iipart];
 936     hBkgExp1_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg1[iipart]] : hSignalExp_py[partbkg1[iipart]];
 937     hBkgExp2_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg2[iipart]] : hSignalExp_py[partbkg2[iipart]];
 938     hBkgExp3_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg3[iipart]] : hSignalExp_py[partbkg3[iipart]];
 939     hBkgExp4_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg4[iipart]] : hSignalExp_py[partbkg4[iipart]];
 940
 941     /* check histos if requested */
 942     if (checkHistoFlag) {
 943       TCanvas *cCheckHisto = new TCanvas("cCheckHisto");
 944       cCheckHisto->Divide(2, 3);
 945       cCheckHisto->cd(1);
 946       hSignal_py->Draw();
 947       cCheckHisto->cd(2);
 948       hSigExp_py->Draw();
 949       cCheckHisto->cd(3);
 950       hBkgExp1_py->Draw();
 951       cCheckHisto->cd(4);
 952       hBkgExp2_py->Draw();
 953       cCheckHisto->cd(5);
 954       hBkgExp3_py->Draw();
 955       cCheckHisto->cd(6);
 956       hMismatch_py->Draw();
 957       return;
 958     }
 959
 960     Double_t rangeMin = fitSigmaMin[iipart], rangeMax = fitSigmaMax[iipart];
 961     Bool_t constrainSignal = kFALSE;
 962     Bool_t constrainBkg1 = kFALSE;
 963     Bool_t constrainBkg2 = kFALSE;
 964     Bool_t forceGaussianSignal = kFALSE;
 965     Bool_t fitBkg1 = kTRUE, fitBkg2 = kTRUE, fitBkg3 = kTRUE, fitBkg4 = kTRUE;
 966
 967     /* check whether we can fit electrons */
 968     if (!fitElectrons) {
 969       printf("INHIBIT FIT ELECTRONS\n");
 970       if (partbkg1[iipart] == AliPID::kElectron) fitBkg1 = kFALSE;
 971       if (partbkg2[iipart] == AliPID::kElectron) fitBkg2 = kFALSE;
 972       if (partbkg3[iipart] == AliPID::kElectron) fitBkg3 = kFALSE;
 973       if (partbkg4[iipart] == AliPID::kElectron) fitBkg4 = kFALSE;
 974     }
 975     /* check whether we can fit muons */
 976     if (!fitMuons) {
 977       printf("INHIBIT FIT MUONS\n");
 978       if (partbkg1[iipart] == AliPID::kMuon) fitBkg1 = kFALSE;
 979       if (partbkg2[iipart] == AliPID::kMuon) fitBkg2 = kFALSE;
 980       if (partbkg3[iipart] == AliPID::kMuon) fitBkg3 = kFALSE;
 981       if (partbkg4[iipart] == AliPID::kMuon) fitBkg4 = kFALSE;
 982     }
 983     /* check whether we can fit pions */
 984     if (!fitPions) {
 985       printf("INHIBIT FIT PIONS\n");
 986       if (partbkg1[iipart] == AliPID::kPion) fitBkg1 = kFALSE;
 987       if (partbkg2[iipart] == AliPID::kPion) fitBkg2 = kFALSE;
 988       if (partbkg3[iipart] == AliPID::kPion) fitBkg3 = kFALSE;
 989       if (partbkg4[iipart] == AliPID::kPion) fitBkg4 = kFALSE;
 990     }
 991
 992
 993     /* fit */
 994     Double_t param[kNFitParams];
 995     Double_t param_err[kNFitParams];
 996     TOFpid_fit(hSignal_py, hSigExp_py, hBkgExp1_py, hBkgExp2_py, hBkgExp3_py, hBkgExp4_py, hMismatch_py, rangeMin, rangeMax, fitBkg1, fitBkg2, fitBkg3, fitBkg4, constrainSignal, constrainBkg1, constrainBkg2, forceGaussianSignal, param, param_err, canvas);
 997
 998
 999     /* if fitting pions add electrons and to the signal */
1000     if (iipart == 2) {
1001       param[kSignalCounts] += param[kBkg3Counts] + param[kBkg4Counts];
1002       param_err[kSignalCounts] = TMath::Sqrt(param_err[kSignalCounts] * param_err[kSignalCounts] + param_err[kBkg3Counts] * param_err[kBkg3Counts] + param_err[kBkg4Counts] * param_err[kBkg3Counts]);
1003     }
1004     /* else add to bkg1 (pions) */
1005     else {
1006       param[kBkg1Counts] += param[kBkg3Counts] + param[kBkg4Counts];
1007       param_err[kBkg1Counts] = TMath::Sqrt(param_err[kBkg1Counts] * param_err[kBkg1Counts] + param_err[kBkg3Counts] * param_err[kBkg3Counts] + param_err[kBkg4Counts] * param_err[kBkg3Counts]);
1008     }
1009
1010     /* check requested pt-range */
1011     if (requestedRange)
1012       return;
1013
1014     /* write canvas */
1015     fitDir->cd();
1016     canvas->Write(Form("fitDisplay_ptBin_%3.2f_%3.2f", ptBin[ipt], ptBin[ipt + 1]));
1017
1018     /* set histo */
1019     for (Int_t iparam = 0; iparam < kNFitParams; iparam++) {
1020       hFitParamHisto[iparam]->SetBinContent(ipt + 1, param[iparam]);
1021       hFitParamHisto[iparam]->SetBinError(ipt + 1, param_err[iparam]);
1022     }
1023
1024     /* delete */
1025     delete hSignal_py;
1026     delete hMismatch_py;
1027     for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
1028       delete hSignalExp_py[iiipart];
1029       //      delete hSignalExpTail_py[iiipart];
1030     }
1031   }
1032
1033   /* check requested pt-range */
1034   if (requestedRange)
1035     return;
1036
1037   /*** POST-ANALYSIS ***/
1038
1039   TDirectory *postDir = fileout->mkdir("PostAnalysis");
1040
1041   /* compute overflows */
1042   TH1D *hOverflowCounts = new TH1D(*hFitParamHisto[kTotalCounts]);
1043   hOverflowCounts->SetNameTitle("hOverflowCounts", "Overflow counts: TotalCounts - IntegralCounts;p_{T} (GeV/c);counts");
1044   hOverflowCounts->Add(hFitParamHisto[kIntegralCounts], -1.);
1045   /* compute total mismatch counts */
1046   TH1D *hTotalMismatchCounts = new TH1D(*hFitParamHisto[kMismatchCounts]);
1047   hTotalMismatchCounts->SetNameTitle("hTotalMismatchCounts", "Total mismatch counts: MismatchCounts + OverflowCounts;p_{T} (GeV/c);counts");
1048   hTotalMismatchCounts->Add(hOverflowCounts);
1049   /* computed mismatch fraction */
1050   TH1D *hTotalMismatchFraction = new TH1D(*hTotalMismatchCounts);
1051   hTotalMismatchFraction->SetNameTitle("hTotalMismatchFraction", "Total mismatch fraction: TotalMismatchCounts / TotalCounts;p_{T} (GeV/c);");
1052   hTotalMismatchFraction->Divide(hFitParamHisto[kTotalCounts]);
1053   /* compute identified counts */
1054   TH1D *hIdentifiedCounts = new TH1D(*hFitParamHisto[kSignalCounts]);
1055   hIdentifiedCounts->SetNameTitle("hIdentifiedCounts", "Identified counts: SignalCounts + sum(BkgCounts);p_{T} (GeV/c);counts");
1056   hIdentifiedCounts->Add(hFitParamHisto[kBkg1Counts]);
1057   hIdentifiedCounts->Add(hFitParamHisto[kBkg2Counts]);
1058   //  hIdentifiedCounts->Add(hFitParamHisto[kBkg3Counts]);
1059   //  hIdentifiedCounts->Add(hFitParamHisto[kBkg4Counts]);
1060   /* compute signal fraction */
1061   TH1D *hSignalFraction = new TH1D(*hFitParamHisto[kSignalCounts]);
1062   hSignalFraction->SetNameTitle("hSignalFraction", "Signal fraction: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1063   hSignalFraction->Divide(hSignalFraction, hIdentifiedCounts, 1., 1., "B");
1064   /* compute bkg1 fraction */
1065   TH1D *hBkg1Fraction = new TH1D(*hFitParamHisto[kBkg1Counts]);
1066   hBkg1Fraction->SetNameTitle("hBkg1Fraction", "Bkg1 fraction: Bkg1Counts / IdentifiedCounts;p_{T} (GeV/c);");
1067   hBkg1Fraction->Divide(hBkg1Fraction, hIdentifiedCounts, 1., 1., "B");
1068   /* compute bkg2 fraction */
1069   TH1D *hBkg2Fraction = new TH1D(*hFitParamHisto[kBkg2Counts]);
1070   hBkg2Fraction->SetNameTitle("hBkg2Fraction", "Bkg2 fraction: Bkg2Counts / IdentifiedCounts;p_{T} (GeV/c);");
1071   hBkg2Fraction->Divide(hBkg2Fraction, hIdentifiedCounts, 1., 1., "B");
1072   /* compute bkg3 fraction */
1073   TH1D *hBkg3Fraction = new TH1D(*hFitParamHisto[kBkg3Counts]);
1074   hBkg3Fraction->SetNameTitle("hBkg3Fraction", "Bkg3 fraction: Bkg3Counts / IdentifiedCounts;p_{T} (GeV/c);");
1075   hBkg3Fraction->Divide(hBkg3Fraction, hIdentifiedCounts, 1., 1., "B");
1076   /* compute bkg4 fraction */
1077   TH1D *hBkg4Fraction = new TH1D(*hFitParamHisto[kBkg4Counts]);
1078   hBkg4Fraction->SetNameTitle("hBkg4Fraction", "Bkg4 fraction: Bkg4Counts / IdentifiedCounts;p_{T} (GeV/c);");
1079   hBkg4Fraction->Divide(hBkg4Fraction, hIdentifiedCounts, 1., 1., "B");
1080   /* rework */
1081   TH1D *hSignalFractionReworked = new TH1D(*hFitParamHisto[kSignalCounts]);
1082   hSignalFractionReworked->SetNameTitle("hSignalFractionReworked", "Signal fraction reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1083   hSignalFractionReworked->Multiply(fMismatchBalanceFunction[iipart][icharge]);
1084   TH1D *hBkg1FractionReworked = new TH1D(*hFitParamHisto[kBkg1Counts]);
1085   hBkg1FractionReworked->SetNameTitle("hBkg1FractionReworked", "Bkg1 fraction reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1086   hBkg1FractionReworked->Multiply(fMismatchBalanceFunction[partbkg1[iipart]][icharge]);
1087   TH1D *hBkg2FractionReworked = new TH1D(*hFitParamHisto[kBkg2Counts]);
1088   hBkg2FractionReworked->SetNameTitle("hBkg2FractionReworked", "Bkg1 fraction reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1089   hBkg2FractionReworked->Multiply(fMismatchBalanceFunction[partbkg2[iipart]][icharge]);
1090   TH1D *hIdentifiedCountsReworked = new TH1D(*hSignalFractionReworked);
1091   hIdentifiedCountsReworked->SetNameTitle("hIdentifiedCountsReworked", "Identified counts reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1092   hIdentifiedCountsReworked->Add(hBkg1FractionReworked);
1093   hIdentifiedCountsReworked->Add(hBkg2FractionReworked);
1094
1095   hSignalFractionReworked->Divide(hSignalFractionReworked, hIdentifiedCountsReworked, 1., 1., "B");
1096
1097   /* compute mismatch-correction counts */
1098   TH1D *hMismatchCorrectionCounts = new TH1D(*hTotalMismatchCounts);
1099   hMismatchCorrectionCounts->SetNameTitle("hMismatchCorrectionCounts", "Mismatch-correction counts: TotalMismatchCounts * SignalFraction;p_{T} (GeV/c);counts");
1100   hMismatchCorrectionCounts->Multiply(hSignalFractionReworked);
1101   /* compute mismatch-corrected signal counts */
1102   TH1D *hMismatchCorrectedSignalCounts = new TH1D(*hFitParamHisto[kSignalCounts]);
1103   hMismatchCorrectedSignalCounts->SetNameTitle("hMismatchCorrectedSignalCounts", "Mismatch-corrected signal counts: SignalCounts + MismatchCorrectionCounts;p_{T} (GeV/c);counts");
1104   hMismatchCorrectedSignalCounts->Add(hMismatchCorrectionCounts);
1105
1106
1107   /* accepted tracks histo */
1108   if (hAcceptedTracks) {
1109     TH1D *hAcceptedTracks_py = hAcceptedTracks->ProjectionY();
1110     hAcceptedTracks_py->SetNameTitle("hAcceptedTracks", "Accepted tracks;p_{T} (GeV/c);");
1111     hAcceptedTracks_py->Sumw2();
1112   }
1113
1114   /*** RAW SPECTRA ***/
1115
1116   TDirectory *rawDir = fileout->mkdir("RawSpectra");
1117
1118   /* compute normalized raw yield */
1119   TH1D *hNormalizedRawYield = new TH1D(*hFitParamHisto[kSignalCounts]);
1120   hNormalizedRawYield->SetNameTitle("hNormalizedRawYield", "Raw yield;p_{T} (GeV/c);#frac{1}{N_{ev}} #frac{d^{2}N}{dy dp_{T}}");
1121   TOFpid_normalize(hNormalizedRawYield, nevents);
1122   /* compute normalized mismatch-corrected raw yield */
1123   TH1D *hNormalizedMismatchCorrectedRawYield = new TH1D(*hMismatchCorrectedSignalCounts);
1124   hNormalizedMismatchCorrectedRawYield->SetNameTitle("hNormalizedMismatchCorrectedRawYield", "Mismatch-corrected raw yield;p_{T} (GeV/c);#frac{1}{N_{ev}} #frac{d^{2}N}{dy dp_{T}}");
1125   TOFpid_normalize(hNormalizedMismatchCorrectedRawYield, nevents);
1126
1127   /*** OUTPUT ***/
1128
1129   /* write fir params histos */
1130   fitDir->cd();
1131   for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
1132     hFitParamHisto[iparam]->Write();
1133   /* write post-analysis histos */
1134   postDir->cd();
1135   hOverflowCounts->Write();
1136   hTotalMismatchCounts->Write();
1137   hTotalMismatchFraction->Write();
1138   hIdentifiedCounts->Write();
1139   hSignalFraction->Write();
1140   hSignalFractionReworked->Write();
1141   hBkg1Fraction->Write();
1142   hBkg2Fraction->Write();
1143   hBkg3Fraction->Write();
1144   hBkg4Fraction->Write();
1145   hMismatchCorrectionCounts->Write();
1146   hMismatchCorrectedSignalCounts->Write();
1147   if (hAcceptedTracks) hAcceptedTracks_py->Write();
1148   /* write raw spectra histos */
1149   rawDir->cd();
1150   hNormalizedRawYield->Write();
1151   hNormalizedMismatchCorrectedRawYield->Write();
1152
1153   /* clean up */
1154   delete canvas;
1155   for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
1156     delete hFitParamHisto[iparam];
1157   delete hOverflowCounts;
1158   delete hTotalMismatchCounts;
1159   delete hTotalMismatchFraction;
1160   delete hIdentifiedCounts;
1161   delete hSignalFraction;
1162   delete hBkg1Fraction;
1163   delete hBkg2Fraction;
1164   delete hBkg3Fraction;
1165   delete hBkg4Fraction;
1166   delete hMismatchCorrectionCounts;
1167   delete hMismatchCorrectedSignalCounts;
1168   if (hAcceptedTracks) {
1169     delete hAcceptedTracks_py;
1170   }
1171   delete hNormalizedRawYield;
1172   delete hNormalizedMismatchCorrectedRawYield;
1173
1174   /* close file */
1175   fileout->Close();
1176 }
1177
1178 //___________________________________________________________________________________
1179
1180 Float_t
1181 TOFpid_histomin(TH1 *h)
1182 {
1183
1184   for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++)
1185     if (h->GetBinContent(ibin + 1) > 0.)
1186       return h->GetXaxis()->GetBinCenter(ibin + 1);
1187   return kMaxInt;
1188 }
1189
1190 //___________________________________________________________________________________
1191
1192 Float_t
1193 TOFpid_histomax(TH1 *h)
1194 {
1195
1196   for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--)
1197     if (h->GetBinContent(ibin) > 0.)
1198       return h->GetXaxis()->GetBinCenter(ibin);
1199   return -kMaxInt;
1200 }
1201
1202 //___________________________________________________________________________________
1203
1204 TOFpid_checkneg(TH1 *h)
1205 {
1206
1207   for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++)
1208     if (h->GetBinContent(ibin + 1) <= 0.) {
1209       h->SetBinContent(ibin + 1, 1.e-300);
1210       //      h->SetBinError(ibin + 1, 0.1);
1211     }
1212 }
1213
1214 //___________________________________________________________________________________
1215
1216 Double_t
1217 TOFpid_fit(TH1 *hSignal, TH1 *hSigExp, TH1 *hBkgExp1, TH1 *hBkgExp2, TH1 *hBkgExp3, TH1 *hBkgExp4, TH1 *hMismatch, Double_t rangeMin, Double_t rangeMax, Bool_t fitBkg1, Bool_t fitBkg2, Bool_t fitBkg3, Bool_t fitBkg4, Bool_t constrainSignal, Bool_t constrainBkg1, Bool_t constrainBkg2, Bool_t forceGaussianSignal, Double_t *param = NULL, Double_t *param_err = NULL, TCanvas *canvas = NULL)
1218 {
1219
1220   /** ROOFIT ***/
1221   gSystem->Load("libRooFit");
1222   using namespace RooFit;
1223  /*** LOAD GAUSSIANTAIL CLASS ***/
1224   gSystem->Load("RooFermiCutoff_cxx");
1225   gSystem->Load("RooGaussianTail_cxx");
1226
1227   /*** DEFINE FIT RANGE ***/
1228
1229   printf("***** FIT RANGE DEFINITION *****\n");
1230
1231   /* check mismatch histogram to define min/max fit range */
1232   //  rangeMin = TMath::Max(rangeMin, TOFpid_histomin(hMismatch));
1233   //  rangeMax = TMath::Min(rangeMax, TOFpid_histomax(hMismatch));
1234   /* fix zeroes */
1235   TOFpid_checkneg(hMismatch);
1236
1237   /* define range */
1238   RooRealVar x("x", "n_{#sigma}", 0., rangeMin, rangeMax, "");
1239   printf("FIT RANGE DEFINED: %f -> %f\n", rangeMin, rangeMax);
1240   printf("********************************\n");
1241
1242   /*** DEFINE HISTOGRAM DATA ***/
1243
1244   /* define data to fit and background from input histogram */
1245   RooDataHist hdata("hdata", "hdata", x, hSignal);
1246   RooDataHist hsig("hsig", "hsig", x, hSigExp);
1247   RooDataHist hbkg1("hbkg1", "hbkg1", x, hBkgExp1);
1248   RooDataHist hbkg2("hbkg2", "hbkg2", x, hBkgExp2);
1249   RooDataHist hbkg3("hbkg3", "hbkg3", x, hBkgExp3);
1250   RooDataHist hbkg4("hbkg4", "hbkg4", x, hBkgExp4);
1251   RooDataHist hmismatch("hmismatch", "hmismatch", x, hMismatch);
1252
1253   /*** DEFINE SIGNAL SHAPE ***/
1254
1255   printf("***** SIGNAL SHAPE DEFINITION *****\n");
1256
1257   /* variables */
1258   if (FIX_SIGNAL_MEAN) {
1259     RooRealVar mean("mean", "mean", DEFAULT_SIGNAL_MEAN, "");
1260     printf("FIXED SIGNAL_MEAN = %f\n", mean.getVal());
1261   }
1262   else {
1263     RooRealVar mean("mean", "mean", DEFAULT_SIGNAL_MEAN, MIN_SIGNAL_MEAN, MAX_SIGNAL_MEAN, "");
1264     printf("FREE SIGNAL_MEAN = %f [%f, %f]\n", mean.getVal(), MIN_SIGNAL_MEAN, MAX_SIGNAL_MEAN);
1265   }
1266   if (FIX_SIGNAL_SIGMA) {
1267     RooRealVar sigma("sigma", "sigma", DEFAULT_SIGNAL_SIGMA, "");
1268     printf("FIXED SIGNAL_SIGMA = %f\n", sigma.getVal());
1269   }
1270   else {
1271     RooRealVar sigma("sigma", "sigma", DEFAULT_SIGNAL_SIGMA, MIN_SIGNAL_SIGMA, MAX_SIGNAL_SIGMA, "");
1272     printf("FREE SIGNAL_SIGMA = %f\n", sigma.getVal());
1273   }
1274   if (FIX_SIGNAL_TAIL) {
1275     RooRealVar tail("tail", "tail", DEFAULT_SIGNAL_TAIL, "");
1276     printf("FIXED SIGNAL_TAIL = %f\n", tail.getVal());
1277   }
1278   else {
1279     RooRealVar tail("tail", "tail", DEFAULT_SIGNAL_TAIL, MIN_SIGNAL_TAIL, MAX_SIGNAL_TAIL, "");
1280     printf("FREE SIGNAL_TAIL = %f\n", tail.getVal());
1281   }
1282   RooRealVar gaussianfrac("gaussianfrac", "gaussianfrac", 1., 0., 1., "");
1283   RooRealVar sigalpha("sigalpha", "sigalpha", 0., -10., 0.);
1284
1285   /* shapes */
1286   if (GAUSSIAN_SIGNAL || forceGaussianSignal) {
1287     printf("USING GAUSSIAN_SIGNAL SHAPE\n");
1288     RooGaussian signal("signal", "signal", x, mean, sigma);
1289   }
1290   else if (GAUSSIANTAIL_SIGNAL) {
1291     printf("USING GAUSSIANTAIL_SIGNAL SHAPE\n");
1292     RooGaussianTail signal("signal", "signal", x, mean, sigma, tail);
1293   }
1294   else if (GAUSSIANTAIL2_SIGNAL) {
1295     printf("USING GAUSSIANTAIL2_SIGNAL SHAPE\n");
1296     RooGaussianTail signal("signal", "signal", x, mean, sigma, sigma);
1297   }
1298   else if (GAUSSIANPLUSGAUSSIANTAIL_SIGNAL) {
1299     printf("USING GAUSSIANPLUSGAUSSIANTAIL_SIGNAL SHAPE\n");
1300     RooGaussian gaussian("gaussian", "gaussian", x, mean, sigma);
1301     RooGaussianTail gaussiantail("gaussiantail", "gaussiantail", x, mean, sigma, tail);
1302     RooAddPdf signal("signal", "signal", RooArgList(gaussian, gaussiantail), gaussianfrac, kTRUE);
1303
1304   }
1305   else if (GAUSSIANPLUSEXPONENTIAL_SIGNAL) {
1306     printf("USING GAUSSIANPLUSEXPONENTIAL_SIGNAL SHAPE\n");
1307     RooGaussian gaussian("gaussian", "gaussian", x, mean, sigma);
1308     RooExponential sigexpo("sigexpo", "sigexpo", x, sigalpha);
1309     RooRealVar sigcutoff("sigcutoff", "sigcutoff", 0.);
1310     RooRealVar sigpower("sigpower", "sigpower", 0.1);
1311     RooFermiCutoff sigfermi("sigfermi", "sigfermi", x, sigcutoff, sigpower);
1312     RooProdPdf exposignal("exposignal", "exposignal", sigfermi, sigexpo, -100.);
1313     RooRealVar gaussianfrac("gaussianfrac", "gaussianfrac", 0.9, 0.7, 1., "");
1314     RooAddPdf signal("signal", "signal", RooArgList(gaussian, exposignal), gaussianfrac, kTRUE);
1315   }
1316   else if (EXPECTED_SIGNAL_TEMPLATE) {
1317     printf("SHAPE OF SIGNAL FROM TEMPLATE HISTOGRAM\n");
1318     RooHistPdf signal("signal", "signal", x, hsig);
1319   }
1320   else if (EXPECTED_SIGNAL_FIT) {
1321     /* fitting bkg1 */
1322     TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1323     hSigExp->Fit(fGaus, "0q");
1324     Double_t sig_mean = fGaus->GetParameter(1);
1325     Double_t sig_sigma = fGaus->GetParameter(2);
1326     mean.setConstant(kFALSE);
1327     mean.setVal(sig_mean);
1328     mean.setRange(sig_mean - 10., sig_mean + 10.);
1329     mean.setConstant(kFALSE);
1330     sigma.setConstant(kFALSE);
1331     sigma.setVal(sig_sigma);
1332     sigma.setRange(sig_sigma * 0.5, sig_sigma * 1.5);
1333     tail.setConstant(kFALSE);
1334     tail.setVal(1.);
1335     tail.setRange(0.5, 5.);
1336     RooGaussianTail signal("signal", "signal", x, mean, sigma, tail);
1337     signal.fitTo(hsig, Range(sig_mean - 5. * sig_sigma, sig_mean + 10. * sig_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1338 #if DISPLAY
1339     TCanvas *cSig_fit = new TCanvas("cSig_fit");
1340     RooPlot *sig_frame = x.frame();
1341     hsig.plotOn(sig_frame);
1342     signal.plotOn(sig_frame, LineColor(kRed));
1343     sig_frame->Draw();
1344     cSig_fit->Update();
1345 #endif
1346     printf("SIGNAL PARAMETERS AFTER FIT OF EXPECTED SIGNAL\n");
1347     printf("mean  = %f +- %f\n", mean.getVal(), mean.getError());
1348     printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError());
1349     printf("tail  = %f +- %f\n", tail.getVal(), tail.getError());
1350     /* scale parameters if requested */
1351     if (SCALE_SIGNAL_SIGMA != 1.) {
1352       printf("SCALE FITTED SIGNAL SIGMA BY %f\n", SCALE_SIGNAL_SIGMA);
1353       sigma.setVal(sigma.getVal() * SCALE_SIGNAL_SIGMA);
1354     }
1355     if (SCALE_SIGNAL_TAIL != 1.) {
1356       printf("SCALE FITTED SIGNAL TAIL BY %f\n", SCALE_SIGNAL_TAIL);
1357       tail.setVal(tail.getVal() * SCALE_SIGNAL_TAIL);
1358     }
1359     /* fix/release parameters if requested */
1360     if (FIX_SIGNAL_MEAN) {
1361       printf("SETTING FITTED SIGNAL MEAN AS CONSTANTS\n");
1362       mean.setConstant(kTRUE);
1363     }
1364     else {
1365       printf("SETTING FITTED SIGNAL MEAN AS FREE\n");
1366       //      mean.setRange(mean.getVal() - 0.25 * TMath::Abs(mean.getVal()), mean.getVal() + 0.25 * TMath::Abs(mean.getVal()));
1367       //      mean.setRange(mean.getVal() - 0.5 * TMath::Abs(mean.getVal()), mean.getVal() + 0.5 * TMath::Abs(mean.getVal()));
1368       mean.setRange(mean.getVal() - 0.2, mean.getVal() + 0.2);
1369     }
1370     if (FIX_SIGNAL_SIGMA) {
1371       printf("SETTING FITTED SIGNAL SIGMA AS CONSTANTS\n");
1372       sigma.setConstant(kTRUE);
1373     }
1374     else {
1375       printf("SETTING FITTED SIGNAL SIGMA AS FREE\n");
1376       sigma.setRange(sigma.getVal() * 0.8, sigma.getVal() * 1.2);
1377       //      sigma.setRange(sigma.getVal() - 0.5, sigma.getVal() + 0.5);
1378     }
1379     if (FIX_SIGNAL_TAIL) {
1380       printf("SETTING FITTED SIGNAL TAIL AS CONSTANTS\n");
1381       tail.setConstant(kTRUE);
1382     }
1383     else {
1384       printf("SETTING FITTED SIGNAL TAIL AS FREE\n");
1385       //      tail.setRange(0.5, 5.0);
1386       tail.setRange(tail.getVal() * 0.8, tail.getVal() * 1.2);
1387       //      tail.setRange(tail.getVal() - 0.5, tail.getVal() + 0.5);
1388     }
1389   }
1390   else {
1391     printf("SHAPE OF SIGNAL NOT DEFINE: using GAUSSIAN_SIGNAL\n");
1392     RooGaussian signal("signal", "signal", x, mean, sigma);
1393   }
1394
1395
1396 #if 0
1397   if (constrainSignal) {
1398 #if 0
1399   /* fix expected signal and constrain parameters if requested */
1400   signal.fitTo(hsig);
1401 #if 0
1402   TCanvas *cConstrainSignal = new TCanvas("cConstrainSignal");
1403   RooPlot *xfr = x.frame();
1404   hsig.plotOn(xfr);
1405   signal.plotOn(xfr, LineColor(kRed));
1406   xfr->Draw();
1407   cConstrainSignal->Update();
1408 #endif
1409   printf("SIGNAL PARAMETERS AFTER FIT OF EXPECTED SIGNAL\n");
1410   printf("mean  = %f +- %f\n", mean.getVal(), mean.getError());
1411   printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError());
1412   printf("tail  = %f +- %f\n", tail.getVal(), tail.getError());
1413   if (constrainSignal) {
1414     mean.setConstant(kTRUE);
1415     sigma.setConstant(kTRUE);
1416     tail.setConstant(kTRUE);
1417   printf("SIGNAL PARAMETERS CONSTRAINED AFTER FIT OF EXPECTED SIGNAL\n");
1418   }
1419 #endif
1420   }
1421 #endif
1422
1423   if (constrainSignal) {
1424     mean.setConstant(kTRUE);
1425     sigma.setConstant(kTRUE);
1426     tail.setConstant(kTRUE);
1427     //    mean.setRange(-0.1, 0.1);
1428     //    sigma.setRange(0.95, 1.05);
1429     //    tail.setRange(0.95, 1.25);
1430   }
1431
1432   printf("***********************************\n");
1433
1434   /*** DEFINE IDENTIFIED BACKGROUND SHAPES ***/
1435
1436   printf("***** IDENTIFIED BACKGROUND SHAPE DEFINITION *****\n");
1437
1438   /* shapes */
1439 if (EXPECTED_BACKGROUND_TEMPLATE) {
1440   printf("USING EXPECTED BACKGROUND TEMPLATE SHAPES FROM HISTOGRAMS\n");
1441     RooHistPdf bkg1("bkg1", "bkg1", x, hbkg1, 0);
1442     RooHistPdf bkg2("bkg2", "bkg2", x, hbkg2, 0);
1443       RooHistPdf bkg3("bkg3", "bkg3", x, hbkg3, 0);
1444       RooHistPdf bkg4("bkg4", "bkg4", x, hbkg4, 0);
1445   }
1446  else if (EXPECTED_BACKGROUND_FIT) {
1447     printf("USING EXPECTED BACKGROUND FITTED SHAPES FROM HISTOGRAMS\n");
1448     /* fitting bkg1 */
1449     TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1450     hBkgExp1->Fit(fGaus, "0q");
1451     Double_t bkgexp1_mean = fGaus->GetParameter(1);
1452     Double_t bkgexp1_sigma = fGaus->GetParameter(2);
1453     RooRealVar mean_bkg1("mean_bkg1", "mean_bkg1", bkgexp1_mean, bkgexp1_mean - 10., bkgexp1_mean + 10., "");
1454     RooRealVar sigma_bkg1("sigma_bkg1", "sigma_bkg1", bkgexp1_sigma, bkgexp1_sigma * 0.5, bkgexp1_sigma * 1.5, "");
1455     RooRealVar tail_bkg1("tail_bkg1", "tail_bkg1", 1.0, 0.5, 5., "");
1456     RooGaussianTail bkg1("bkg1", "bkg1", x, mean_bkg1, sigma_bkg1, tail_bkg1);
1457     bkg1.fitTo(hbkg1, Range(bkgexp1_mean - 5. * bkgexp1_sigma, bkgexp1_mean + 10. * bkgexp1_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1458 #if DISPLAY
1459     TCanvas *cBkg1_fit = new TCanvas("cBkg1_fit");
1460     RooPlot *bkg1_frame = x.frame();
1461     hbkg1.plotOn(bkg1_frame);
1462     bkg1.plotOn(bkg1_frame, LineColor(kCyan+1));
1463     bkg1_frame->Draw();
1464     cBkg1_fit->Update();
1465 #endif
1466     printf("BACKGROUND-1 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-1\n");
1467     printf("mean_bkg1  = %f +- %f\n", mean_bkg1.getVal(), mean_bkg1.getError());
1468     printf("sigma_bkg1 = %f +- %f\n", sigma_bkg1.getVal(), sigma_bkg1.getError());
1469     printf("tail_bkg1  = %f +- %f\n", tail_bkg1.getVal(), tail_bkg1.getError());
1470     /* scale parameters if requested */
1471     if (SCALE_BACKGROUND_SIGMA != 1.) {
1472       printf("SCALE FITTED BACKGROUND-1 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA);
1473       sigma_bkg1.setVal(sigma_bkg1.getVal() * SCALE_BACKGROUND_SIGMA);
1474     }
1475     if (SCALE_BACKGROUND_TAIL != 1.) {
1476       printf("SCALE FITTED BACKGROUND-1 TAIL BY %f\n", SCALE_BACKGROUND_TAIL);
1477       tail_bkg1.setVal(tail_bkg1.getVal() * SCALE_BACKGROUND_TAIL);
1478     }
1479     /* fix/release parameters if requested */
1480     if (FIX_BACKGROUND_MEAN) {
1481       printf("SETTING BACKGROUND-1 FITTED MEAN AS CONSTANTS\n");
1482       mean_bkg1.setConstant(kTRUE);
1483     }
1484     else {
1485       printf("SETTING BACKGROUND-1 FITTED MEAN AS FREE\n");
1486       mean_bkg1.setRange(mean_bkg1.getVal() - 0.25 * TMath::Abs(mean_bkg1.getVal()), mean_bkg1.getVal() + 0.25 * TMath::Abs(mean_bkg1.getVal()));
1487     }
1488     if (FIX_BACKGROUND_SIGMA) {
1489       printf("SETTING BACKGROUND-1 FITTED SIGMA AS CONSTANTS\n");
1490       sigma_bkg1.setConstant(kTRUE);
1491     }
1492     else {
1493       printf("SETTING BACKGROUND-1 FITTED SIGMA AS FREE\n");
1494       sigma_bkg1.setRange(sigma_bkg1.getVal() * 0.75, sigma_bkg1.getVal() * 1.25);
1495     }
1496     if (FIX_BACKGROUND_TAIL) {
1497       printf("SETTING BACKGROUND-1 FITTED TAIL AS CONSTANTS\n");
1498       tail_bkg1.setConstant(kTRUE);
1499     }
1500     else {
1501       printf("SETTING BACKGROUND-1 FITTED TAIL AS FREE\n");
1502       tail_bkg1.setRange(tail_bkg1.getVal() * 0.75, tail_bkg1.getVal() * 1.25);
1503     }
1504     /* fitting bkg2 */
1505     TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1506     hBkgExp2->Fit(fGaus, "0q");
1507     Double_t bkgexp2_mean = fGaus->GetParameter(1);
1508     Double_t bkgexp2_sigma = fGaus->GetParameter(2);
1509     RooRealVar mean_bkg2("mean_bkg2", "mean_bkg2", bkgexp2_mean, bkgexp2_mean - 10., bkgexp2_mean + 10., "");
1510     RooRealVar sigma_bkg2("sigma_bkg2", "sigma_bkg2", bkgexp2_sigma, bkgexp2_sigma * 0.5, bkgexp2_sigma * 1.5, "");
1511     RooRealVar tail_bkg2("tail_bkg2", "tail_bkg2", 1.0, 0.5, 5., "");
1512     RooGaussianTail bkg2("bkg2", "bkg2", x, mean_bkg2, sigma_bkg2, tail_bkg2);
1513     bkg2.fitTo(hbkg2, Range(bkgexp2_mean - 5. * bkgexp2_sigma, bkgexp2_mean + 10. * bkgexp2_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1514 #if DISPLAY
1515     TCanvas *cBkg2_fit = new TCanvas("cBkg2_fit");
1516     RooPlot *bkg2_frame = x.frame();
1517     hbkg2.plotOn(bkg2_frame);
1518     bkg2.plotOn(bkg2_frame,  LineColor(kGreen+1));
1519     bkg2_frame->Draw();
1520     cBkg2_fit->Update();
1521 #endif
1522     printf("BACKGROUND-2 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-2\n");
1523     printf("mean_bkg2  = %f +- %f\n", mean_bkg2.getVal(), mean_bkg2.getError());
1524     printf("sigma_bkg2 = %f +- %f\n", sigma_bkg2.getVal(), sigma_bkg2.getError());
1525     printf("tail_bkg2  = %f +- %f\n", tail_bkg2.getVal(), tail_bkg2.getError());
1526     /* scale parameters if requested */
1527     if (SCALE_BACKGROUND_SIGMA != 1.) {
1528       printf("SCALE FITTED BACKGROUND-2 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA);
1529       sigma_bkg2.setVal(sigma_bkg2.getVal() * SCALE_BACKGROUND_SIGMA);
1530     }
1531     if (SCALE_BACKGROUND_TAIL != 1.) {
1532       printf("SCALE FITTED BACKGROUND-2 TAIL BY %f\n", SCALE_BACKGROUND_TAIL);
1533       tail_bkg2.setVal(tail_bkg2.getVal() * SCALE_BACKGROUND_TAIL);
1534     }
1535     /* fix/release parameters if requested */
1536     if (FIX_BACKGROUND_MEAN) {
1537       printf("SETTING BACKGROUND-2 FITTED MEAN AS CONSTANTS\n");
1538       mean_bkg2.setConstant(kTRUE);
1539     }
1540     else {
1541       printf("SETTING BACKGROUND-2 FITTED MEAN AS FREE\n");
1542       mean_bkg2.setRange(mean_bkg2.getVal() - 0.25 * TMath::Abs(mean_bkg2.getVal()), mean_bkg2.getVal() + 0.25 * TMath::Abs(mean_bkg2.getVal()));
1543     }
1544     if (FIX_BACKGROUND_SIGMA) {
1545       printf("SETTING BACKGROUND-2 FITTED SIGMA AS CONSTANTS\n");
1546       sigma_bkg2.setConstant(kTRUE);
1547     }
1548     else {
1549       printf("SETTING BACKGROUND-2 FITTED SIGMA AS FREE\n");
1550       sigma_bkg2.setRange(sigma_bkg2.getVal() * 0.75, sigma_bkg2.getVal() * 1.25);
1551     }
1552     if (FIX_BACKGROUND_TAIL) {
1553       printf("SETTING BACKGROUND-2 FITTED TAIL AS CONSTANTS\n");
1554       tail_bkg2.setConstant(kTRUE);
1555     }
1556     else {
1557       printf("SETTING BACKGROUND-2 FITTED TAIL AS FREE\n");
1558       tail_bkg2.setRange(tail_bkg2.getVal() * 0.75, tail_bkg2.getVal() * 1.25);
1559     }
1560     /* fitting bkg3 */
1561     TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1562     hBkgExp3->Fit(fGaus, "0q");
1563     Double_t bkgexp3_mean = fGaus->GetParameter(1);
1564     Double_t bkgexp3_sigma = fGaus->GetParameter(2);
1565     RooRealVar mean_bkg3("mean_bkg3", "mean_bkg3", bkgexp3_mean, bkgexp3_mean - 10., bkgexp3_mean + 10., "");
1566     RooRealVar sigma_bkg3("sigma_bkg3", "sigma_bkg3", bkgexp3_sigma, bkgexp3_sigma * 0.5, bkgexp3_sigma * 1.5, "");
1567     RooRealVar tail_bkg3("tail_bkg3", "tail_bkg3", 1., 0.5, 5., "");
1568     RooGaussianTail bkg3("bkg3", "bkg3", x, mean_bkg3, sigma_bkg3, tail_bkg3);
1569     bkg3.fitTo(hbkg3, Range(bkgexp3_mean - 5. * bkgexp3_sigma, bkgexp3_mean + 10. * bkgexp3_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1570 #if DISPLAY
1571     TCanvas *cBkg3_fit = new TCanvas("cBkg3_fit");
1572     RooPlot *bkg3_frame = x.frame();
1573     hbkg3.plotOn(bkg3_frame);
1574     bkg3.plotOn(bkg3_frame,  LineColor(kYellow+1));
1575     bkg3_frame->Draw();
1576     cBkg3_fit->Update();
1577 #endif
1578     printf("BACKGROUND-3 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-3\n");
1579     printf("mean_bkg3  = %f +- %f\n", mean_bkg3.getVal(), mean_bkg3.getError());
1580     printf("sigma_bkg3 = %f +- %f\n", sigma_bkg3.getVal(), sigma_bkg3.getError());
1581     printf("tail_bkg3  = %f +- %f\n", tail_bkg3.getVal(), tail_bkg3.getError());
1582     /* scale parameters if requested */
1583     if (SCALE_BACKGROUND_SIGMA != 1.) {
1584       printf("SCALE FITTED BACKGROUND-3 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA);
1585       sigma_bkg3.setVal(sigma_bkg3.getVal() * SCALE_BACKGROUND_SIGMA);
1586     }
1587     if (SCALE_BACKGROUND_TAIL != 1.) {
1588       printf("SCALE FITTED BACKGROUND-3 TAIL BY %f\n", SCALE_BACKGROUND_TAIL);
1589       tail_bkg3.setVal(tail_bkg3.getVal() * SCALE_BACKGROUND_TAIL);
1590     }
1591     /* fix/release parameters if requested */
1592     if (FIX_BACKGROUND_MEAN) {
1593       printf("SETTING BACKGROUND-3 FITTED MEAN AS CONSTANTS\n");
1594       mean_bkg3.setConstant(kTRUE);
1595     }
1596     else {
1597       printf("SETTING BACKGROUND-3 FITTED MEAN AS FREE\n");
1598       mean_bkg3.setRange(mean_bkg3.getVal() - 0.25 * TMath::Abs(mean_bkg3.getVal()), mean_bkg3.getVal() + 0.25 * TMath::Abs(mean_bkg3.getVal()));
1599     }
1600     if (FIX_BACKGROUND_SIGMA) {
1601       printf("SETTING BACKGROUND-3 FITTED SIGMA AS CONSTANTS\n");
1602       sigma_bkg3.setConstant(kTRUE);
1603     }
1604     else {
1605       printf("SETTING BACKGROUND-3 FITTED SIGMA AS FREE\n");
1606       sigma_bkg3.setRange(sigma_bkg3.getVal() * 0.75, sigma_bkg3.getVal() * 1.25);
1607     }
1608     if (FIX_BACKGROUND_TAIL) {
1609       printf("SETTING BACKGROUND-3 FITTED TAIL AS CONSTANTS\n");
1610       tail_bkg3.setConstant(kTRUE);
1611     }
1612     else {
1613       printf("SETTING BACKGROUND-3 FITTED TAIL AS FREE\n");
1614       tail_bkg3.setRange(tail_bkg3.getVal() * 0.75, tail_bkg3.getVal() * 1.25);
1615     }
1616     /* fitting bkg4 */
1617     TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1618     hBkgExp4->Fit(fGaus, "0q");
1619     Double_t bkgexp4_mean = fGaus->GetParameter(1);
1620     Double_t bkgexp4_sigma = fGaus->GetParameter(2);
1621     RooRealVar mean_bkg4("mean_bkg4", "mean_bkg4", bkgexp4_mean, bkgexp4_mean - 10., bkgexp4_mean + 10., "");
1622     RooRealVar sigma_bkg4("sigma_bkg4", "sigma_bkg4", bkgexp4_sigma, bkgexp4_sigma * 0.5, bkgexp4_sigma * 1.5, "");
1623     RooRealVar tail_bkg4("tail_bkg4", "tail_bkg4", 1., 0.5, 5., "");
1624     RooGaussianTail bkg4("bkg4", "bkg4", x, mean_bkg4, sigma_bkg4, tail_bkg4);
1625     bkg4.fitTo(hbkg4, Range(bkgexp4_mean - 5. * bkgexp4_sigma, bkgexp4_mean + 10. * bkgexp4_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1626 #if DISPLAY
1627     TCanvas *cBkg4_fit = new TCanvas("cBkg4_fit");
1628     RooPlot *bkg4_frame = x.frame();
1629     hbkg4.plotOn(bkg4_frame);
1630     bkg4.plotOn(bkg4_frame,  LineColor(kYellow+2));
1631     bkg4_frame->Draw();
1632     cBkg4_fit->Update();
1633 #endif
1634     printf("BACKGROUND-4 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-4\n");
1635     printf("mean_bkg4  = %f +- %f\n", mean_bkg4.getVal(), mean_bkg4.getError());
1636     printf("sigma_bkg4 = %f +- %f\n", sigma_bkg4.getVal(), sigma_bkg4.getError());
1637     printf("tail_bkg4  = %f +- %f\n", tail_bkg4.getVal(), tail_bkg4.getError());
1638     /* scale parameters if requested */
1639     if (SCALE_BACKGROUND_SIGMA != 1.) {
1640       printf("SCALE FITTED BACKGROUND-4 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA);
1641       sigma_bkg4.setVal(sigma_bkg4.getVal() * SCALE_BACKGROUND_SIGMA);
1642     }
1643     if (SCALE_BACKGROUND_TAIL != 1.) {
1644       printf("SCALE FITTED BACKGROUND-4 TAIL BY %f\n", SCALE_BACKGROUND_TAIL);
1645       tail_bkg4.setVal(tail_bkg4.getVal() * SCALE_BACKGROUND_TAIL);
1646     }
1647     /* fix/release parameters if requested */
1648     if (FIX_BACKGROUND_MEAN) {
1649       printf("SETTING BACKGROUND-4 FITTED MEAN AS CONSTANTS\n");
1650       mean_bkg4.setConstant(kTRUE);
1651     }
1652     else {
1653       printf("SETTING BACKGROUND-4 FITTED MEAN AS FREE\n");
1654       mean_bkg4.setRange(mean_bkg4.getVal() - 0.25 * TMath::Abs(mean_bkg4.getVal()), mean_bkg4.getVal() + 0.25 * TMath::Abs(mean_bkg4.getVal()));
1655     }
1656     if (FIX_BACKGROUND_SIGMA) {
1657       printf("SETTING BACKGROUND-4 FITTED SIGMA AS CONSTANTS\n");
1658       sigma_bkg4.setConstant(kTRUE);
1659     }
1660     else {
1661       printf("SETTING BACKGROUND-4 FITTED SIGMA AS FREE\n");
1662       sigma_bkg4.setRange(sigma_bkg4.getVal() * 0.75, sigma_bkg4.getVal() * 1.25);
1663     }
1664     if (FIX_BACKGROUND_TAIL) {
1665       printf("SETTING BACKGROUND-4 FITTED TAIL AS CONSTANTS\n");
1666       tail_bkg4.setConstant(kTRUE);
1667     }
1668     else {
1669       printf("SETTING BACKGROUND-4 FITTED TAIL AS FREE\n");
1670       tail_bkg4.setRange(tail_bkg4.getVal() * 0.75, tail_bkg4.getVal() * 1.25);
1671     }
1672   }
1673   else if (GAUSSIAN_BACKGROUND) {
1674     printf("USING GAUSSIAN BACKGROUND SHAPES (not reccomended)\n");
1675     RooRealVar mean1("mean1", "mean1", hBkgExp1->GetMean(), hBkgExp1->GetMean() * 0.95, hBkgExp1->GetMean() * 1.05, "");
1676     RooRealVar sigma1("sigma1", "sigma1", hBkgExp1->GetRMS(), hBkgExp1->GetRMS() * 0.5, hBkgExp1->GetRMS() * 1.5, "");
1677     RooGaussian bkg1("bkg1", "bkg1", x, mean1, sigma1);
1678
1679     RooRealVar mean2("mean2", "mean2", hBkgExp2->GetMean(), hBkgExp2->GetMean() * 0.95, hBkgExp2->GetMean() * 1.05, "");
1680     RooRealVar sigma2("sigma2", "sigma2", hBkgExp2->GetRMS(), hBkgExp2->GetRMS() * 0.5, hBkgExp2->GetRMS() * 1.5, "");
1681     RooGaussian bkg2("bkg2", "bkg2", x, mean2, sigma2);
1682
1683     RooRealVar mean3("mean3", "mean3", hBkgExp3->GetMean(), hBkgExp3->GetMean() * 0.95, hBkgExp3->GetMean() * 1.05, "");
1684     RooRealVar sigma3("sigma3", "sigma3", hBkgExp3->GetRMS(), hBkgExp3->GetRMS() * 0.5, hBkgExp3->GetRMS() * 1.5, "");
1685     RooGaussian bkg3("bkg3", "bkg3", x, mean3, sigma3);
1686   }
1687   else {
1688     printf("SHAPE OF BACKGROUND NOT DEFINE: using EXPECTED_BACKGROUND\n");
1689     RooHistPdf bkg1("bkg1", "bkg1", x, hbkg1, 0);
1690     RooHistPdf bkg2("bkg2", "bkg2", x, hbkg2, 0);
1691     RooHistPdf bkg3("bkg3", "bkg3", x, hbkg3, 0);
1692     RooHistPdf bkg4("bkg4", "bkg4", x, hbkg4, 0);
1693   }
1694   printf("**************************************************\n");
1695
1696   /*** DEFINE MISMATCH BACKGROUND SHAPE ***/
1697
1698   printf("***** MISMATCH BACKGROUND SHAPE DEFINITION *****\n");
1699
1700   /* variables and generic shapes */
1701   Double_t expectedCutoff = hBkgExp3->GetMean();
1702   Double_t expectedCutoffRMS = hBkgExp3->GetRMS();
1703   //    RooRealVar cutoff("cutoff", "cutoff", -30., -50., 0., "");
1704   RooRealVar cutoff("cutoff", "cutoff", expectedCutoff, expectedCutoff - 3. * expectedCutoffRMS, expectedCutoff, "");
1705   //  RooRealVar cutoff("cutoff", "cutoff", expectedCutoff, "");
1706   //  RooRealVar power("power", "power", 1., 0.5, 5.0, "");
1707   RooRealVar power("power", "power", 1., "");
1708   RooFermiCutoff fermi("fermi", "fermi", x, cutoff, power);
1709   RooRealVar alpha("alpha", "alpha", 0., -10., 0., "");
1710   RooExponential expo("expo", "expo", x, alpha);
1711   RooUniform uniform("uniform", "uniform", x);
1712
1713   /* mismatch shape */
1714   if (EXPECTED_MISMATCH) {
1715     printf("USING EXPECTED MISMATCH SHAPE FROM HISTOGRAMS\n");
1716     RooHistPdf mismatch("mismatch", "mismatch", x, hmismatch, 0);
1717   }
1718   else if (UNIFORM_MISMATCH) {
1719     printf("USING UNIFORM BACKGROUND SHAPE WITH CUTOFF\n");
1720     RooProdPdf mismatch("mismatch", "mismatch", fermi, uniform, -100.);
1721   }
1722   else if (EXPONENTIAL_MISMATCH) {
1723     printf("USING EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n");
1724     RooProdPdf mismatch("mismatch", "mismatch", fermi, expo, -100.);
1725   }
1726   else if (DOUBLEEXPONENTIAL_MISMATCH) {
1727     printf("USING DOUBLE-EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n");
1728     RooRealVar alpha1("alpha1", "alpha1", 0., -10., 0., "");
1729     RooRealVar alpha2("alpha2", "alpha2", 0., -10., 0., "");
1730     RooRealVar frac("frac", "frac", 0.5, 0., 1., "");
1731     RooGenericPdf doubleexpo("doubleexpo", "TMath::Exp(alpha1 * x) + frac * TMath::Exp(alpha2 * x)", RooArgSet(x, alpha1, alpha2, frac));
1732     RooProdPdf mismatch("mismatch", "mismatch", fermi, doubleexpo, -100.);
1733   }
1734   else if (UNIFORMPLUSEXPONENTIAL_MISMATCH) {
1735     printf("USING UNIFORM-PLUS-EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n");
1736     RooRealVar funiform("funiform", "funiform", 100., 0., 100000., "");
1737     RooRealVar fexpo("fexpo", "fexpo", 100., 0., 100000., "");
1738     RooAddPdf uniformplusexpo("uniformplusexpo", "uniformplusexpo", RooArgList(uniform, expo), RooArgList(funiform, fexpo), kFALSE);
1739     RooProdPdf mismatch("mismatch", "mismatch", fermi, uniformplusexpo, -100.);
1740   }
1741   else {
1742     printf("SHAPE OF MISMATCH NOT DEFINE: using EXPECTED_MISMATCH\n");
1743     RooHistPdf mismatch("mismatch", "mismatch", x, hmismatch, 0);
1744   }
1745   printf("************************************************\n");
1746
1747   /*** DEFINE THE MODEL ***/
1748
1749   printf("***** MODEL DEFINITION *****\n");
1750
1751   /* variables */
1752   Double_t integral = hdata.sumEntries();
1753   RooRealVar nsignal("nsignal", "nsignal", 0.3 * integral, 0., integral);
1754   RooRealVar nbkg1("nbkg1", "nbkg1", 0.3 * integral, 0., integral);
1755   RooRealVar nbkg2("nbkg2", "nbkg2", 0.3 * integral, 0., integral);
1756   RooRealVar nbkg3("nbkg3", "nbkg3", 0.3 * integral, 0., integral);
1757   RooRealVar nbkg4("nbkg4", "nbkg4", 0.3 * integral, 0., integral);
1758   RooRealVar nmismatch("nmismatch", "nmismatch", 0.1 * integral, 0., integral);
1759
1760 if (!fitBkg1) {
1761   nbkg1.setVal(0.);
1762   nbkg1.setConstant(kTRUE);
1763  }
1764 if (!fitBkg2) {
1765   nbkg2.setVal(0.);
1766   nbkg2.setConstant(kTRUE);
1767  }
1768 if (!fitBkg3) {
1769   nbkg3.setVal(0.);
1770   nbkg3.setConstant(kTRUE);
1771  }
1772 if (!fitBkg4) {
1773   nbkg4.setVal(0.);
1774   nbkg4.setConstant(kTRUE);
1775  }
1776
1777   RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3, bkg4, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nbkg3, nbkg4, nmismatch));
1778
1779 #if 0
1780   /* the model */
1781   if (USE_ELECTRON_BACKGROUND && fitElectrons && fitPions) {
1782     printf("USING ELECTRON BACKGROUND\n");
1783     if (NO_MISMATCH) {
1784       printf("NOT USING MISMATCH BACKGROUND\n");
1785       nmismatch.setVal(0.);
1786       RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3/*, bkg4*/), RooArgList(nsignal, nbkg1, nbkg2, nbkg3/*, nbkg4*/));
1787     }
1788     else {
1789       printf("USING MISMATCH BACKGROUND\n");
1790       RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3/*, bkg4*/, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nbkg3/*, nbkg4*/, nmismatch));
1791     }
1792   }
1793   else if (!USE_ELECTRON_BACKGROUND || !fitElectrons) {
1794     printf("NOT USING ELECTRON BACKGROUND\n");
1795     nbkg3.setVal(0.);
1796     nbkg4.setVal(0.);
1797     if (NO_MISMATCH) {
1798       printf("NOT USING MISMATCH BACKGROUND\n");
1799       nmismatch.setVal(0.);
1800       RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2), RooArgList(nsignal, nbkg1, nbkg2));
1801     }
1802     else {
1803       printf("USING MISMATCH BACKGROUND\n");
1804       RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nmismatch));
1805     }
1806   }
1807   printf("****************************\n");
1808 #endif
1809
1810
1811
1812
1813   /*** FIT ***/
1814
1815   printf("***** FIT *****\n");
1816
1817   printf("SIGNAL PARAMETERS BEFORE GLOBAL FIT\n");
1818   printf("mean  = %f +- %f\n", mean.getVal(), mean.getError());
1819   printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError());
1820   printf("tail  = %f +- %f\n", tail.getVal(), tail.getError());
1821
1822   /* fit and draw */
1823   if (canvas) canvas->cd();
1824   //  model.chi2FitTo(hdata, Extended(kTRUE), Verbose(kFALSE), SumW2Error(kFALSE), Range(-40., 140.), Binned(kTRUE));
1825 //  model.fitTo(hdata, Extended(kTRUE), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10), Range(-10., 10.));
1826 model.fitTo(hdata, Range(rangeMin, rangeMax), Extended(kTRUE), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1827
1828   printf("***************\n");
1829
1830   /*** DRAW ***/
1831 #if DISPLAY
1832   RooPlot *xframe = x.frame();
1833 hdata.plotOn(xframe, XErrorSize(0), DrawOption("PZ"));
1834  printf("plotting model...\n");
1835 model.plotOn(xframe, LineWidth(2)/*, Precision(1.e-4)*/);
1836  printf("plotting signal...\n");
1837 model.plotOn(xframe, Components(signal), LineWidth(2), LineColor(kRed)/*, Precision(1.e-4)*/);
1838  printf("plotting bkg1...\n");
1839   model.plotOn(xframe, Components(bkg1), LineWidth(2), LineStyle(kDashed), LineColor(kCyan+1));
1840  printf("plotting bkg2...\n");
1841   model.plotOn(xframe, Components(bkg2), LineWidth(2), LineStyle(kDashed), LineColor(kGreen+1));
1842   if (USE_ELECTRON_BACKGROUND) {
1843     model.plotOn(xframe, Components(bkg3), LineWidth(2), LineStyle(kDashed), LineColor(kYellow+1));
1844     model.plotOn(xframe, Components(bkg4), LineWidth(2), LineStyle(kDashed), LineColor(kYellow+2));
1845   }
1846   if (!NO_MISMATCH)
1847     model.plotOn(xframe, Components(mismatch), LineWidth(2), LineStyle(kDashed), LineColor(kGray+1));
1848   hSignal->SetFillColor(kYellow);
1849   hSignal->SetLineWidth(0);
1850   hSignal->SetFillStyle(0);
1851   hSignal->SetMinimum(0.1);
1852   hSignal->GetXaxis()->SetRangeUser(rangeMin, rangeMax);
1853 //  hSignal->Draw();
1854 xframe->SetMinimum(0.1);
1855   xframe->Draw();
1856 #endif
1857   if (canvas) canvas->Update();
1858
1859   /*** COMPUTE CHI2 ***/
1860   Double_t datax, datapoint, datapoint_err, modelpoint;
1861   Double_t chi2 = 0.;
1862   Int_t ndf = 0;
1863   for (Int_t ibin = 0; ibin < hSignal->GetNbinsX(); ibin++) {
1864     datax = hSignal->GetBinCenter(ibin + 1);
1865     if (datax < rangeMin || datax > rangeMax) continue;
1866     datapoint = hSignal->GetBinContent(ibin + 1);
1867     datapoint_err = hSignal->GetBinError(ibin + 1);
1868     if (datapoint_err == 0.) continue;
1869     x.setVal(datax);
1870     modelpoint = model.getVal();
1871     chi2 += (datapoint - modelpoint) * (datapoint - modelpoint) / (datapoint_err * datapoint_err);
1872     ndf++;
1873   }
1874
1875
1876 /*** PRINT FIT OUTPUT ***/
1877 printf("***** CHI-SQUARE *****\n");
1878   printf("chi-square = %f\n", chi2);
1879   printf("NDF        = %d\n", ndf);
1880   printf("chi2/NDF   = %f\n", ndf > 0 ? chi2 / ndf : 0.);
1881   printf("***** SIGNAL-SHAPE INFO *****\n");
1882   printf("mean      = %f +- %f\n", mean.getVal(), mean.getError());
1883   printf("sigma     = %f +- %f\n", sigma.getVal(), sigma.getError());
1884   printf("tail      = %f +- %f\n", tail.getVal(), tail.getError());
1885   printf("pure gaus = %f +- %f\n", gaussianfrac.getVal(), gaussianfrac.getError());
1886 printf("*****************************\n");
1887 printf("***** COUNTS *****\n");
1888 printf("total     = %f\n", hSignal->Integral(-1, -1));
1889   printf("integral  = %f\n", hdata.sumEntries());
1890   printf("nsignal   = %f +- %f\n", nsignal.getVal(), nsignal.getError());
1891   printf("nbkg1     = %f +- %f\n", nbkg1.getVal(), nbkg1.getError());
1892   printf("nbkg2     = %f +- %f\n", nbkg2.getVal(), nbkg2.getError());
1893   printf("nbkg3     = %f +- %f\n", nbkg3.getVal(), nbkg3.getError());
1894   printf("nbkg4     = %f +- %f\n", nbkg4.getVal(), nbkg4.getError());
1895   printf("nmismatch = %f +- %f\n", nmismatch.getVal(), nmismatch.getError());
1896   printf("******************\n");
1897
1898   /*** OUTPUT FIT PARAMS ***/
1899
1900   param[kMean] = mean.getVal();
1901   param_err[kMean] = mean.getError();
1902   param[kSigma] = sigma.getVal();
1903   param_err[kSigma] = sigma.getError();
1904   param[kTail] = tail.getVal();
1905   param_err[kTail] = tail.getError();
1906   param[kTotalCounts] = hSignal->GetEntries();
1907   param_err[kTotalCounts] = sqrt(hSignal->GetEntries());
1908   param[kIntegralCounts] = hdata.sumEntries();
1909   param_err[kIntegralCounts] = sqrt(hdata.sumEntries());
1910   param[kSignalCounts] = nsignal.getVal();
1911   param_err[kSignalCounts] = nsignal.getError();
1912   param[kBkg1Counts] = nbkg1.getVal();
1913   param_err[kBkg1Counts] = nbkg1.getError();
1914   param[kBkg2Counts] = nbkg2.getVal();
1915   param_err[kBkg2Counts] = nbkg2.getError();
1916   param[kBkg3Counts] = nbkg3.getVal();
1917   param_err[kBkg3Counts] = nbkg3.getError();
1918   param[kBkg4Counts] = nbkg4.getVal();
1919   param_err[kBkg4Counts] = nbkg4.getError();
1920   param[kMismatchCounts] = nmismatch.getVal();
1921   param_err[kMismatchCounts] = nmismatch.getError();
1922
1923   return;
1924 }
1925
1926 //___________________________________________________________________________________
1927
1928 TOFpid_efficiency(Int_t ipart, Int_t icharge, Int_t icent, Int_t useTPCcut = -1, Double_t minsignalFrac = 0., Int_t nTrials = 1000)
1929 {
1930
1931   /* prepare centrality name */
1932   Char_t centName[1024];
1933   if (icent < 0 || icent >= NcentralityBins)
1934     sprintf(centName, "cent0090");
1935   else
1936     sprintf(centName, "cent%02d%02d", centralityBin[icent], centralityBin[icent + 1]);
1937
1938   /* fraction names */
1939   const Char_t *fractionName[AliPID::kSPECIES][AliPID::kSPECIES] = {
1940     "hSignalFraction", "hBkg4Fraction", "hBkg3Fraction", "hBkg1Fraction", "hBkg2Fraction",
1941     "hBkg4Fraction", "hSignalFraction", "hBkg3Fraction", "hBkg1Fraction", "hBkg2Fraction",
1942     "hBkg3Fraction", "hBkg4Fraction", "hSignalFraction", "hBkg1Fraction", "hBkg2Fraction",
1943     "hBkg3Fraction", "hBkg4Fraction", "hBkg1Fraction", "hSignalFraction", "hBkg2Fraction",
1944     "hBkg3Fraction", "hBkg4Fraction", "hBkg1Fraction", "hBkg2Fraction", "hSignalFraction"
1945   };
1946
1947   enum ETPCcut_t {
1948     kCurrentDir,
1949     k11Cut,
1950     k12Cut,
1951     k21Cut,
1952     k22Cut,
1953     k23Cut,
1954     k32Cut,
1955     k33Cut,
1956     kNTPCcuts
1957   };
1958   const Char_t *tpccutdir[8] = {
1959     ".",
1960     "TOFpid_cutOnTPC[-1.0,1.0]",
1961     "TOFpid_cutOnTPC[1.0,2.0]",
1962     "TOFpid_cutOnTPC[-2.0,-1.0]",
1963     "TOFpid_cutOnTPC[-2.0,2.0]",
1964     "TOFpid_cutOnTPC[2.0,3.0]",
1965     "TOFpid_cutOnTPC[-3.0,-2.0]",
1966     "TOFpid_cutOnTPC[-3.0,3.0]"
1967   };
1968
1969   /* get data */
1970   Char_t filename[1024];
1971   TH1D *hAccepted[AliPID::kSPECIES][kNTPCcuts];
1972   TH1D *hIdentified[AliPID::kSPECIES][kNTPCcuts];
1973   TH1D *hEfficiencyIn[AliPID::kSPECIES][kNTPCcuts];
1974   TH1D *hFraction[AliPID::kSPECIES][AliPID::kSPECIES][kNTPCcuts];
1975   for (Int_t itpccut = 0; itpccut < kNTPCcuts; itpccut++) {
1976
1977     /* check whether we use this cut */
1978     if (useTPCcut >= 0 && itpccut != useTPCcut)
1979       continue;
1980
1981     for (Int_t iipart = 0; iipart < AliPID::kSPECIES; iipart++) {
1982       /* skip electrons and muons */
1983       if (iipart == AliPID::kMuon)
1984         continue;
1985
1986       /* open file */
1987       sprintf(filename, "%s/TOFspectrum_%s_%s_%s_%sID.root", tpccutdir[itpccut], centName, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart));
1988       TFile *filein = TFile::Open(filename);
1989       if (!filein || !filein->IsOpen()) {
1990         printf("cannot open %s\n", filename);
1991         return;
1992       }
1993
1994       /* get accepted tracks */
1995       hAccepted[iipart][itpccut] = (TH1D *)filein->Get("PostAnalysis/hAcceptedTracks");
1996       if (!hAccepted[iipart][itpccut]) {
1997         printf("cannot find PostAnalysis/hAcceptedTracks in %s\n", filename);
1998         return;
1999       }
2000
2001       /* get identified tracks */
2002       hIdentified[iipart][itpccut] = (TH1D *)filein->Get("PostAnalysis/hIdentifiedCounts");
2003       if (!hIdentified[iipart][itpccut]) {
2004         printf("cannot find PostAnalysis/hIdentifiedCounts in %s\n", filename);
2005         return;
2006       }
2007
2008       /* compute efficiency */
2009       hEfficiencyIn[iipart][itpccut] = new TH1D(*hIdentified[iipart][itpccut]);
2010       hEfficiencyIn[iipart][itpccut]->Divide(hEfficiencyIn[iipart][itpccut], hAccepted[iipart][itpccut], 1., 1., "B");
2011
2012       /* get fractions */
2013       for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
2014         /* skip muons */
2015         if (iiipart == AliPID::kMuon)
2016           continue;
2017
2018         hFraction[iipart][iiipart][itpccut] = (TH1D *)filein->Get(Form("PostAnalysis/%s", fractionName[iipart][iiipart]));
2019         if (!hFraction[iipart][iiipart][itpccut]) {
2020           printf("cannot find PostAnalysis/%s in %s\n", fractionName[iipart][iiipart], filename);
2021           return;
2022         }
2023       }
2024     }
2025   }
2026
2027   /* prepare output efficiency histos */
2028   TH1D *hEfficiencyOut[AliPID::kSPECIES];
2029   for (Int_t iipart = 0; iipart < AliPID::kSPECIES; iipart++) {
2030     hEfficiencyOut[iipart] = new TH1D(Form("hPIDEff_%d_%s_%s_%sID", icent, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NptBins, ptBin);
2031   }
2032
2033
2034   /*** MATRIX ***/
2035
2036   const Int_t nPart = 4;
2037   Int_t partIndex[4] = {AliPID::kElectron, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
2038
2039   TMatrixD Meffin(nPart, 1);
2040   TMatrixD Meffout(nPart, 1);
2041   TMatrixD Mfrac(nPart, nPart);
2042
2043   Double_t eff[4], effe[4], frac[4][4], frace[4][4], geneff, genfrac;
2044   Bool_t gotbestcut[4];
2045
2046   TH1F *hEffTemp[4];
2047   for (Int_t iipart = 0; iipart < nPart; iipart++)
2048     hEffTemp[iipart] = new TH1F(Form("hEffTemp_%d", iipart), "", 20000, -10., 10.);
2049
2050   /* loop over pt-bins */
2051   for (Int_t ibin = 0; ibin < NptBins; ibin++) {
2052
2053     /* reset temp histos */
2054     for (Int_t iipart = 0; iipart < nPart; iipart++)
2055       hEffTemp[iipart]->Reset();
2056
2057     /* get measured data */
2058     for (Int_t iipart = 0; iipart < nPart; iipart++) {
2059
2060       /* select the best set of cuts */
2061       Double_t signalFrac, maxsignalFrac = minsignalFrac;
2062       Int_t bestcut = -1;
2063       gotbestcut[iipart] = kFALSE;
2064       for (Int_t itpccut = 0; itpccut < kNTPCcuts; itpccut++) {
2065
2066         /* check whether we use this cut */
2067         if (useTPCcut >= 0 && itpccut != useTPCcut)
2068           continue;
2069
2070         /* maximize the signal fraction */
2071         signalFrac = hFraction[partIndex[iipart]][partIndex[iipart]][itpccut]->GetBinContent(ibin + 1);
2072         if (signalFrac > maxsignalFrac) {
2073           maxsignalFrac = signalFrac;
2074           bestcut = itpccut;
2075           gotbestcut[iipart] = kTRUE;
2076         }
2077       }
2078       if (bestcut < 0)
2079         continue;
2080
2081       eff[iipart] = hEfficiencyIn[partIndex[iipart]][bestcut]->GetBinContent(ibin + 1);
2082       effe[iipart] = hEfficiencyIn[partIndex[iipart]][bestcut]->GetBinError(ibin + 1);
2083       for (Int_t iiipart = 0; iiipart < nPart; iiipart++) {
2084         frac[iipart][iiipart] = hFraction[partIndex[iipart]][partIndex[iiipart]][bestcut]->GetBinContent(ibin + 1);
2085         frace[iipart][iiipart] = hFraction[partIndex[iipart]][partIndex[iiipart]][bestcut]->GetBinError(ibin + 1);
2086       }
2087     }
2088
2089     /* check best cuts */
2090     Bool_t skip = kFALSE;
2091     for (Int_t iipart = 0; iipart < nPart; iipart++)
2092       if (!gotbestcut[iipart])
2093         skip = kTRUE;
2094     if (skip) continue;
2095
2096     /* loop over trials */
2097     for (Int_t itry = 0; itry < nTrials; itry++) {
2098
2099       /* setup matrix */
2100       for (Int_t iipart = 0; iipart < nPart; iipart++) {
2101         geneff = gRandom->Gaus(eff[iipart], effe[iipart]);
2102         if (geneff < 0.) geneff = 0.;
2103         if (geneff > 1.) geneff = 1.;
2104         Meffin[iipart] = geneff != 0. ? 1. / geneff : 0.;
2105         for (Int_t iiipart = 0; iiipart < nPart; iiipart++) {
2106           genfrac = gRandom->Gaus(frac[iipart][iiipart], frace[iipart][iiipart]);
2107           if (genfrac < 0.) genfrac = 0.;
2108           if (genfrac > 1.) genfrac = 1.;
2109           Mfrac[iipart][iiipart] = genfrac;
2110         }
2111       }
2112
2113       /* invert matrix */
2114       TDecompLU lu(Mfrac);
2115       TMatrixD Minv(nPart, nPart);
2116       if (!lu.Decompose()) continue;
2117       lu.Invert(Minv);
2118       Meffout = Minv * Meffin;
2119
2120       /* fill histos */
2121       for (Int_t iipart = 0; iipart < nPart; iipart++) {
2122         if (Meffout.GetMatrixArray()[iipart] > 0.)
2123           hEffTemp[iipart]->Fill(1. / Meffout.GetMatrixArray()[iipart]);
2124       }
2125
2126     } /* end of loop over trials */
2127
2128
2129     /* get average and RMS */
2130     for (Int_t iipart = 0; iipart < nPart; iipart++) {
2131       hEfficiencyOut[partIndex[iipart]]->SetBinContent(ibin + 1, hEffTemp[iipart]->GetMean());
2132       hEfficiencyOut[partIndex[iipart]]->SetBinError(ibin + 1, hEffTemp[iipart]->GetRMS());
2133     }
2134
2135   } /* end of loop over pt-bins */
2136
2137
2138   hEfficiencyOut[AliPID::kPion]->SetMarkerStyle(20);
2139   hEfficiencyOut[AliPID::kKaon]->SetMarkerStyle(20);
2140   hEfficiencyOut[AliPID::kProton]->SetMarkerStyle(20);
2141   hEfficiencyOut[AliPID::kElectron]->SetMarkerStyle(20);
2142
2143   hEfficiencyOut[AliPID::kPion]->SetMarkerColor(4);
2144   hEfficiencyOut[AliPID::kKaon]->SetMarkerColor(8);
2145   hEfficiencyOut[AliPID::kProton]->SetMarkerColor(2);
2146   hEfficiencyOut[AliPID::kElectron]->SetMarkerColor(1);
2147
2148   hEfficiencyOut[AliPID::kPion]->Draw();
2149   hEfficiencyOut[AliPID::kKaon]->Draw("same");
2150   hEfficiencyOut[AliPID::kProton]->Draw("same");
2151   hEfficiencyOut[AliPID::kElectron]->Draw("same");
2152
2153
2154   /* output */
2155   TFile *fileout = TFile::Open(Form("TOFpid_efficiency_cent%d_%s_%s.root", icent, AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
2156   hEfficiencyOut[AliPID::kPion]->Write();
2157   hEfficiencyOut[AliPID::kKaon]->Write();
2158   hEfficiencyOut[AliPID::kProton]->Write();
2159   hEfficiencyOut[AliPID::kElectron]->Write();
2160   fileout->Close();
2161
2162 }
2163
2164 TOFpid_normalize(TH1D *h, Double_t nevents = 8.42446600000000000e+06)
2165 {
2166
2167   Double_t counts, counts_err;
2168   for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++) {
2169     counts = h->GetBinContent(ibin + 1);
2170     counts_err = h->GetBinError(ibin + 1);
2171     counts /= h->GetBinWidth(ibin + 1);
2172     counts_err /= h->GetBinWidth(ibin + 1);
2173     h->SetBinContent(ibin + 1, counts);
2174     h->SetBinError(ibin + 1, counts_err);
2175   }
2176
2177   h->Scale(1. / nevents);
2178
2179 }
2180
2181 TOFpid_normalizeAndwrite(const Char_t *fileoutname, const Char_t *corrstring = "")
2182 {
2183
2184   TFile *fpiplus = TFile::Open("TOFpid_spectrum_pion_plus.root");
2185   TFile *fpiminus = TFile::Open("TOFpid_spectrum_pion_minus.root");
2186   TFile *fkaplus = TFile::Open("TOFpid_spectrum_kaon_plus.root");
2187   TFile *fkaminus = TFile::Open("TOFpid_spectrum_kaon_minus.root");
2188   TFile *fprplus = TFile::Open("TOFpid_spectrum_proton_plus.root");
2189   TFile *fprminus = TFile::Open("TOFpid_spectrum_proton_minus.root");
2190
2191   hpiplus = (TH1F *)fpiplus->Get(Form("hSignal%sCounts", corrstring));
2192   hpiminus = (TH1F *)fpiminus->Get(Form("hSignal%sCounts", corrstring));
2193   hkaplus = (TH1F *)fkaplus->Get(Form("hSignal%sCounts", corrstring));
2194   hkaminus = (TH1F *)fkaminus->Get(Form("hSignal%sCounts", corrstring));
2195   hprplus = (TH1F *)fprplus->Get(Form("hSignal%sCounts", corrstring));
2196   hprminus = (TH1F *)fprminus->Get(Form("hSignal%sCounts", corrstring));
2197
2198   hpiplus->SetName("hpiplus");
2199   hpiminus->SetName("hpiminus");
2200   hkaplus->SetName("hkaplus");
2201   hkaminus->SetName("hkaminus");
2202   hprplus->SetName("hprplus");
2203   hprminus->SetName("hprminus");
2204
2205   TOFpid_normalize(hpiplus);
2206   TOFpid_normalize(hpiminus);
2207   TOFpid_normalize(hkaplus);
2208   TOFpid_normalize(hkaminus);
2209   TOFpid_normalize(hprplus);
2210   TOFpid_normalize(hprminus);
2211
2212   TFile *fileout = TFile::Open(fileoutname, "RECREATE");
2213   hpiplus->Write("hpiplus");
2214   hpiminus->Write("hpiminus");
2215   hkaplus->Write("hkaplus");
2216   hkaminus->Write("hkaminus");
2217   hprplus->Write("hprplus");
2218   hprminus->Write("hprminus");
2219   fileout->Close();
2220
2221 }
2222
2223 /**************************************************************/
2224
2225 TOFpid_rawSpectra(const Char_t *destdir = ".")
2226 {
2227
2228   Int_t marker[2] = {20, 21};
2229   Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2230   Char_t *partLatex[AliPID::kSPECIES][2] = {
2231     "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2232   };
2233
2234   Float_t deltaRapidity = rapidityMaxCut - rapidityMinCut;
2235
2236   TFile *fileout = TFile::Open("TOF_rawSpectra.root", "RECREATE");
2237   TH1D *hRaw;
2238   Char_t title[1024];
2239   for (Int_t icent = -1; icent < NcentralityBins; icent++)
2240     for (Int_t icharge = 0; icharge < kNCharges; icharge++)
2241       for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
2242         hRaw = TOFpid_rawSpectra(destdir, ipart, icharge, icent);
2243         if (!hRaw) continue;
2244         hRaw->Scale(1. / deltaRapidity);
2245         hRaw->SetMarkerStyle(marker[icharge]);
2246         hRaw->SetMarkerColor(color[ipart]);
2247         hRaw->SetLineColor(1);
2248         hRaw->SetLineWidth(1);
2249         if (icent == -1)
2250           sprintf(title, "%s (MB);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge]);
2251         else
2252           sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1]);
2253         hRaw->SetTitle(title);
2254         if (icent == -1)
2255           hRaw->SetName(Form("hRaw_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
2256         else
2257           hRaw->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2258         fileout->cd();
2259         hRaw->Write();
2260       }
2261
2262   fileout->Close();
2263 }
2264
2265 TH1D *
2266 TOFpid_rawSpectra(const Char_t *dirname, Int_t ipart, Int_t icharge, Int_t icent)
2267 {
2268
2269   /* open data */
2270   Char_t outfilename[1024];
2271   if (icent < 0 || icent >= NcentralityBins) {
2272     sprintf(outfilename, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2273   }
2274   else {
2275     sprintf(outfilename, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", dirname, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2276   }
2277   TFile *filein = TFile::Open(outfilename);
2278   if (!filein || !filein->IsOpen()) {
2279     printf("cannot open %s\n", outfilename);
2280     return NULL;
2281   }
2282   /* get data */
2283   TH1D *h = (TH1D *)filein->Get("RawSpectra/hNormalizedRawYield");
2284   if (!h) {
2285     printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename);
2286     return NULL;
2287   }
2288
2289   return h;
2290
2291 }
2292
2293 TOFpid_rawSpectra_mismatchCorrected(const Char_t *destdir = ".")
2294 {
2295
2296   Int_t marker[2] = {20, 21};
2297   Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2298   Char_t *partLatex[AliPID::kSPECIES][2] = {
2299     "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2300   };
2301
2302   Float_t deltaRapidity = rapidityMaxCut - rapidityMinCut;
2303
2304   TFile *fileout = TFile::Open("TOF_rawSpectra_mismatchCorrected.root", "RECREATE");
2305   TH1D *hRaw;
2306   Char_t title[1024];
2307   for (Int_t icent = -1; icent < NcentralityBins; icent++)
2308     for (Int_t icharge = 0; icharge < kNCharges; icharge++)
2309       for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
2310         hRaw = TOFpid_rawSpectra_mismatchCorrected(destdir, ipart, icharge, icent);
2311         if (!hRaw) continue;
2312         hRaw->Scale(1. / deltaRapidity);
2313         hRaw->SetMarkerStyle(marker[icharge]);
2314         hRaw->SetMarkerColor(color[ipart]);
2315         hRaw->SetLineColor(1);
2316         hRaw->SetLineWidth(1);
2317         if (icent == -1)
2318           sprintf(title, "%s (MB);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge]);
2319         else
2320           sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1]);
2321         hRaw->SetTitle(title);
2322         if (icent == -1)
2323           hRaw->SetName(Form("hRaw_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
2324         else
2325           hRaw->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2326         fileout->cd();
2327         hRaw->Write();
2328       }
2329
2330   fileout->Close();
2331 }
2332
2333 TH1D *
2334 TOFpid_rawSpectra_mismatchCorrected(const Char_t *dirname, Int_t ipart, Int_t icharge, Int_t icent)
2335 {
2336
2337   /* open data */
2338   Char_t outfilename[1024];
2339   if (icent < 0 || icent >= NcentralityBins) {
2340     sprintf(outfilename, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2341   }
2342   else {
2343     sprintf(outfilename, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", dirname, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2344   }
2345   TFile *filein = TFile::Open(outfilename);
2346   if (!filein || !filein->IsOpen()) {
2347     printf("cannot open %s\n", outfilename);
2348     return NULL;
2349   }
2350   /* get data */
2351   TH1D *h = (TH1D *)filein->Get("RawSpectra/hNormalizedMismatchCorrectedRawYield");
2352   if (!h) {
2353     printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename);
2354     return NULL;
2355   }
2356
2357   return h;
2358
2359 }
2360
2361 /**************************************************************/
2362
2363 TOFpid_electronCorrection()
2364 {
2365
2366     Int_t marker[2] = {20, 21};
2367   Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2368   Char_t *partLatex[AliPID::kSPECIES][2] = {
2369     "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2370   };
2371
2372
2373   TFile *fileout = TFile::Open("TOF_electronCorrection.root", "RECREATE");
2374   TH1D *hCorr;
2375   TProfile *pCorrAv = new TProfile("pCorrAv", "", NptBins, ptBin, "s");
2376   Char_t title[1024];
2377   for (Int_t icent = 0; icent < NcentralityBins; icent++)
2378     for (Int_t icharge = 0; icharge < kNCharges; icharge++)
2379       for (Int_t ipart = 2; ipart < 3; ipart++) {
2380         hCorr = TOFpid_electronCorrection(ipart, icharge, icent);
2381         hCorr->SetMarkerStyle(marker[icharge]);
2382         hCorr->SetMarkerColor(color[ipart]);
2383         hCorr->SetLineColor(1);
2384         hCorr->SetLineWidth(1);
2385         sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);electron-subtracted pion fraction;", partLatex[ipart][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1]);
2386         hCorr->SetTitle(title);
2387         hCorr->SetName(Form("hElectronCorr_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2388         fileout->cd();
2389         hCorr->Write();
2390         /* fill for average correction */
2391         for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2392           pCorrAv->Fill(hCorr->GetBinCenter(ipt + 1), hCorr->GetBinContent(ipt + 1));
2393         }
2394       }
2395   hCorr = new TH1D("hElectronCorr_average", ";p_{T} (GeV/c);electron-subtracted pion fraction;", NptBins, ptBin);
2396   for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2397     hCorr->SetBinContent(ipt + 1, pCorrAv->GetBinContent(ipt + 1));
2398     hCorr->SetBinError(ipt + 1, pCorrAv->GetBinError(ipt + 1));
2399   }
2400   hCorr->Write();
2401   fileout->Close();
2402 }
2403
2404 TH1D *
2405 TOFpid_electronCorrection(Int_t ipart, Int_t icharge, Int_t icent)
2406 {
2407
2408   Int_t marker[2] = {20, 21};
2409   Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2410   Char_t *partLatex[AliPID::kSPECIES][2] = {
2411     "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2412   };
2413
2414   TH1D *hr = TOFpid_systematics_ratio("defaultFit_electronCut", "defaultFit", ipart, icharge, icent, "electronCorrection", "", marker[icharge], color[ipart], kTRUE);
2415   TF1 fOne("fOne", "1.", fitPtMin[ipart], fitPtMax[ipart]);
2416   hr->Add(&fOne, -1.);
2417   hr->Scale(1. / 0.866);
2418   hr->Add(&fOne, 1.);
2419
2420   return hr;
2421 }
2422
2423 /**************************************************************/
2424
2425 SystematicCheck(const Char_t *defaultname = "defaultFit")
2426 {
2427
2428   const Int_t ndata = 6;
2429   const Char_t *name[ndata] = {
2430     "signalFit_fixed_scaleSigma_09",
2431     "signalFit_fixed_scaleSigma_11",
2432     "signalFit_fixed_scaleTail_09",
2433     "signalFit_fixed_scaleTail_11",
2434     "signalFit_fixed_scaleSigma_09_scaleTail_11",
2435     "signalFit_fixed_scaleSigma_11_scaleTail_09"
2436   };
2437   for (Int_t idata = 0; idata < ndata; idata++)
2438     SystematicCheck(name[idata], defaultname);
2439
2440 }
2441
2442 SystematicCheck(const Char_t *checkname, const Char_t *defaultname = "defaultFit")
2443 {
2444
2445   gROOT->LoadMacro("HistoUtils.C");
2446
2447   Char_t filename1[1024];
2448   Char_t filename2[1024];
2449
2450   Int_t marker[2] = {20, 25};
2451   Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2452   Char_t *partLatex[AliPID::kSPECIES][2] = {
2453     "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2454   };
2455
2456
2457   TFile *fileout = TFile::Open(Form("SystematicCheck_%s.root", checkname), "RECREATE");
2458   Char_t title[1024];
2459   TH1 *hd;
2460   for (Int_t icent = 0; icent < NcentralityBins; icent++)
2461     for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
2462       for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
2463
2464         if (icent < 0 || icent >= NcentralityBins) {
2465           sprintf(filename1, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", checkname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2466           sprintf(filename2, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", defaultname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2467         }
2468         else {
2469           sprintf(filename1, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", checkname, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2470           sprintf(filename2, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", defaultname, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2471         }
2472
2473         hd = HistoUtils_systematics(filename1, filename2, "RawSpectra/hNormalizedRawYield", NULL);
2474         hd->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2475         sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge], centralityBin[icent], centralityBin[icent + 1]);
2476         hd->SetTitle(title);
2477         hd->SetLineColor(1);
2478         hd->SetLineWidth(1);
2479         hd->SetMarkerStyle(marker[icharge]);
2480         hd->SetMarkerColor(color[ipart]);
2481         fileout->cd();
2482         hd->Write();
2483         delete hd;
2484       }
2485   fileout->Close();
2486
2487 }
2488
2489 /**************************************************************/
2490
2491 TOFpid_systematics()
2492 {
2493
2494   TFile *fileout = TFile::Open("TOFpid_systematics.root", "RECREATE");
2495   TH1D *hSys;
2496   for (Int_t icent = 0; icent < NcentralityBins; icent++)
2497     for (Int_t icharge = 0; icharge < kNCharges; icharge++)
2498       for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
2499         hSys = TOFpid_systematics(ipart, icharge, icent);
2500         fileout->cd();
2501         hSys->Write(Form("hRawSys_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2502         delete hSys;
2503       }
2504
2505   fileout->Close();
2506 }
2507
2508 TH1D *
2509 TOFpid_systematics(Int_t ipart, Int_t icharge, Int_t icent)
2510 {
2511
2512   TH1D *hSignalFit = TOFpid_systematics_signalFit(ipart, icharge, icent);
2513   TH1D *hBkgFit = TOFpid_systematics_bkgFit(ipart, icharge, icent);
2514   TH1D *hFitRange = TOFpid_systematics_fitRange(ipart, icharge, icent);
2515
2516   TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin);
2517   Double_t sigsys, bkgsys, rangesys, totsys = 0.;
2518   for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2519     sigsys = hSignalFit->GetBinError(ipt + 1);
2520     bkgsys = hBkgFit->GetBinError(ipt + 1);
2521     rangesys = hFitRange->GetBinError(ipt + 1);
2522     totsys = TMath::Sqrt(sigsys * sigsys + bkgsys * bkgsys + rangesys * rangesys);
2523     hSys->SetBinContent(ipt + 1, totsys);
2524   }
2525
2526   hSys->Draw();
2527
2528   delete hSignalFit;
2529   delete hBkgFit;
2530   delete hFitRange;
2531
2532   return hSys;
2533 }
2534
2535 TH1D *
2536 TOFpid_systematics_fitRange(Int_t ipart, Int_t icharge, Int_t icent)
2537 {
2538
2539   TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin);
2540   hArea->SetMinimum(0.8);
2541   hArea->SetMaximum(1.2);
2542   hArea->Draw();
2543
2544   TH1D *hr = TOFpid_systematics_ratio("default_wideRange", "default", ipart, icharge, icent, "wide", "wide-range fit;p_{T} (GeV/c);ratio wrt. default;", 20, 4);
2545   hr->Draw("same");
2546
2547   TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin);
2548   hSys->SetFillStyle(0);
2549   hSys->SetMarkerSize(0);
2550   for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2551     Double_t val, max = 0.;
2552     if (hr->GetBinContent(ipt + 1) == 0.) continue;
2553     val = TMath::Abs(hr->GetBinContent(ipt + 1) - 1.);
2554     hSys->SetBinContent(ipt + 1, 1.);
2555     hSys->SetBinError(ipt + 1, val);
2556   }
2557   hSys->Draw("same, E2");
2558
2559   delete hr;
2560
2561   return hSys;
2562 }
2563
2564 TH1D *
2565 TOFpid_systematics_signalFit(Int_t ipart, Int_t icharge, Int_t icent)
2566 {
2567
2568   const Int_t ndata = 6;
2569   const Char_t *name[ndata] = {
2570     "signalFit_fixed_scaleSigma_09",
2571     "signalFit_fixed_scaleSigma_11",
2572     "signalFit_fixed_scaleTail_09",
2573     "signalFit_fixed_scaleTail_11",
2574     "signalFit_fixed_scaleSigma_09_scaleTail_11",
2575     "signalFit_fixed_scaleSigma_11_scaleTail_09"
2576   };
2577   const Char_t *title[ndata] = {
2578     "-10% #sigma;p_{T} (GeV/c); ratio",
2579     "+10% #sigma;p_{T} (GeV/c); ratio",
2580     "-10% #tau;p_{T} (GeV/c); ratio",
2581     "+10% #tau;p_{T} (GeV/c); ratio",
2582     "-10% #sigma, +10% #tau;p_{T} (GeV/c); ratio",
2583     "+10% #sigma, -10% #tau;p_{T} (GeV/c); ratio"
2584   };
2585   Int_t marker[ndata] = {22, 28, 22, 28, 22, 28};
2586   Int_t color[ndata] = {2, 2, 8, 8, 4, 4};
2587
2588   TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin);
2589   hArea->SetMinimum(0.8);
2590   hArea->SetMaximum(1.2);
2591   hArea->Draw();
2592
2593   TH1D *hr[ndata];
2594   for (Int_t idata = 0; idata < ndata; idata++) {
2595     hr[idata] = TOFpid_systematics_ratio(name[idata], "defaultFit", ipart, icharge, icent, name[idata], title[idata], marker[idata], color[idata]);
2596     hr[idata]->Draw("same");
2597   }
2598
2599   TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin);
2600   hSys->SetFillStyle(0);
2601   hSys->SetMarkerSize(0);
2602   for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2603     Double_t val, max = 0.;
2604     for (Int_t idata = 0; idata < ndata; idata++) {
2605       if (hr[idata]->GetBinContent(ipt + 1) == 0.) continue;
2606       val = TMath::Abs(hr[idata]->GetBinContent(ipt + 1) - 1.);
2607       if (val > 0.9) continue;
2608       if (val > max)
2609         max = val;
2610     }
2611     hSys->SetBinContent(ipt + 1, 1.);
2612     hSys->SetBinError(ipt + 1, max);
2613   }
2614   hSys->Draw("same, E2");
2615
2616   //  delete hArea;
2617   //  for (Int_t idata = 0; idata < ndata; idata++)
2618   //    delete hr[idata];
2619
2620   return hSys;
2621 }
2622
2623 TH1D *
2624 TOFpid_systematics_bkgFit(Int_t ipart, Int_t icharge, Int_t icent)
2625 {
2626
2627   const Int_t ndata = 6;
2628   const Char_t *name[ndata] = {
2629     "bkgFit_fixed_scaleSigma_09",
2630     "bkgFit_fixed_scaleSigma_11",
2631     "bkgFit_fixed_scaleTail_09",
2632     "bkgFit_fixed_scaleTail_11",
2633     "bkgFit_fixed_scaleSigma_09_scaleTail_11",
2634     "bkgFit_fixed_scaleSigma_11_scaleTail_09"
2635   };
2636   const Char_t *title[ndata] = {
2637     "-10% #sigma;p_{T} (GeV/c); ratio",
2638     "+10% #sigma;p_{T} (GeV/c); ratio",
2639     "-10% #tau;p_{T} (GeV/c); ratio",
2640     "+10% #tau;p_{T} (GeV/c); ratio",
2641     "-10% #sigma, +10% #tau;p_{T} (GeV/c); ratio",
2642     "+10% #sigma, -10% #tau;p_{T} (GeV/c); ratio"
2643   };
2644   Int_t marker[ndata] = {22, 28, 22, 28, 22, 28};
2645   Int_t color[ndata] = {2, 2, 8, 8, 4, 4};
2646
2647   TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin);
2648   hArea->SetMinimum(0.5);
2649   hArea->SetMaximum(1.5);
2650   hArea->Draw();
2651
2652   TH1D *hr[ndata];
2653   for (Int_t idata = 0; idata < ndata; idata++) {
2654     hr[idata] = TOFpid_systematics_ratio(name[idata], "defaultFit", ipart, icharge, icent, name[idata], title[idata], marker[idata], color[idata]);
2655     hr[idata]->Draw("same");
2656   }
2657
2658   TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin);
2659   hSys->SetFillStyle(0);
2660   hSys->SetMarkerSize(0);
2661   for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2662     Double_t val, max = 0.;
2663     for (Int_t idata = 0; idata < ndata; idata++) {
2664       if (hr[idata]->GetBinContent(ipt + 1) == 0.) continue;
2665       val = TMath::Abs(hr[idata]->GetBinContent(ipt + 1) - 1.);
2666       if (val > 0.9) continue;
2667       if (val > max)
2668         max = val;
2669     }
2670     hSys->SetBinContent(ipt + 1, 1.);
2671     hSys->SetBinError(ipt + 1, max);
2672   }
2673   hSys->Draw("same, E2");
2674
2675   //  delete hArea;
2676   //  for (Int_t idata = 0; idata < ndata; idata++)
2677   //    delete hr[idata];
2678
2679   return hSys;
2680 }
2681
2682 TH1D *
2683 TOFpid_systematics_ratio(const Char_t *dirname1, const Char_t *dirname2, Int_t ipart, Int_t icharge, Int_t icent, const Char_t *name = "rawRatio", const Char_t *title = ";p_{T} (GeV/c);raw yield ratio;", Int_t marker = 20, Int_t color = 2, Bool_t correlated = kFALSE)
2684 {
2685
2686   TH1D *hr = new TH1D("hr", "", NptBins, ptBin);
2687
2688   /* open data */
2689   Char_t outfilename1[1024];
2690   Char_t outfilename2[1024];
2691   if (icent < 0 || icent >= NcentralityBins) {
2692     sprintf(outfilename1, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname1, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2693     sprintf(outfilename2, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname2, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2694   }
2695   else {
2696     sprintf(outfilename1, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", dirname1, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2697     sprintf(outfilename2, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", dirname2, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2698   }
2699   TFile *filein1 = TFile::Open(outfilename1);
2700   if (!filein1 || !filein1->IsOpen()) {
2701     printf("cannot open %s\n", outfilename1);
2702     return;
2703   }
2704   TFile *filein2 = TFile::Open(outfilename2);
2705   if (!filein2 || !filein2->IsOpen()) {
2706     printf("cannot open %s\n", outfilename2);
2707     return;
2708   }
2709   /* get data */
2710   TH1D *h1 = (TH1D *)filein1->Get("RawSpectra/hNormalizedRawYield");
2711   if (!h1) {
2712     printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename1);
2713     return;
2714   }
2715   TH1D *h2 = (TH1D *)filein2->Get("RawSpectra/hNormalizedRawYield");
2716   if (!h2) {
2717     printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename2);
2718     return;
2719   }
2720   /* ratio */
2721   if (correlated) hr->Divide(h1, h2, 1., 1., "B");
2722   else hr->Divide(h1, h2);
2723   hr->SetNameTitle(name, title);
2724   hr->SetMarkerStyle(marker);
2725   hr->SetMarkerColor(color);
2726
2727   filein1->Close();
2728   filein2->Close();
2729
2730   return hr;
2731
2732 }
2733