#include "CommonDefs.C" #define DISPLAY 0 /* SIGNAL SHAPE */ Bool_t GAUSSIAN_SIGNAL = kFALSE; Bool_t GAUSSIANTAIL_SIGNAL = kFALSE; Bool_t GAUSSIANTAIL2_SIGNAL = kFALSE; Bool_t GAUSSIANPLUSGAUSSIANTAIL_SIGNAL = kFALSE; Bool_t GAUSSIANPLUSEXPONENTIAL_SIGNAL = kFALSE; Bool_t EXPECTED_SIGNAL_TAIL = kTRUE; Bool_t EXPECTED_SIGNAL_TEMPLATE = kTRUE; Bool_t EXPECTED_SIGNAL_FIT = kFALSE; /* SIGNAL PARAMETERS */ Bool_t FIX_SIGNAL_MEAN = kFALSE; Bool_t FIX_SIGNAL_SIGMA = kFALSE; Bool_t FIX_SIGNAL_TAIL = kFALSE; Float_t SCALE_SIGNAL_SIGMA = 1.; Float_t SCALE_SIGNAL_TAIL = 1.; /* OTHER STUFF */ Char_t *SIGNAL_PARAM_FILE = NULL;//"signalParamFile.root"; Float_t DEFAULT_SIGNAL_MEAN = 0.; Float_t MIN_SIGNAL_MEAN = -0.2; Float_t MAX_SIGNAL_MEAN = 0.2.; Float_t DEFAULT_SIGNAL_SIGMA = 1.; Float_t MIN_SIGNAL_SIGMA = 0.8; Float_t MAX_SIGNAL_SIGMA = 1.2; Float_t DEFAULT_SIGNAL_TAIL = 1.; Float_t MIN_SIGNAL_TAIL = 0.5; Float_t MAX_SIGNAL_TAIL = 1.5; /* BACKGROUND */ Bool_t EXPECTED_BACKGROUND_TAIL = kTRUE; Bool_t EXPECTED_BACKGROUND_TEMPLATE = kTRUE; Bool_t EXPECTED_BACKGROUND_FIT = kFALSE; Bool_t GAUSSIAN_BACKGROUND = kFALSE; Bool_t USE_ELECTRON_BACKGROUND = kTRUE; /* BACKGROUND PARAMETERS */ Bool_t FIX_BACKGROUND_MEAN = kTRUE; Bool_t FIX_BACKGROUND_SIGMA = kTRUE; Bool_t FIX_BACKGROUND_TAIL = kTRUE; Float_t SCALE_BACKGROUND_SIGMA = 1.; Float_t SCALE_BACKGROUND_TAIL = 1.; /* MISMATCH */ Bool_t NO_MISMATCH = kFALSE; Bool_t EXPECTED_MISMATCH = kTRUE; Bool_t EXPONENTIAL_MISMATCH = kFALSE; Bool_t UNIFORM_MISMATCH = kFALSE; Bool_t DOUBLEEXPONENTIAL_MISMATCH = kFALSE; Bool_t UNIFORMPLUSEXPONENTIAL_MISMATCH = kFALSE; /* AUTOMATIC CONSTRAINS */ Float_t FIT_ELECTRONS_PT_MIN = 0.; Float_t FIT_ELECTRONS_PT_MAX = 1.0; Float_t FIT_MUONS_PT_MIN = 0.7; Float_t FIT_MUONS_PT_MAX = 0.; Float_t FIT_PIONS_PT_MIN = 0.; Float_t FIT_PIONS_PT_MAX = 5.; Float_t CONSTRAINSIGNAL_LIMIT = 0.; Float_t CONSTRAINTAIL_LIMIT = 0.; enum EFitParams_t { kMean, kSigma, kTail, kTotalCounts, kIntegralCounts, kSignalCounts, kBkg1Counts, kBkg2Counts, kBkg3Counts, kBkg4Counts, kMismatchCounts, kNFitParams }; /* fit output params name */ const Char_t *fitParamName[kNFitParams] = { "Mean", "Sigma", "Tail", "TotalCounts", "IntegralCounts", "SignalCounts", "Bkg1Counts", "Bkg2Counts", "Bkg3Counts", "Bkg4Counts", "MismatchCounts" }; /* fit output params title */ const Char_t *fitParamTitle[kNFitParams] = { "Signal mean;p_{T} (GeV/c);#mu (ps)", "Signal sigma;p_{T} (GeV/c);#sigma (ps)", "Signal tail;p_{T} (GeV/c);#sigma_{tail} (ps)", "Total counts;p_{T} (GeV/c);counts", "Total counts within fit range;p_{T} (GeV/c);counts", "Signal counts;p_{T} (GeV/c);counts", "Background-1 counts;p_{T} (GeV/c);counts", "Background-2 counts;p_{T} (GeV/c);counts", "Background-3 counts;p_{T} (GeV/c);counts", "Background-4 counts;p_{T} (GeV/c);counts", "Mismatch counts within fit range;p_{T} (GeV/c);counts" }; /**************************************************************/ /*** GENERATION OF TEMPLATE HISTOS ****************************/ /**************************************************************/ const Int_t NmismatchTrials = 1; const Int_t NexpectedTrials = 1; /**************************************************************/ /*** HISTOS AND BINNING ***************************************/ /**************************************************************/ /**************************************************************/ /**************************************************************/ enum EHistoParam_t { kCentrality, kTOFsigma, kPt, kTPCsigma, kNHistoParams }; /**************************************************************/ /**************************************************************/ const Int_t NtofsigmaBins = 1750; Double_t tofsigmaBin[NtofsigmaBins + 1]; Double_t tofsigmaMin = -100., tofsigmaMax = 250., tofsigmaStep = (tofsigmaMax - tofsigmaMin) / NtofsigmaBins; /**************************************************************/ const Int_t NtofsignalBins = 2000; Double_t tofsignalBin[NtofsignalBins + 1]; Double_t tofsignalMin = -24400., tofsignalMax = 24400., tofsignalStep = (tofsignalMax - tofsignalMin) / NtofsignalBins; /**************************************************************/ const Int_t NtofresoBins = 500; Double_t tofresoBin[NtofsignalBins + 1]; Double_t tofresoMin = 0., tofresoMax = 500., tofresoStep = (tofresoMax - tofresoMin) / NtofresoBins; /**************************************************************/ /**************************************************************/ const Int_t NtpcsigmaBins = 10; Double_t tpcsigmaBin[NtpcsigmaBins + 1]; Double_t tpcsigmaMin = -5., tpcsigmaMax = 5., tpcsigmaStep = (tpcsigmaMax - tpcsigmaMin) / NtpcsigmaBins; /**************************************************************/ Int_t NparamsBins[kNHistoParams] = {NcentralityBins, NtofsigmaBins, NptBins, NtpcsigmaBins}; Double_t *paramBin[kNHistoParams] = {centralityBin, tofsigmaBin, ptBin, tpcsigmaBin}; /**************************************************************/ /**************************************************************/ /**************************************************************/ /**************************************************************/ /**************************************************************/ AliTOFGeometry tofGeo; Float_t c = TMath::C() * 1.e2 / 1.e12; /* cm/ps */ Float_t c_1 = 1. / c; Double_t GenerateRandomHit(TH1F *hT0Fill, Double_t t0fill, Int_t index) { Int_t det[5]; Float_t length, timeexp, pos[3]; /* compute length and expected time */ tofGeo.GetVolumeIndices(index, det); tofGeo.GetPosPar(det, pos); length = 0.; for (Int_t i = 0; i < 3; i++) length += pos[i] * pos[i]; length = TMath::Sqrt(length); timeexp = length * c_1; Double_t hittime = hT0Fill->GetRandom() - t0fill + timeexp; return hittime; } /**************************************************************/ 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) { printf("****************************************\n"); printf("RUNNING TOF PID:\n"); if (dorapidityCut) printf("RAPIDITY-CUT: %s\n", AliPID::ParticleName(ipart)); else printf("NO RAPIDITY-CUT\n"); printf("CHARGE: %s\n", chargeName[icharge]); printf("PARTICLE-ID: %s\n", AliPID::ParticleName(iipart)); if (electronCut) printf("-> ELECTRON CUT REQUESTED\n"); if (cutOnTPC) printf(" -> CUT-ON-TPC [%3.1f,%3.1f]\n", tpcsignalMin, tpcsignalMax); printf("****************************************\n"); /* include path for ACLic */ gSystem->AddIncludePath("-I$ALICE_ROOT/include"); gSystem->AddIncludePath("-I$ALICE_ROOT/TOF"); /* load libraries */ gSystem->Load("libANALYSIS"); gSystem->Load("libANALYSISalice"); /* build analysis task class */ gROOT->LoadMacro("AliAnalysisParticle.cxx+g"); gROOT->LoadMacro("AliAnalysisEvent.cxx+g"); gROOT->LoadMacro("AliAnalysisTrack.cxx+g"); /* create TOF response with tail */ // gROOT->LoadMacro("~/ALICE.2011/ANALYSIS/TOFSpectraPbPb/macros/TOFsignal.C"); gROOT->LoadMacro("TOFsignal.C"); TF1 *fTOFsignal = new TF1("fTOFsignal", TOFsignal, -2440., 2440., 4); fTOFsignal->SetParameter(0, 1.); fTOFsignal->SetParameter(1, 0.); fTOFsignal->SetParameter(2, tofReso); fTOFsignal->SetParameter(3, tofTail); /* open file, get tree and connect */ TFile *filein = TFile::Open(filename); TTree *treein = (TTree *)filein->Get("aodTree"); printf("got \"aodTree\": %d entries\n", treein->GetEntries()); AliAnalysisEvent *analysisEvent = new AliAnalysisEvent(); TClonesArray *analysisTrackArray = new TClonesArray("AliAnalysisTrack"); AliAnalysisTrack *analysisTrack = NULL; treein->SetBranchAddress("AnalysisEvent", &analysisEvent); treein->SetBranchAddress("AnalysisTrack", &analysisTrackArray); /* open hT0fill for mismatch */ TFile *filein_T0Fill = TFile::Open(t0FillOnlineFileName); TH1F *hT0Fill = (TH1F *)filein_T0Fill->Get("hT0Fill"); Double_t t0fill = t0Fill_offset; /* open enabled flag map */ TH1F *hEnabledFlag = NULL; if (enabledChannelsFileName) { TFile *enabledfile = TFile::Open(enabledChannelsFileName); hEnabledFlag = (TH1F *)enabledfile->Get("hEnabledFlag"); } /**************************************************************/ /*** HISTOS ***************************************************/ /**************************************************************/ /* run-time binning */ for (Int_t ibin = 0; ibin < NtofsigmaBins + 1; ibin++) tofsigmaBin[ibin] = tofsigmaMin + ibin * tofsigmaStep; for (Int_t ibin = 0; ibin < NtofsignalBins + 1; ibin++) tofsignalBin[ibin] = tofsignalMin + ibin * tofsignalStep; for (Int_t ibin = 0; ibin < NtpcsigmaBins + 1; ibin++) tpcsigmaBin[ibin] = tpcsigmaMin + ibin * tpcsigmaStep; for (Int_t ibin = 0; ibin < NtofresoBins + 1; ibin++) tofresoBin[ibin] = tofresoMin + ibin * tofresoStep; /* histos */ TH1F *hAcceptedEvents = new TH1F(Form("hAcceptedEvents_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin); TH2I *hAcceptedTracks = new TH2I(Form("hAcceptedTracks_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin); TH3I *hTOFreso = new TH3I(Form("hTOFreso_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofresoBins, tofresoBin); TH3I *hTOFpid = new TH3I(Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin); 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); TH3I *hTOFmismatch = new TH3I(Form("hTOFmismatch_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin); TH3I *hTOFexpected[AliPID::kSPECIES]; for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) { 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); } 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); 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); TH3I *hTOFexpected_delta[AliPID::kSPECIES]; for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) { 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); } /**************************************************************/ /**************************************************************/ /**************************************************************/ /* TOF PID response */ AliTOFPIDResponse tofResponse; tofResponse.SetTimeResolution(tofReso); /* TPC PID response */ AliTPCPIDResponse *tpcResponse = AliAnalysisTrack::GetTPCResponse(); /* start stopwatch */ TStopwatch timer; timer.Start(); /* verbose */ printf("***** RUNNING for %s %s *****\n", chargeName[icharge], AliPID::ParticleName(ipart)); if (!dorapidityCut) printf("***** NO RAPIDITY CUT REQUESTED *****\n"); if (cutOnTPC) { printf("***** CUT-ON-TPC requested %3.1f-%3.1f *****\n", tpcsignalMin, tpcsignalMax); } /* loop over events */ Bool_t hastofpid; Int_t charge, index; UShort_t dedxN; Double_t cent, p, pt, mt, tofsignal, tpcsignal; Double_t dedx, bethe, deltadedx, dedx_sigma, ptpc; 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; /* open TZERO calibfile */ TH1 *hCentrality_TZEROA_mean = NULL; TH1 *hCentrality_TZEROA_sigma = NULL; TH1 *hCentrality_TZEROC_mean = NULL; TH1 *hCentrality_TZEROC_sigma = NULL; TH1 *hCentrality_TOF_mean = NULL; TH1 *hCentrality_TOF_TZEROA_mean = NULL; TH1 *hCentrality_TOF_TZEROC_mean = NULL; TH1 *hCentrality_TOF_TZEROTOF_mean = NULL; TH1 *hCentrality_TZEROA_reso = NULL; TH1 *hCentrality_TZEROC_reso = NULL; TH1 *hCentrality_TZEROAND_reso = NULL; if (1) { TFile *calibfile = TFile::Open("TZEROcalibration.root"); hCentrality_TZEROA_mean = (TH1 *)calibfile->Get("hCentrality_TZEROA_mean"); hCentrality_TZEROA_sigma = (TH1 *)calibfile->Get("hCentrality_TZEROA_sigma"); hCentrality_TZEROC_mean = (TH1 *)calibfile->Get("hCentrality_TZEROC_calib"); hCentrality_TZEROC_sigma = (TH1 *)calibfile->Get("hCentrality_TZEROC_sigma"); hCentrality_TOF_mean = (TH1 *)calibfile->Get("hCentrality_TOF_mean"); hCentrality_TOF_TZEROA_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROA_mean"); hCentrality_TOF_TZEROC_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROC_mean"); hCentrality_TOF_TZEROTOF_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROTOF_mean"); TFile *resofile = TFile::Open("TZEROresolution.root"); hCentrality_TZEROA_reso = (TH1 *)resofile->Get("hTZEROA_reso"); hCentrality_TZEROC_reso = (TH1 *)resofile->Get("hTZEROC_reso"); hCentrality_TZEROAND_reso = (TH1 *)resofile->Get("hTZEROAND_reso"); } Double_t TZEROA_mean; Double_t TZEROA_sigma; Double_t TZEROC_mean; Double_t TZEROC_sigma; Double_t TOF_mean; Double_t TOF_TZEROA_mean; Double_t TOF_TZEROC_mean; Double_t TOF_TZEROTOF_mean; Double_t TZEROA; Double_t TZEROA_reso; Bool_t hasTZEROA; Double_t TZEROC; Double_t TZEROC_reso; Bool_t hasTZEROC; Double_t TZEROAND; Double_t TZEROAND_reso; Bool_t hasTZEROAND; Double_t TZEROTOF; Double_t TZEROTOF_reso; Bool_t hasTZEROTOF; Double_t TZEROMEAN; Double_t TZEROMEAN_weight; Double_t TZEROBEST; Double_t TZEROBEST_reso; Double_t param[kNHistoParams]; for (Int_t iev = startEv; iev < treein->GetEntries() && iev < evMax; iev++) { /* get event */ treein->GetEvent(iev); if (iev % 10000 == 0) printf("iev = %d\n", iev); /* check event */ if (!analysisEvent->AcceptEvent(acceptEventType)) continue; /*** ACCEPTED EVENT ***/ /* apply time-zero TOF correction */ // analysisEvent->ApplyTimeZeroTOFCorrection(); /* get centrality */ cent = analysisEvent->GetCentralityPercentile(centralityEstimator); /* vertex */ Double_t vertexz = analysisEvent->GetVertexZ(); /* fill histos */ hAcceptedEvents->Fill(cent); /* TZERO corrections */ Int_t icent; for (icent = 0; icent < NcentralityBins; icent++) if (cent < centralityBin[icent + 1]) break; TZEROA_mean = hCentrality_TZEROA_mean ? hCentrality_TZEROA_mean->GetBinContent(icent + 1) : 0.; TZEROA_sigma = hCentrality_TZEROA_sigma ? hCentrality_TZEROA_sigma->GetBinContent(icent + 1) : 1000.; TZEROC_mean = hCentrality_TZEROC_mean ? hCentrality_TZEROC_mean->GetBinContent(icent + 1) : 0.; TZEROC_sigma = hCentrality_TZEROC_sigma ? hCentrality_TZEROC_sigma->GetBinContent(icent + 1) : 1000.; TOF_mean = hCentrality_TOF_mean ? hCentrality_TOF_mean->GetBinContent(icent + 1) : 0.; TOF_TZEROA_mean = hCentrality_TOF_TZEROA_mean ? hCentrality_TOF_TZEROA_mean->GetBinContent(icent + 1) : 0.; TOF_TZEROC_mean = hCentrality_TOF_TZEROC_mean ? hCentrality_TOF_TZEROC_mean->GetBinContent(icent + 1) : 0.; TOF_TZEROTOF_mean = hCentrality_TOF_TZEROTOF_mean ? hCentrality_TOF_TZEROTOF_mean->GetBinContent(icent + 1) : 0.; TZEROA_reso = hCentrality_TZEROA_reso ? hCentrality_TZEROA_reso->GetBinContent(icent + 1) : 70.; TZEROC_reso = hCentrality_TZEROC_reso ? hCentrality_TZEROC_reso->GetBinContent(icent + 1) : 70.; TZEROAND_reso = hCentrality_TZEROAND_reso ? hCentrality_TZEROAND_reso->GetBinContent(icent + 1) : 50.; /* simulate inefficient TZERO */ Bool_t resetTZERO = kFALSE; if (forcetimezeroineff > 0.) if (gRandom->Uniform() < forcetimezeroineff) resetTZERO = kTRUE; /* loop over tracks */ for (Int_t itrk = 0; itrk < analysisTrackArray->GetEntries(); itrk++) { /* get track */ analysisTrack = (AliAnalysisTrack *)analysisTrackArray->At(itrk); if (!analysisTrack) continue; /* check accepted track */ if (!analysisTrack->AcceptTrack()) continue; /* check rapidity */ if (dorapidityCut) if (analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift > rapidityMaxCut || analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift < rapidityMinCut) continue; /* check charge */ charge = analysisTrack->GetSign() > 0. ? kPositive : kNegative; if (charge != icharge) continue; /*** ACCEPTED TRACK ***/ /* get track info */ p = analysisTrack->GetP(); pt = analysisTrack->GetPt(); /* compute track mt */ mt = TMath::Sqrt(pt * pt + AliPID::ParticleMass(ipart) * AliPID::ParticleMass(ipart)); /* get TPC info */ dedx = analysisTrack->GetTPCdEdx(); dedxN = analysisTrack->GetTPCdEdxN(); ptpc = analysisTrack->GetTPCmomentum(); /* TPC signal */ bethe = tpcResponse->GetExpectedSignal(ptpc, iipart); /* fix electron signal */ if (iipart == AliPID::kElectron) bethe += 23.; deltadedx = dedx - bethe; dedx_sigma = 0.07 * bethe; tpcsignal = deltadedx / dedx_sigma; /* electronCut requested, remove electrons */ if (electronCut) { /* TPC signal */ bethe = tpcResponse->GetExpectedSignal(ptpc, AliPID::kElectron); /* fix electron signal */ bethe += 23.; deltadedx = dedx - bethe; dedx_sigma = 0.07 * bethe; tpcsignal = deltadedx / dedx_sigma; if (TMath::Abs(tpcsignal) < 1.5) continue; } /* cut on TPC signal if requested */ if (cutOnTPC && (tpcsignal < tpcsignalMin || tpcsignal > tpcsignalMax)) continue; /* fill histos */ hAcceptedTracks->Fill(cent, pt); /* set TOF pid flag */ hastofpid = analysisTrack->HasTOFPID(hEnabledFlag); /* check TOF pid */ if (!hastofpid) continue; /*** ACCEPTED TRACK WITH TOF PID ***/ /* apply expected time correction */ // analysisTrack->ApplyTOFExpectedTimeCorrection(); /* get TOF info */ index = analysisTrack->GetTOFIndex(); time = analysisTrack->GetTOFTime() - TOF_mean; time_sigma = tofReso; /* TZEROTOF */ TZEROTOF = analysisEvent->GetTimeZeroTOF(analysisTrack->GetP()); TZEROTOF_reso = analysisEvent->GetTimeZeroTOFSigma(analysisTrack->GetP()); hasTZEROTOF = TZEROTOF_reso < 150.; if (hasTZEROTOF) { TZEROTOF += TOF_TZEROTOF_mean - TOF_mean; TZEROTOF_reso *= TZEROTOF_resoScaleFactor; } /* TZEROA */ TZEROA = analysisEvent->GetTimeZeroT0(1) - TZEROA_shift; // TZEROA_reso = TZEROA_resolution[icent]; hasTZEROA = TMath::Abs(TZEROA - TZEROA_mean) < 3. * TZEROA_sigma; TZEROA += TOF_TZEROA_mean - TOF_mean; /* TZEROC */ TZEROC = analysisEvent->GetTimeZeroT0(2) - TZEROC_shift; // TZEROC_reso = TZEROC_resolution[icent]; hasTZEROC = TMath::Abs(TZEROC - TZEROC_mean) < 3. * TZEROC_sigma; TZEROC += TOF_TZEROC_mean - TOF_mean; /* TZEROAND */ TZEROAND = (TZEROA + TZEROC) * 0.5; // TZEROAND_reso = TZEROAND_resolution[icent]; hasTZEROAND = hasTZEROA && hasTZEROC; /* TZEROMEAN */ TZEROMEAN = TZEROTOF / TZEROTOF_reso / TZEROTOF_reso; TZEROMEAN_weight = 1. / TZEROTOF_reso / TZEROTOF_reso; if (hasTZEROAND) { // printf("TZEROAND\n"); TZEROMEAN += TZEROAND / TZEROAND_reso / TZEROAND_reso; TZEROMEAN_weight = 1. / TZEROAND_reso / TZEROAND_reso; } else if (hasTZEROA) { // printf("TZEROA\n"); TZEROMEAN += TZEROA / TZEROA_reso / TZEROA_reso; TZEROMEAN_weight = 1. / TZEROA_reso / TZEROA_reso; } else if (hasTZEROC) { // printf("TZEROC\n"); TZEROMEAN += TZEROC / TZEROC_reso / TZEROC_reso; TZEROMEAN_weight = 1. / TZEROC_reso / TZEROC_reso; } timezero = TZEROMEAN / TZEROMEAN_weight; timezero_sigma = TMath::Sqrt(1. / TZEROMEAN_weight); /* TZEROBEST */ TZEROBEST = TZEROTOF; TZEROBEST_reso = TZEROTOF_reso; if (hasTZEROAND && TZEROAND_reso < TZEROBEST_reso) { TZEROBEST = TZEROAND; TZEROBEST_reso = TZEROAND_reso; } else if (hasTZEROA && TZEROA_reso < TZEROBEST_reso) { TZEROBEST = TZEROA; TZEROBEST_reso = TZEROA_reso; } if (hasTZEROC && TZEROC_reso < TZEROBEST_reso) { TZEROBEST = TZEROC; TZEROBEST_reso = TZEROC_reso; } timezero = TZEROBEST; timezero_sigma = TZEROBEST_reso; /* DEBUG */ // timezero = 0.;//TZEROTOF; // timezero_sigma = 203.854691;//TZEROTOF_reso; // if (timezero == 0.) // printf("%f %f\n", timezero, timezero_sigma); timezero_sigma *= scaletimezerosigma; if (resetTZERO) { timezero = 0.; timezero_sigma = timezero_spread; } tof = time - timezero; tof_sigma = TMath::Sqrt(time_sigma * time_sigma + timezero_sigma * timezero_sigma); /* TOF expected time */ texp = analysisTrack->GetTOFExpTime(iipart); texp_sigma = analysisTrack->GetTOFExpTimeSigma(iipart) * scaletexpreso[iipart]; /* TOF signal */ deltat = tof - texp; deltat_sigma = TMath::Sqrt(tof_sigma * tof_sigma + texp_sigma * texp_sigma); tofsignal = deltat / deltat_sigma; /* fill histo */ hTOFpid->Fill(cent, pt, tofsignal); hTOFpid_delta->Fill(cent, p, deltat); hTOFreso->Fill(cent, pt, deltat_sigma); if (hasTZEROTOF || hasTZEROA || hasTZEROC || hasTZEROAND) hTOFpid_withTZERO->Fill(cent, pt, tofsignal); /*** EXPECTED MISMATCH ***/ /* loop to generated random hits */ for (Int_t irnd = 0; irnd < NmismatchTrials; irnd++) { /* generate ramdom tof values */ tof_rnd = GenerateRandomHit(hT0Fill, t0fill, index); /* TOF signal */ deltat = tof_rnd - texp; tofsignal = deltat / deltat_sigma; /* fill histo */ hTOFmismatch->Fill(cent, pt, tofsignal); hTOFmismatch_delta->Fill(cent, p, deltat); } /* end of loop over generated random hits */ /*** EXPECTED SIGNALS ***/ /* loop over other particles */ for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) { /* generate expected tof value */ tof_th = analysisTrack->GetTOFExpTime(iiipart); texp_sigma = analysisTrack->GetTOFExpTimeSigma(iiipart) * scaletexpreso[iiipart]; /* loop over many trials */ for (Int_t irnd = 0; irnd < NexpectedTrials; irnd++) { /* tof response smearing */ signal_smear = fTOFsignal->GetRandom(); /* timezero resolution smearing */ timezero_smear = gRandom->Gaus(0., timezero_sigma); /* texp resolution smearing */ texp_smear = gRandom->Gaus(0., texp_sigma); /* deltat and sigma */ deltat = tof_th - texp + signal_smear + timezero_smear + texp_smear; tofsignal = deltat / deltat_sigma; /* fill histo */ hTOFexpected[iiipart]->Fill(cent, pt, tofsignal); hTOFexpected_delta[iiipart]->Fill(cent, p, deltat); } /* end of loop over many trials */ } /* end of loop over other particle */ } /* end of loop over tracks */ } /* end of loop over events */ /* stop stopwatch */ timer.Stop(); timer.Print(); /* output */ TString outputstring = "TOFpid"; if (!dorapidityCut) outputstring += "_noRapidityCut"; if (electronCut) outputstring += "_electronCut"; if (cutOnTPC) outputstring += Form("_cutOnTPC[%3.1f,%3.1f]", , tpcsignalMin, tpcsignalMax); outputstring += Form("_%s_%s_%sID.%s", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), filename); TFile *fileout = TFile::Open(outputstring.Data(), "RECREATE"); hAcceptedEvents->Write(); hAcceptedTracks->Write(); hTOFpid->Write(); hTOFreso->Write(); hTOFpid_withTZERO->Write(); hTOFmismatch->Write(); for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) hTOFexpected[iiipart]->Write(); hTOFpid_delta->Write(); hTOFmismatch_delta->Write(); for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) hTOFexpected_delta[iiipart]->Write(); fileout->Close(); } //___________________________________________________________________________________ TOFspectra_defaultFit(const Char_t *filename) { Bool_t EXPECTED_SIGNAL_TEMPLATE = kTRUE; Bool_t EXPECTED_SIGNAL_FIT = kFALSE; Bool_t EXPECTED_BACKGROUND_TEMPLATE = kFALSE; Bool_t EXPECTED_BACKGROUND_FIT = kTRUE; } //___________________________________________________________________________________ TOFspectra_defaultFit_fitElectrons(const Char_t *filename, Float_t electronLimit = 5.) { TOFspectra(filename, electronLimit); } //___________________________________________________________________________________ TOFspectra_signalFit(Bool_t fixParams = kTRUE, Float_t scaleSigma = 1., Float_t scaleTail = 1.) { EXPECTED_SIGNAL_TEMPLATE = kFALSE; EXPECTED_SIGNAL_FIT = kTRUE; FIX_SIGNAL_MEAN = fixParams; FIX_SIGNAL_SIGMA = fixParams; FIX_SIGNAL_TAIL = fixParams; SCALE_SIGNAL_SIGMA = scaleSigma; SCALE_SIGNAL_TAIL = scaleTail; } //___________________________________________________________________________________ TOFspectra_bkgFit(Bool_t fixParams = kTRUE, Float_t scaleSigma = 1., Float_t scaleTail = 1.) { EXPECTED_BACKGROUND_TEMPLATE = kFALSE; EXPECTED_BACKGROUND_FIT = kTRUE; FIX_BACKGROUND_MEAN = fixParams; FIX_BACKGROUND_SIGMA = fixParams; FIX_BACKGROUND_TAIL = fixParams; SCALE_BACKGROUND_SIGMA = scaleSigma; SCALE_BACKGROUND_TAIL = scaleTail; TOFspectra(filename); } //___________________________________________________________________________________ void TOFspectra(const Char_t *filename, Float_t electronLimit = 0.) { for (Int_t icent = 0; icent < NcentralityBins; icent++) for (Int_t icharge = 0; icharge < kNCharges; icharge++) for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) TOFspectrum(filename, ipart, icharge, ipart, icent, -1., -1., electronLimit); } //___________________________________________________________________________________ #define USENSIGMA 1 /* fit ranges */ Double_t fitPtMin[AliPID::kSPECIES] = {0.5, 0.5, 0.3, 0.4, 0.5}; Double_t fitPtMax[AliPID::kSPECIES] = {3.0, 3.0, 3.0, 3.0, 5.0}; #if USENSIGMA Double_t fitSigmaMin[AliPID::kSPECIES] = {tofsigmaMin, tofsigmaMin, -25., -75., -65.}; Double_t fitSigmaMax[AliPID::kSPECIES] = {tofsigmaMax, tofsigmaMax, 225., 200., 100.}; #else Double_t fitSigmaMin[AliPID::kSPECIES] = {-24400., -24400., -24400., -24400., -24400.}; Double_t fitSigmaMax[AliPID::kSPECIES] = {24400., 24400., 24400., 24400., 24400.}; #endif void 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) { printf("****************************************\n"); printf("RUNNING SPECTRA FIT:\n"); printf("RAPIDITY-CUT: %s\n", AliPID::ParticleName(ipart)); printf("RAPIDITY-CUT: %s\n", AliPID::ParticleName(ipart)); printf("CHARGE: %s\n", chargeName[icharge]); printf("PARTICLE: %s\n", AliPID::ParticleName(iipart)); printf("CENTRALITY BIN: %d\n", icent); printf("****************************************\n"); /* mismatch balance functions */ TF1 *fMismatchBalanceFunction[5][2]; Double_t mismatchBalanceParam0[5][2] = { 0., 0., 0., 0., 0.669488, 0.599374, -1.2582, -2.53613, 5.48850e-01, 10.1622 }; Double_t mismatchBalanceParam1[5][2] = { 0., 0., 0., 0., -1.64894, -0.825764, -1.63683, -2.55156, -1.64894, -7.09531 }; for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) { for (Int_t iiicharge = 0; iiicharge < 2; iiicharge++) { fMismatchBalanceFunction[iiipart][iiicharge] = new TF1(Form(""), "1. - [0]*TMath::Exp(x*[1])", 0., 5.); fMismatchBalanceFunction[iiipart][iiicharge]->SetParameter(0, mismatchBalanceParam0[iiipart][iiicharge]); fMismatchBalanceFunction[iiipart][iiicharge]->SetParameter(1, mismatchBalanceParam1[iiipart][iiicharge]); } } /* open data */ TFile *filein = TFile::Open(filename); /* get number of events */ TH1F *hAcceptedEvents = (TH1F *)filein->Get(Form("hAcceptedEvents_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); if (!hAcceptedEvents) { printf("cannot find %s\n", Form("hAcceptedEvents_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); return; } Double_t nevents; if (icent < 0 || icent >= NcentralityBins) nevents = hAcceptedEvents->GetEntries(); else nevents = hAcceptedEvents->Integral(icent + 1, icent + 1); printf("N_EVENTS : %d\n", nevents); printf("****************************************\n"); /* get histos */ TH2I *hAcceptedTracks = (TH2I *)filein->Get(Form("hAcceptedTracks_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); if (!hAcceptedTracks) { printf("cannot find %s\n", Form("hAcceptedTracks_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); // return; } #if USENSIGMA TH3I *hTOFpid = (TH3I *)filein->Get(Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); #else TH3I *hTOFpid = (TH3I *)filein->Get(Form("hTOFpid_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); #endif if (!hTOFpid) { printf("cannot find %s\n", Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); return; } #if USENSIGMA TH3I *hTOFmismatch = (TH3I *)filein->Get(Form("hTOFmismatch_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); #else TH3I *hTOFmismatch = (TH3I *)filein->Get(Form("hTOFmismatch_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); #endif if (!hTOFmismatch) { printf("cannot find %s\n", Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart))); return; } TH3I *hTOFexpected[AliPID::kSPECIES]; for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) { #if USENSIGMA hTOFexpected[iiipart] = (TH3I *)filein->Get(Form("hTOFexpected_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart))); #else hTOFexpected[iiipart] = (TH3I *)filein->Get(Form("hTOFexpected_delta_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart))); #endif if (!hTOFexpected[iiipart]) { printf("cannot find %s\n", Form("hTOFexpected_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart))); return; } } /* setup centrality range */ if (icent < 0 || icent >= NcentralityBins) { printf("WARNING: undefined centrality -> using 00-90\% range\n"); if (hAcceptedTracks) hAcceptedTracks->GetXaxis()->SetRange(1, NcentralityBins); hTOFpid->GetXaxis()->SetRange(1, NcentralityBins); hTOFmismatch->GetXaxis()->SetRange(1, NcentralityBins); for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) hTOFexpected[iiipart]->GetXaxis()->SetRange(1, NcentralityBins); } else { printf("***** FITTING CENTRALITY-BIN [%02d, %02d] %% *****\n", centralityBin[icent], centralityBin[icent + 1]); if (hAcceptedTracks) hAcceptedTracks->GetXaxis()->SetRange(icent + 1, icent + 1); hTOFpid->GetXaxis()->SetRange(icent + 1, icent + 1); hTOFmismatch->GetXaxis()->SetRange(icent + 1, icent + 1); for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) hTOFexpected[iiipart]->GetXaxis()->SetRange(icent + 1, icent + 1); } /* init flags */ Bool_t requestedRange = kFALSE; Bool_t fitElectrons = kTRUE; Bool_t fitMuons = kTRUE; Bool_t fitPions = kTRUE; /* setup pt range if requested */ if (ptMin > -0.001 && ptMax > -0.001 && ptMax > ptMin) { printf("***** FITTING PT-BIN [%f, %f] GeV/c *****\n", ptMin, ptMax); requestedRange = kTRUE; /* check electron-fit is allowed */ fitElectrons = kTRUE; if ((ptMin + 0.001) < FIT_ELECTRONS_PT_MIN || (ptMax - 0.001) > FIT_ELECTRONS_PT_MAX) fitElectrons = kFALSE; /* check muon-fit is allowed */ fitMuons = kTRUE; if ((ptMin + 0.001) < FIT_MUONS_PT_MIN || (ptMax - 0.001) > FIT_MUONS_PT_MAX) fitMuons = kFALSE; /* check pion-fit is allowed */ fitPions = kTRUE; if ((ptMin + 0.001) < FIT_PIONS_PT_MIN || (ptMax - 0.001) > FIT_PIONS_PT_MAX) fitPions = kFALSE; hTOFpid->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001); hTOFmismatch->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001); for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) hTOFexpected[iiipart]->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001); } /* output */ Char_t outfilename[1024]; if (icent < 0 || icent >= NcentralityBins) sprintf(outfilename, "TOFspectrum_cent0090_%s_%s_%sID.root", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)); else { sprintf(outfilename, "TOFspectrum_cent%02d%02d_%s_%s_%sID.root", centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)); } TFile *fileout = TFile::Open(outfilename, "RECREATE"); TDirectory *fitDir = fileout->mkdir("FitParams"); /* canvas */ TCanvas *canvas = new TCanvas("canvas"); canvas->SetLogy(); /* histo */ TH1D *hFitParamHisto[kNFitParams]; for (Int_t iparam = 0; iparam < kNFitParams; iparam++) hFitParamHisto[iparam] = new TH1D(Form("h%s", fitParamName[iparam]), fitParamTitle[iparam], NptBins, ptBin); /* loop over ptBins */ for (Int_t ipt = 0; ipt < NptBins; ipt++) { if (!requestedRange) { if ((ptBin[ipt] + 0.001) < fitPtMin[ipart] || (ptBin[ipt + 1] - 0.001) > fitPtMax[ipart]) continue; printf("***** FITTING PT-BIN [%f, %f] GeV/c *****\n", ptBin[ipt], ptBin[ipt + 1]); /* check electron-fit is allowed */ fitElectrons = kTRUE; if ((ptBin[ipt] + 0.001) < FIT_ELECTRONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_ELECTRONS_PT_MAX) fitElectrons = kFALSE; /* check muon-fit is allowed */ fitMuons = kTRUE; if ((ptBin[ipt] + 0.001) < FIT_MUONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_MUONS_PT_MAX) fitMuons = kFALSE; /* check pion-fit is allowed */ fitPions = kTRUE; if ((ptBin[ipt] + 0.001) < FIT_PIONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_PIONS_PT_MAX) fitPions = kFALSE; hTOFpid->GetYaxis()->SetRange(ipt + 1, ipt + 1); hTOFmismatch->GetYaxis()->SetRange(ipt + 1, ipt + 1); for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) hTOFexpected[iiipart]->GetYaxis()->SetRange(ipt + 1, ipt + 1); } /* nsigma projections */ TH1 *hSignal_py = hTOFpid->Project3D("z"); TH1 *hMismatch_py = hTOFmismatch->Project3D("z"); TH1 *hSignalExp_py[AliPID::kSPECIES]; TH1 *hSignalExpTail_py[AliPID::kSPECIES]; for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) { hSignalExp_py[iiipart] = hTOFexpected[iiipart]->Project3D("z"); hSignalExpTail_py[iiipart] = hTOFexpected[iiipart]->Project3D("z"); } /* prepare histos for the fitter */ Int_t partbkg1[AliPID::kSPECIES] = {AliPID::kKaon, AliPID::kKaon, AliPID::kKaon, AliPID::kPion, AliPID::kPion}; Int_t partbkg2[AliPID::kSPECIES] = {AliPID::kProton, AliPID::kProton, AliPID::kProton, AliPID::kProton, AliPID::kKaon}; Int_t partbkg3[AliPID::kSPECIES] = {AliPID::kPion, AliPID::kPion, AliPID::kElectron, AliPID::kElectron, AliPID::kElectron}; Int_t partbkg4[AliPID::kSPECIES] = {AliPID::kMuon, AliPID::kElectron, AliPID::kMuon, AliPID::kMuon, AliPID::kMuon}; TH1 *hSigExp_py, *hBkgExp1_py, *hBkgExp2_py, *hBkgExp3_py, *hBkgExp4_py; hSigExp_py = EXPECTED_SIGNAL_TAIL ? hSignalExpTail_py[iipart] : hSignalExp_py[iipart]; hBkgExp1_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg1[iipart]] : hSignalExp_py[partbkg1[iipart]]; hBkgExp2_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg2[iipart]] : hSignalExp_py[partbkg2[iipart]]; hBkgExp3_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg3[iipart]] : hSignalExp_py[partbkg3[iipart]]; hBkgExp4_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg4[iipart]] : hSignalExp_py[partbkg4[iipart]]; /* check histos if requested */ if (checkHistoFlag) { TCanvas *cCheckHisto = new TCanvas("cCheckHisto"); cCheckHisto->Divide(2, 3); cCheckHisto->cd(1); hSignal_py->Draw(); cCheckHisto->cd(2); hSigExp_py->Draw(); cCheckHisto->cd(3); hBkgExp1_py->Draw(); cCheckHisto->cd(4); hBkgExp2_py->Draw(); cCheckHisto->cd(5); hBkgExp3_py->Draw(); cCheckHisto->cd(6); hMismatch_py->Draw(); return; } Double_t rangeMin = fitSigmaMin[iipart], rangeMax = fitSigmaMax[iipart]; Bool_t constrainSignal = kFALSE; Bool_t constrainBkg1 = kFALSE; Bool_t constrainBkg2 = kFALSE; Bool_t forceGaussianSignal = kFALSE; Bool_t fitBkg1 = kTRUE, fitBkg2 = kTRUE, fitBkg3 = kTRUE, fitBkg4 = kTRUE; /* check whether we can fit electrons */ if (!fitElectrons) { printf("INHIBIT FIT ELECTRONS\n"); if (partbkg1[iipart] == AliPID::kElectron) fitBkg1 = kFALSE; if (partbkg2[iipart] == AliPID::kElectron) fitBkg2 = kFALSE; if (partbkg3[iipart] == AliPID::kElectron) fitBkg3 = kFALSE; if (partbkg4[iipart] == AliPID::kElectron) fitBkg4 = kFALSE; } /* check whether we can fit muons */ if (!fitMuons) { printf("INHIBIT FIT MUONS\n"); if (partbkg1[iipart] == AliPID::kMuon) fitBkg1 = kFALSE; if (partbkg2[iipart] == AliPID::kMuon) fitBkg2 = kFALSE; if (partbkg3[iipart] == AliPID::kMuon) fitBkg3 = kFALSE; if (partbkg4[iipart] == AliPID::kMuon) fitBkg4 = kFALSE; } /* check whether we can fit pions */ if (!fitPions) { printf("INHIBIT FIT PIONS\n"); if (partbkg1[iipart] == AliPID::kPion) fitBkg1 = kFALSE; if (partbkg2[iipart] == AliPID::kPion) fitBkg2 = kFALSE; if (partbkg3[iipart] == AliPID::kPion) fitBkg3 = kFALSE; if (partbkg4[iipart] == AliPID::kPion) fitBkg4 = kFALSE; } /* fit */ Double_t param[kNFitParams]; Double_t param_err[kNFitParams]; 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); /* if fitting pions add electrons and to the signal */ if (iipart == 2) { param[kSignalCounts] += param[kBkg3Counts] + param[kBkg4Counts]; param_err[kSignalCounts] = TMath::Sqrt(param_err[kSignalCounts] * param_err[kSignalCounts] + param_err[kBkg3Counts] * param_err[kBkg3Counts] + param_err[kBkg4Counts] * param_err[kBkg3Counts]); } /* else add to bkg1 (pions) */ else { param[kBkg1Counts] += param[kBkg3Counts] + param[kBkg4Counts]; param_err[kBkg1Counts] = TMath::Sqrt(param_err[kBkg1Counts] * param_err[kBkg1Counts] + param_err[kBkg3Counts] * param_err[kBkg3Counts] + param_err[kBkg4Counts] * param_err[kBkg3Counts]); } /* check requested pt-range */ if (requestedRange) return; /* write canvas */ fitDir->cd(); canvas->Write(Form("fitDisplay_ptBin_%3.2f_%3.2f", ptBin[ipt], ptBin[ipt + 1])); /* set histo */ for (Int_t iparam = 0; iparam < kNFitParams; iparam++) { hFitParamHisto[iparam]->SetBinContent(ipt + 1, param[iparam]); hFitParamHisto[iparam]->SetBinError(ipt + 1, param_err[iparam]); } /* delete */ delete hSignal_py; delete hMismatch_py; for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) { delete hSignalExp_py[iiipart]; // delete hSignalExpTail_py[iiipart]; } } /* check requested pt-range */ if (requestedRange) return; /*** POST-ANALYSIS ***/ TDirectory *postDir = fileout->mkdir("PostAnalysis"); /* compute overflows */ TH1D *hOverflowCounts = new TH1D(*hFitParamHisto[kTotalCounts]); hOverflowCounts->SetNameTitle("hOverflowCounts", "Overflow counts: TotalCounts - IntegralCounts;p_{T} (GeV/c);counts"); hOverflowCounts->Add(hFitParamHisto[kIntegralCounts], -1.); /* compute total mismatch counts */ TH1D *hTotalMismatchCounts = new TH1D(*hFitParamHisto[kMismatchCounts]); hTotalMismatchCounts->SetNameTitle("hTotalMismatchCounts", "Total mismatch counts: MismatchCounts + OverflowCounts;p_{T} (GeV/c);counts"); hTotalMismatchCounts->Add(hOverflowCounts); /* computed mismatch fraction */ TH1D *hTotalMismatchFraction = new TH1D(*hTotalMismatchCounts); hTotalMismatchFraction->SetNameTitle("hTotalMismatchFraction", "Total mismatch fraction: TotalMismatchCounts / TotalCounts;p_{T} (GeV/c);"); hTotalMismatchFraction->Divide(hFitParamHisto[kTotalCounts]); /* compute identified counts */ TH1D *hIdentifiedCounts = new TH1D(*hFitParamHisto[kSignalCounts]); hIdentifiedCounts->SetNameTitle("hIdentifiedCounts", "Identified counts: SignalCounts + sum(BkgCounts);p_{T} (GeV/c);counts"); hIdentifiedCounts->Add(hFitParamHisto[kBkg1Counts]); hIdentifiedCounts->Add(hFitParamHisto[kBkg2Counts]); // hIdentifiedCounts->Add(hFitParamHisto[kBkg3Counts]); // hIdentifiedCounts->Add(hFitParamHisto[kBkg4Counts]); /* compute signal fraction */ TH1D *hSignalFraction = new TH1D(*hFitParamHisto[kSignalCounts]); hSignalFraction->SetNameTitle("hSignalFraction", "Signal fraction: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);"); hSignalFraction->Divide(hSignalFraction, hIdentifiedCounts, 1., 1., "B"); /* compute bkg1 fraction */ TH1D *hBkg1Fraction = new TH1D(*hFitParamHisto[kBkg1Counts]); hBkg1Fraction->SetNameTitle("hBkg1Fraction", "Bkg1 fraction: Bkg1Counts / IdentifiedCounts;p_{T} (GeV/c);"); hBkg1Fraction->Divide(hBkg1Fraction, hIdentifiedCounts, 1., 1., "B"); /* compute bkg2 fraction */ TH1D *hBkg2Fraction = new TH1D(*hFitParamHisto[kBkg2Counts]); hBkg2Fraction->SetNameTitle("hBkg2Fraction", "Bkg2 fraction: Bkg2Counts / IdentifiedCounts;p_{T} (GeV/c);"); hBkg2Fraction->Divide(hBkg2Fraction, hIdentifiedCounts, 1., 1., "B"); /* compute bkg3 fraction */ TH1D *hBkg3Fraction = new TH1D(*hFitParamHisto[kBkg3Counts]); hBkg3Fraction->SetNameTitle("hBkg3Fraction", "Bkg3 fraction: Bkg3Counts / IdentifiedCounts;p_{T} (GeV/c);"); hBkg3Fraction->Divide(hBkg3Fraction, hIdentifiedCounts, 1., 1., "B"); /* compute bkg4 fraction */ TH1D *hBkg4Fraction = new TH1D(*hFitParamHisto[kBkg4Counts]); hBkg4Fraction->SetNameTitle("hBkg4Fraction", "Bkg4 fraction: Bkg4Counts / IdentifiedCounts;p_{T} (GeV/c);"); hBkg4Fraction->Divide(hBkg4Fraction, hIdentifiedCounts, 1., 1., "B"); /* rework */ TH1D *hSignalFractionReworked = new TH1D(*hFitParamHisto[kSignalCounts]); hSignalFractionReworked->SetNameTitle("hSignalFractionReworked", "Signal fraction reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);"); hSignalFractionReworked->Multiply(fMismatchBalanceFunction[iipart][icharge]); TH1D *hBkg1FractionReworked = new TH1D(*hFitParamHisto[kBkg1Counts]); hBkg1FractionReworked->SetNameTitle("hBkg1FractionReworked", "Bkg1 fraction reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);"); hBkg1FractionReworked->Multiply(fMismatchBalanceFunction[partbkg1[iipart]][icharge]); TH1D *hBkg2FractionReworked = new TH1D(*hFitParamHisto[kBkg2Counts]); hBkg2FractionReworked->SetNameTitle("hBkg2FractionReworked", "Bkg1 fraction reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);"); hBkg2FractionReworked->Multiply(fMismatchBalanceFunction[partbkg2[iipart]][icharge]); TH1D *hIdentifiedCountsReworked = new TH1D(*hSignalFractionReworked); hIdentifiedCountsReworked->SetNameTitle("hIdentifiedCountsReworked", "Identified counts reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);"); hIdentifiedCountsReworked->Add(hBkg1FractionReworked); hIdentifiedCountsReworked->Add(hBkg2FractionReworked); hSignalFractionReworked->Divide(hSignalFractionReworked, hIdentifiedCountsReworked, 1., 1., "B"); /* compute mismatch-correction counts */ TH1D *hMismatchCorrectionCounts = new TH1D(*hTotalMismatchCounts); hMismatchCorrectionCounts->SetNameTitle("hMismatchCorrectionCounts", "Mismatch-correction counts: TotalMismatchCounts * SignalFraction;p_{T} (GeV/c);counts"); hMismatchCorrectionCounts->Multiply(hSignalFractionReworked); /* compute mismatch-corrected signal counts */ TH1D *hMismatchCorrectedSignalCounts = new TH1D(*hFitParamHisto[kSignalCounts]); hMismatchCorrectedSignalCounts->SetNameTitle("hMismatchCorrectedSignalCounts", "Mismatch-corrected signal counts: SignalCounts + MismatchCorrectionCounts;p_{T} (GeV/c);counts"); hMismatchCorrectedSignalCounts->Add(hMismatchCorrectionCounts); /* accepted tracks histo */ if (hAcceptedTracks) { TH1D *hAcceptedTracks_py = hAcceptedTracks->ProjectionY(); hAcceptedTracks_py->SetNameTitle("hAcceptedTracks", "Accepted tracks;p_{T} (GeV/c);"); hAcceptedTracks_py->Sumw2(); } /*** RAW SPECTRA ***/ TDirectory *rawDir = fileout->mkdir("RawSpectra"); /* compute normalized raw yield */ TH1D *hNormalizedRawYield = new TH1D(*hFitParamHisto[kSignalCounts]); hNormalizedRawYield->SetNameTitle("hNormalizedRawYield", "Raw yield;p_{T} (GeV/c);#frac{1}{N_{ev}} #frac{d^{2}N}{dy dp_{T}}"); TOFpid_normalize(hNormalizedRawYield, nevents); /* compute normalized mismatch-corrected raw yield */ TH1D *hNormalizedMismatchCorrectedRawYield = new TH1D(*hMismatchCorrectedSignalCounts); hNormalizedMismatchCorrectedRawYield->SetNameTitle("hNormalizedMismatchCorrectedRawYield", "Mismatch-corrected raw yield;p_{T} (GeV/c);#frac{1}{N_{ev}} #frac{d^{2}N}{dy dp_{T}}"); TOFpid_normalize(hNormalizedMismatchCorrectedRawYield, nevents); /*** OUTPUT ***/ /* write fir params histos */ fitDir->cd(); for (Int_t iparam = 0; iparam < kNFitParams; iparam++) hFitParamHisto[iparam]->Write(); /* write post-analysis histos */ postDir->cd(); hOverflowCounts->Write(); hTotalMismatchCounts->Write(); hTotalMismatchFraction->Write(); hIdentifiedCounts->Write(); hSignalFraction->Write(); hSignalFractionReworked->Write(); hBkg1Fraction->Write(); hBkg2Fraction->Write(); hBkg3Fraction->Write(); hBkg4Fraction->Write(); hMismatchCorrectionCounts->Write(); hMismatchCorrectedSignalCounts->Write(); if (hAcceptedTracks) hAcceptedTracks_py->Write(); /* write raw spectra histos */ rawDir->cd(); hNormalizedRawYield->Write(); hNormalizedMismatchCorrectedRawYield->Write(); /* clean up */ delete canvas; for (Int_t iparam = 0; iparam < kNFitParams; iparam++) delete hFitParamHisto[iparam]; delete hOverflowCounts; delete hTotalMismatchCounts; delete hTotalMismatchFraction; delete hIdentifiedCounts; delete hSignalFraction; delete hBkg1Fraction; delete hBkg2Fraction; delete hBkg3Fraction; delete hBkg4Fraction; delete hMismatchCorrectionCounts; delete hMismatchCorrectedSignalCounts; if (hAcceptedTracks) { delete hAcceptedTracks_py; } delete hNormalizedRawYield; delete hNormalizedMismatchCorrectedRawYield; /* close file */ fileout->Close(); } //___________________________________________________________________________________ Float_t TOFpid_histomin(TH1 *h) { for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++) if (h->GetBinContent(ibin + 1) > 0.) return h->GetXaxis()->GetBinCenter(ibin + 1); return kMaxInt; } //___________________________________________________________________________________ Float_t TOFpid_histomax(TH1 *h) { for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) if (h->GetBinContent(ibin) > 0.) return h->GetXaxis()->GetBinCenter(ibin); return -kMaxInt; } //___________________________________________________________________________________ TOFpid_checkneg(TH1 *h) { for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++) if (h->GetBinContent(ibin + 1) <= 0.) { h->SetBinContent(ibin + 1, 1.e-300); // h->SetBinError(ibin + 1, 0.1); } } //___________________________________________________________________________________ Double_t 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) { /** ROOFIT ***/ gSystem->Load("libRooFit"); using namespace RooFit; /*** LOAD GAUSSIANTAIL CLASS ***/ gSystem->Load("RooFermiCutoff_cxx.so"); gSystem->Load("RooGaussianTail_cxx.so"); /*** DEFINE FIT RANGE ***/ printf("***** FIT RANGE DEFINITION *****\n"); /* check mismatch histogram to define min/max fit range */ // rangeMin = TMath::Max(rangeMin, TOFpid_histomin(hMismatch)); // rangeMax = TMath::Min(rangeMax, TOFpid_histomax(hMismatch)); /* fix zeroes */ TOFpid_checkneg(hMismatch); /* define range */ RooRealVar x("x", "n_{#sigma}", 0., rangeMin, rangeMax, ""); printf("FIT RANGE DEFINED: %f -> %f\n", rangeMin, rangeMax); printf("********************************\n"); /*** DEFINE HISTOGRAM DATA ***/ /* define data to fit and background from input histogram */ RooDataHist hdata("hdata", "hdata", x, hSignal); RooDataHist hsig("hsig", "hsig", x, hSigExp); RooDataHist hbkg1("hbkg1", "hbkg1", x, hBkgExp1); RooDataHist hbkg2("hbkg2", "hbkg2", x, hBkgExp2); RooDataHist hbkg3("hbkg3", "hbkg3", x, hBkgExp3); RooDataHist hbkg4("hbkg4", "hbkg4", x, hBkgExp4); RooDataHist hmismatch("hmismatch", "hmismatch", x, hMismatch); /*** DEFINE SIGNAL SHAPE ***/ printf("***** SIGNAL SHAPE DEFINITION *****\n"); /* variables */ if (FIX_SIGNAL_MEAN) { RooRealVar mean("mean", "mean", DEFAULT_SIGNAL_MEAN, ""); printf("FIXED SIGNAL_MEAN = %f\n", mean.getVal()); } else { RooRealVar mean("mean", "mean", DEFAULT_SIGNAL_MEAN, MIN_SIGNAL_MEAN, MAX_SIGNAL_MEAN, ""); printf("FREE SIGNAL_MEAN = %f [%f, %f]\n", mean.getVal(), MIN_SIGNAL_MEAN, MAX_SIGNAL_MEAN); } if (FIX_SIGNAL_SIGMA) { RooRealVar sigma("sigma", "sigma", DEFAULT_SIGNAL_SIGMA, ""); printf("FIXED SIGNAL_SIGMA = %f\n", sigma.getVal()); } else { RooRealVar sigma("sigma", "sigma", DEFAULT_SIGNAL_SIGMA, MIN_SIGNAL_SIGMA, MAX_SIGNAL_SIGMA, ""); printf("FREE SIGNAL_SIGMA = %f\n", sigma.getVal()); } if (FIX_SIGNAL_TAIL) { RooRealVar tail("tail", "tail", DEFAULT_SIGNAL_TAIL, ""); printf("FIXED SIGNAL_TAIL = %f\n", tail.getVal()); } else { RooRealVar tail("tail", "tail", DEFAULT_SIGNAL_TAIL, MIN_SIGNAL_TAIL, MAX_SIGNAL_TAIL, ""); printf("FREE SIGNAL_TAIL = %f\n", tail.getVal()); } RooRealVar gaussianfrac("gaussianfrac", "gaussianfrac", 1., 0., 1., ""); RooRealVar sigalpha("sigalpha", "sigalpha", 0., -10., 0.); /* shapes */ if (GAUSSIAN_SIGNAL || forceGaussianSignal) { printf("USING GAUSSIAN_SIGNAL SHAPE\n"); RooGaussian signal("signal", "signal", x, mean, sigma); } else if (GAUSSIANTAIL_SIGNAL) { printf("USING GAUSSIANTAIL_SIGNAL SHAPE\n"); RooGaussianTail signal("signal", "signal", x, mean, sigma, tail); } else if (GAUSSIANTAIL2_SIGNAL) { printf("USING GAUSSIANTAIL2_SIGNAL SHAPE\n"); RooGaussianTail signal("signal", "signal", x, mean, sigma, sigma); } else if (GAUSSIANPLUSGAUSSIANTAIL_SIGNAL) { printf("USING GAUSSIANPLUSGAUSSIANTAIL_SIGNAL SHAPE\n"); RooGaussian gaussian("gaussian", "gaussian", x, mean, sigma); RooGaussianTail gaussiantail("gaussiantail", "gaussiantail", x, mean, sigma, tail); RooAddPdf signal("signal", "signal", RooArgList(gaussian, gaussiantail), gaussianfrac, kTRUE); } else if (GAUSSIANPLUSEXPONENTIAL_SIGNAL) { printf("USING GAUSSIANPLUSEXPONENTIAL_SIGNAL SHAPE\n"); RooGaussian gaussian("gaussian", "gaussian", x, mean, sigma); RooExponential sigexpo("sigexpo", "sigexpo", x, sigalpha); RooRealVar sigcutoff("sigcutoff", "sigcutoff", 0.); RooRealVar sigpower("sigpower", "sigpower", 0.1); RooFermiCutoff sigfermi("sigfermi", "sigfermi", x, sigcutoff, sigpower); RooProdPdf exposignal("exposignal", "exposignal", sigfermi, sigexpo, -100.); RooRealVar gaussianfrac("gaussianfrac", "gaussianfrac", 0.9, 0.7, 1., ""); RooAddPdf signal("signal", "signal", RooArgList(gaussian, exposignal), gaussianfrac, kTRUE); } else if (EXPECTED_SIGNAL_TEMPLATE) { printf("SHAPE OF SIGNAL FROM TEMPLATE HISTOGRAM\n"); RooHistPdf signal("signal", "signal", x, hsig); } else if (EXPECTED_SIGNAL_FIT) { /* fitting bkg1 */ TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus"); hSigExp->Fit(fGaus, "0q"); Double_t sig_mean = fGaus->GetParameter(1); Double_t sig_sigma = fGaus->GetParameter(2); mean.setConstant(kFALSE); mean.setVal(sig_mean); mean.setRange(sig_mean - 10., sig_mean + 10.); mean.setConstant(kFALSE); sigma.setConstant(kFALSE); sigma.setVal(sig_sigma); sigma.setRange(sig_sigma * 0.5, sig_sigma * 1.5); tail.setConstant(kFALSE); tail.setVal(1.); tail.setRange(0.5, 5.); RooGaussianTail signal("signal", "signal", x, mean, sigma, tail); signal.fitTo(hsig, Range(sig_mean - 5. * sig_sigma, sig_mean + 10. * sig_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10)); #if DISPLAY TCanvas *cSig_fit = new TCanvas("cSig_fit"); RooPlot *sig_frame = x.frame(); hsig.plotOn(sig_frame); signal.plotOn(sig_frame, LineColor(kRed)); sig_frame->Draw(); cSig_fit->Update(); #endif printf("SIGNAL PARAMETERS AFTER FIT OF EXPECTED SIGNAL\n"); printf("mean = %f +- %f\n", mean.getVal(), mean.getError()); printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError()); printf("tail = %f +- %f\n", tail.getVal(), tail.getError()); /* scale parameters if requested */ if (SCALE_SIGNAL_SIGMA != 1.) { printf("SCALE FITTED SIGNAL SIGMA BY %f\n", SCALE_SIGNAL_SIGMA); sigma.setVal(sigma.getVal() * SCALE_SIGNAL_SIGMA); } if (SCALE_SIGNAL_TAIL != 1.) { printf("SCALE FITTED SIGNAL TAIL BY %f\n", SCALE_SIGNAL_TAIL); tail.setVal(tail.getVal() * SCALE_SIGNAL_TAIL); } /* fix/release parameters if requested */ if (FIX_SIGNAL_MEAN) { printf("SETTING FITTED SIGNAL MEAN AS CONSTANTS\n"); mean.setConstant(kTRUE); } else { printf("SETTING FITTED SIGNAL MEAN AS FREE\n"); // mean.setRange(mean.getVal() - 0.25 * TMath::Abs(mean.getVal()), mean.getVal() + 0.25 * TMath::Abs(mean.getVal())); // mean.setRange(mean.getVal() - 0.5 * TMath::Abs(mean.getVal()), mean.getVal() + 0.5 * TMath::Abs(mean.getVal())); mean.setRange(mean.getVal() - 0.2, mean.getVal() + 0.2); } if (FIX_SIGNAL_SIGMA) { printf("SETTING FITTED SIGNAL SIGMA AS CONSTANTS\n"); sigma.setConstant(kTRUE); } else { printf("SETTING FITTED SIGNAL SIGMA AS FREE\n"); sigma.setRange(sigma.getVal() * 0.8, sigma.getVal() * 1.2); // sigma.setRange(sigma.getVal() - 0.5, sigma.getVal() + 0.5); } if (FIX_SIGNAL_TAIL) { printf("SETTING FITTED SIGNAL TAIL AS CONSTANTS\n"); tail.setConstant(kTRUE); } else { printf("SETTING FITTED SIGNAL TAIL AS FREE\n"); // tail.setRange(0.5, 5.0); tail.setRange(tail.getVal() * 0.8, tail.getVal() * 1.2); // tail.setRange(tail.getVal() - 0.5, tail.getVal() + 0.5); } } else { printf("SHAPE OF SIGNAL NOT DEFINE: using GAUSSIAN_SIGNAL\n"); RooGaussian signal("signal", "signal", x, mean, sigma); } #if 0 if (constrainSignal) { #if 0 /* fix expected signal and constrain parameters if requested */ signal.fitTo(hsig); #if 0 TCanvas *cConstrainSignal = new TCanvas("cConstrainSignal"); RooPlot *xfr = x.frame(); hsig.plotOn(xfr); signal.plotOn(xfr, LineColor(kRed)); xfr->Draw(); cConstrainSignal->Update(); #endif printf("SIGNAL PARAMETERS AFTER FIT OF EXPECTED SIGNAL\n"); printf("mean = %f +- %f\n", mean.getVal(), mean.getError()); printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError()); printf("tail = %f +- %f\n", tail.getVal(), tail.getError()); if (constrainSignal) { mean.setConstant(kTRUE); sigma.setConstant(kTRUE); tail.setConstant(kTRUE); printf("SIGNAL PARAMETERS CONSTRAINED AFTER FIT OF EXPECTED SIGNAL\n"); } #endif } #endif if (constrainSignal) { mean.setConstant(kTRUE); sigma.setConstant(kTRUE); tail.setConstant(kTRUE); // mean.setRange(-0.1, 0.1); // sigma.setRange(0.95, 1.05); // tail.setRange(0.95, 1.25); } printf("***********************************\n"); /*** DEFINE IDENTIFIED BACKGROUND SHAPES ***/ printf("***** IDENTIFIED BACKGROUND SHAPE DEFINITION *****\n"); /* shapes */ if (EXPECTED_BACKGROUND_TEMPLATE) { printf("USING EXPECTED BACKGROUND TEMPLATE SHAPES FROM HISTOGRAMS\n"); RooHistPdf bkg1("bkg1", "bkg1", x, hbkg1, 0); RooHistPdf bkg2("bkg2", "bkg2", x, hbkg2, 0); RooHistPdf bkg3("bkg3", "bkg3", x, hbkg3, 0); RooHistPdf bkg4("bkg4", "bkg4", x, hbkg4, 0); } else if (EXPECTED_BACKGROUND_FIT) { printf("USING EXPECTED BACKGROUND FITTED SHAPES FROM HISTOGRAMS\n"); /* fitting bkg1 */ TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus"); hBkgExp1->Fit(fGaus, "0q"); Double_t bkgexp1_mean = fGaus->GetParameter(1); Double_t bkgexp1_sigma = fGaus->GetParameter(2); RooRealVar mean_bkg1("mean_bkg1", "mean_bkg1", bkgexp1_mean, bkgexp1_mean - 10., bkgexp1_mean + 10., ""); RooRealVar sigma_bkg1("sigma_bkg1", "sigma_bkg1", bkgexp1_sigma, bkgexp1_sigma * 0.5, bkgexp1_sigma * 1.5, ""); RooRealVar tail_bkg1("tail_bkg1", "tail_bkg1", 1.0, 0.5, 5., ""); RooGaussianTail bkg1("bkg1", "bkg1", x, mean_bkg1, sigma_bkg1, tail_bkg1); bkg1.fitTo(hbkg1, Range(bkgexp1_mean - 5. * bkgexp1_sigma, bkgexp1_mean + 10. * bkgexp1_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10)); #if DISPLAY TCanvas *cBkg1_fit = new TCanvas("cBkg1_fit"); RooPlot *bkg1_frame = x.frame(); hbkg1.plotOn(bkg1_frame); bkg1.plotOn(bkg1_frame, LineColor(kCyan+1)); bkg1_frame->Draw(); cBkg1_fit->Update(); #endif printf("BACKGROUND-1 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-1\n"); printf("mean_bkg1 = %f +- %f\n", mean_bkg1.getVal(), mean_bkg1.getError()); printf("sigma_bkg1 = %f +- %f\n", sigma_bkg1.getVal(), sigma_bkg1.getError()); printf("tail_bkg1 = %f +- %f\n", tail_bkg1.getVal(), tail_bkg1.getError()); /* scale parameters if requested */ if (SCALE_BACKGROUND_SIGMA != 1.) { printf("SCALE FITTED BACKGROUND-1 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA); sigma_bkg1.setVal(sigma_bkg1.getVal() * SCALE_BACKGROUND_SIGMA); } if (SCALE_BACKGROUND_TAIL != 1.) { printf("SCALE FITTED BACKGROUND-1 TAIL BY %f\n", SCALE_BACKGROUND_TAIL); tail_bkg1.setVal(tail_bkg1.getVal() * SCALE_BACKGROUND_TAIL); } /* fix/release parameters if requested */ if (FIX_BACKGROUND_MEAN) { printf("SETTING BACKGROUND-1 FITTED MEAN AS CONSTANTS\n"); mean_bkg1.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-1 FITTED MEAN AS FREE\n"); mean_bkg1.setRange(mean_bkg1.getVal() - 0.25 * TMath::Abs(mean_bkg1.getVal()), mean_bkg1.getVal() + 0.25 * TMath::Abs(mean_bkg1.getVal())); } if (FIX_BACKGROUND_SIGMA) { printf("SETTING BACKGROUND-1 FITTED SIGMA AS CONSTANTS\n"); sigma_bkg1.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-1 FITTED SIGMA AS FREE\n"); sigma_bkg1.setRange(sigma_bkg1.getVal() * 0.75, sigma_bkg1.getVal() * 1.25); } if (FIX_BACKGROUND_TAIL) { printf("SETTING BACKGROUND-1 FITTED TAIL AS CONSTANTS\n"); tail_bkg1.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-1 FITTED TAIL AS FREE\n"); tail_bkg1.setRange(tail_bkg1.getVal() * 0.75, tail_bkg1.getVal() * 1.25); } /* fitting bkg2 */ TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus"); hBkgExp2->Fit(fGaus, "0q"); Double_t bkgexp2_mean = fGaus->GetParameter(1); Double_t bkgexp2_sigma = fGaus->GetParameter(2); RooRealVar mean_bkg2("mean_bkg2", "mean_bkg2", bkgexp2_mean, bkgexp2_mean - 10., bkgexp2_mean + 10., ""); RooRealVar sigma_bkg2("sigma_bkg2", "sigma_bkg2", bkgexp2_sigma, bkgexp2_sigma * 0.5, bkgexp2_sigma * 1.5, ""); RooRealVar tail_bkg2("tail_bkg2", "tail_bkg2", 1.0, 0.5, 5., ""); RooGaussianTail bkg2("bkg2", "bkg2", x, mean_bkg2, sigma_bkg2, tail_bkg2); bkg2.fitTo(hbkg2, Range(bkgexp2_mean - 5. * bkgexp2_sigma, bkgexp2_mean + 10. * bkgexp2_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10)); #if DISPLAY TCanvas *cBkg2_fit = new TCanvas("cBkg2_fit"); RooPlot *bkg2_frame = x.frame(); hbkg2.plotOn(bkg2_frame); bkg2.plotOn(bkg2_frame, LineColor(kGreen+1)); bkg2_frame->Draw(); cBkg2_fit->Update(); #endif printf("BACKGROUND-2 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-2\n"); printf("mean_bkg2 = %f +- %f\n", mean_bkg2.getVal(), mean_bkg2.getError()); printf("sigma_bkg2 = %f +- %f\n", sigma_bkg2.getVal(), sigma_bkg2.getError()); printf("tail_bkg2 = %f +- %f\n", tail_bkg2.getVal(), tail_bkg2.getError()); /* scale parameters if requested */ if (SCALE_BACKGROUND_SIGMA != 1.) { printf("SCALE FITTED BACKGROUND-2 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA); sigma_bkg2.setVal(sigma_bkg2.getVal() * SCALE_BACKGROUND_SIGMA); } if (SCALE_BACKGROUND_TAIL != 1.) { printf("SCALE FITTED BACKGROUND-2 TAIL BY %f\n", SCALE_BACKGROUND_TAIL); tail_bkg2.setVal(tail_bkg2.getVal() * SCALE_BACKGROUND_TAIL); } /* fix/release parameters if requested */ if (FIX_BACKGROUND_MEAN) { printf("SETTING BACKGROUND-2 FITTED MEAN AS CONSTANTS\n"); mean_bkg2.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-2 FITTED MEAN AS FREE\n"); mean_bkg2.setRange(mean_bkg2.getVal() - 0.25 * TMath::Abs(mean_bkg2.getVal()), mean_bkg2.getVal() + 0.25 * TMath::Abs(mean_bkg2.getVal())); } if (FIX_BACKGROUND_SIGMA) { printf("SETTING BACKGROUND-2 FITTED SIGMA AS CONSTANTS\n"); sigma_bkg2.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-2 FITTED SIGMA AS FREE\n"); sigma_bkg2.setRange(sigma_bkg2.getVal() * 0.75, sigma_bkg2.getVal() * 1.25); } if (FIX_BACKGROUND_TAIL) { printf("SETTING BACKGROUND-2 FITTED TAIL AS CONSTANTS\n"); tail_bkg2.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-2 FITTED TAIL AS FREE\n"); tail_bkg2.setRange(tail_bkg2.getVal() * 0.75, tail_bkg2.getVal() * 1.25); } /* fitting bkg3 */ TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus"); hBkgExp3->Fit(fGaus, "0q"); Double_t bkgexp3_mean = fGaus->GetParameter(1); Double_t bkgexp3_sigma = fGaus->GetParameter(2); RooRealVar mean_bkg3("mean_bkg3", "mean_bkg3", bkgexp3_mean, bkgexp3_mean - 10., bkgexp3_mean + 10., ""); RooRealVar sigma_bkg3("sigma_bkg3", "sigma_bkg3", bkgexp3_sigma, bkgexp3_sigma * 0.5, bkgexp3_sigma * 1.5, ""); RooRealVar tail_bkg3("tail_bkg3", "tail_bkg3", 1., 0.5, 5., ""); RooGaussianTail bkg3("bkg3", "bkg3", x, mean_bkg3, sigma_bkg3, tail_bkg3); bkg3.fitTo(hbkg3, Range(bkgexp3_mean - 5. * bkgexp3_sigma, bkgexp3_mean + 10. * bkgexp3_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10)); #if DISPLAY TCanvas *cBkg3_fit = new TCanvas("cBkg3_fit"); RooPlot *bkg3_frame = x.frame(); hbkg3.plotOn(bkg3_frame); bkg3.plotOn(bkg3_frame, LineColor(kYellow+1)); bkg3_frame->Draw(); cBkg3_fit->Update(); #endif printf("BACKGROUND-3 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-3\n"); printf("mean_bkg3 = %f +- %f\n", mean_bkg3.getVal(), mean_bkg3.getError()); printf("sigma_bkg3 = %f +- %f\n", sigma_bkg3.getVal(), sigma_bkg3.getError()); printf("tail_bkg3 = %f +- %f\n", tail_bkg3.getVal(), tail_bkg3.getError()); /* scale parameters if requested */ if (SCALE_BACKGROUND_SIGMA != 1.) { printf("SCALE FITTED BACKGROUND-3 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA); sigma_bkg3.setVal(sigma_bkg3.getVal() * SCALE_BACKGROUND_SIGMA); } if (SCALE_BACKGROUND_TAIL != 1.) { printf("SCALE FITTED BACKGROUND-3 TAIL BY %f\n", SCALE_BACKGROUND_TAIL); tail_bkg3.setVal(tail_bkg3.getVal() * SCALE_BACKGROUND_TAIL); } /* fix/release parameters if requested */ if (FIX_BACKGROUND_MEAN) { printf("SETTING BACKGROUND-3 FITTED MEAN AS CONSTANTS\n"); mean_bkg3.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-3 FITTED MEAN AS FREE\n"); mean_bkg3.setRange(mean_bkg3.getVal() - 0.25 * TMath::Abs(mean_bkg3.getVal()), mean_bkg3.getVal() + 0.25 * TMath::Abs(mean_bkg3.getVal())); } if (FIX_BACKGROUND_SIGMA) { printf("SETTING BACKGROUND-3 FITTED SIGMA AS CONSTANTS\n"); sigma_bkg3.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-3 FITTED SIGMA AS FREE\n"); sigma_bkg3.setRange(sigma_bkg3.getVal() * 0.75, sigma_bkg3.getVal() * 1.25); } if (FIX_BACKGROUND_TAIL) { printf("SETTING BACKGROUND-3 FITTED TAIL AS CONSTANTS\n"); tail_bkg3.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-3 FITTED TAIL AS FREE\n"); tail_bkg3.setRange(tail_bkg3.getVal() * 0.75, tail_bkg3.getVal() * 1.25); } /* fitting bkg4 */ TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus"); hBkgExp4->Fit(fGaus, "0q"); Double_t bkgexp4_mean = fGaus->GetParameter(1); Double_t bkgexp4_sigma = fGaus->GetParameter(2); RooRealVar mean_bkg4("mean_bkg4", "mean_bkg4", bkgexp4_mean, bkgexp4_mean - 10., bkgexp4_mean + 10., ""); RooRealVar sigma_bkg4("sigma_bkg4", "sigma_bkg4", bkgexp4_sigma, bkgexp4_sigma * 0.5, bkgexp4_sigma * 1.5, ""); RooRealVar tail_bkg4("tail_bkg4", "tail_bkg4", 1., 0.5, 5., ""); RooGaussianTail bkg4("bkg4", "bkg4", x, mean_bkg4, sigma_bkg4, tail_bkg4); bkg4.fitTo(hbkg4, Range(bkgexp4_mean - 5. * bkgexp4_sigma, bkgexp4_mean + 10. * bkgexp4_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10)); #if DISPLAY TCanvas *cBkg4_fit = new TCanvas("cBkg4_fit"); RooPlot *bkg4_frame = x.frame(); hbkg4.plotOn(bkg4_frame); bkg4.plotOn(bkg4_frame, LineColor(kYellow+2)); bkg4_frame->Draw(); cBkg4_fit->Update(); #endif printf("BACKGROUND-4 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-4\n"); printf("mean_bkg4 = %f +- %f\n", mean_bkg4.getVal(), mean_bkg4.getError()); printf("sigma_bkg4 = %f +- %f\n", sigma_bkg4.getVal(), sigma_bkg4.getError()); printf("tail_bkg4 = %f +- %f\n", tail_bkg4.getVal(), tail_bkg4.getError()); /* scale parameters if requested */ if (SCALE_BACKGROUND_SIGMA != 1.) { printf("SCALE FITTED BACKGROUND-4 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA); sigma_bkg4.setVal(sigma_bkg4.getVal() * SCALE_BACKGROUND_SIGMA); } if (SCALE_BACKGROUND_TAIL != 1.) { printf("SCALE FITTED BACKGROUND-4 TAIL BY %f\n", SCALE_BACKGROUND_TAIL); tail_bkg4.setVal(tail_bkg4.getVal() * SCALE_BACKGROUND_TAIL); } /* fix/release parameters if requested */ if (FIX_BACKGROUND_MEAN) { printf("SETTING BACKGROUND-4 FITTED MEAN AS CONSTANTS\n"); mean_bkg4.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-4 FITTED MEAN AS FREE\n"); mean_bkg4.setRange(mean_bkg4.getVal() - 0.25 * TMath::Abs(mean_bkg4.getVal()), mean_bkg4.getVal() + 0.25 * TMath::Abs(mean_bkg4.getVal())); } if (FIX_BACKGROUND_SIGMA) { printf("SETTING BACKGROUND-4 FITTED SIGMA AS CONSTANTS\n"); sigma_bkg4.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-4 FITTED SIGMA AS FREE\n"); sigma_bkg4.setRange(sigma_bkg4.getVal() * 0.75, sigma_bkg4.getVal() * 1.25); } if (FIX_BACKGROUND_TAIL) { printf("SETTING BACKGROUND-4 FITTED TAIL AS CONSTANTS\n"); tail_bkg4.setConstant(kTRUE); } else { printf("SETTING BACKGROUND-4 FITTED TAIL AS FREE\n"); tail_bkg4.setRange(tail_bkg4.getVal() * 0.75, tail_bkg4.getVal() * 1.25); } } else if (GAUSSIAN_BACKGROUND) { printf("USING GAUSSIAN BACKGROUND SHAPES (not reccomended)\n"); RooRealVar mean1("mean1", "mean1", hBkgExp1->GetMean(), hBkgExp1->GetMean() * 0.95, hBkgExp1->GetMean() * 1.05, ""); RooRealVar sigma1("sigma1", "sigma1", hBkgExp1->GetRMS(), hBkgExp1->GetRMS() * 0.5, hBkgExp1->GetRMS() * 1.5, ""); RooGaussian bkg1("bkg1", "bkg1", x, mean1, sigma1); RooRealVar mean2("mean2", "mean2", hBkgExp2->GetMean(), hBkgExp2->GetMean() * 0.95, hBkgExp2->GetMean() * 1.05, ""); RooRealVar sigma2("sigma2", "sigma2", hBkgExp2->GetRMS(), hBkgExp2->GetRMS() * 0.5, hBkgExp2->GetRMS() * 1.5, ""); RooGaussian bkg2("bkg2", "bkg2", x, mean2, sigma2); RooRealVar mean3("mean3", "mean3", hBkgExp3->GetMean(), hBkgExp3->GetMean() * 0.95, hBkgExp3->GetMean() * 1.05, ""); RooRealVar sigma3("sigma3", "sigma3", hBkgExp3->GetRMS(), hBkgExp3->GetRMS() * 0.5, hBkgExp3->GetRMS() * 1.5, ""); RooGaussian bkg3("bkg3", "bkg3", x, mean3, sigma3); } else { printf("SHAPE OF BACKGROUND NOT DEFINE: using EXPECTED_BACKGROUND\n"); RooHistPdf bkg1("bkg1", "bkg1", x, hbkg1, 0); RooHistPdf bkg2("bkg2", "bkg2", x, hbkg2, 0); RooHistPdf bkg3("bkg3", "bkg3", x, hbkg3, 0); RooHistPdf bkg4("bkg4", "bkg4", x, hbkg4, 0); } printf("**************************************************\n"); /*** DEFINE MISMATCH BACKGROUND SHAPE ***/ printf("***** MISMATCH BACKGROUND SHAPE DEFINITION *****\n"); /* variables and generic shapes */ Double_t expectedCutoff = hBkgExp3->GetMean(); Double_t expectedCutoffRMS = hBkgExp3->GetRMS(); // RooRealVar cutoff("cutoff", "cutoff", -30., -50., 0., ""); RooRealVar cutoff("cutoff", "cutoff", expectedCutoff, expectedCutoff - 3. * expectedCutoffRMS, expectedCutoff, ""); // RooRealVar cutoff("cutoff", "cutoff", expectedCutoff, ""); // RooRealVar power("power", "power", 1., 0.5, 5.0, ""); RooRealVar power("power", "power", 1., ""); RooFermiCutoff fermi("fermi", "fermi", x, cutoff, power); RooRealVar alpha("alpha", "alpha", 0., -10., 0., ""); RooExponential expo("expo", "expo", x, alpha); RooUniform uniform("uniform", "uniform", x); /* mismatch shape */ if (EXPECTED_MISMATCH) { printf("USING EXPECTED MISMATCH SHAPE FROM HISTOGRAMS\n"); RooHistPdf mismatch("mismatch", "mismatch", x, hmismatch, 0); } else if (UNIFORM_MISMATCH) { printf("USING UNIFORM BACKGROUND SHAPE WITH CUTOFF\n"); RooProdPdf mismatch("mismatch", "mismatch", fermi, uniform, -100.); } else if (EXPONENTIAL_MISMATCH) { printf("USING EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n"); RooProdPdf mismatch("mismatch", "mismatch", fermi, expo, -100.); } else if (DOUBLEEXPONENTIAL_MISMATCH) { printf("USING DOUBLE-EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n"); RooRealVar alpha1("alpha1", "alpha1", 0., -10., 0., ""); RooRealVar alpha2("alpha2", "alpha2", 0., -10., 0., ""); RooRealVar frac("frac", "frac", 0.5, 0., 1., ""); RooGenericPdf doubleexpo("doubleexpo", "TMath::Exp(alpha1 * x) + frac * TMath::Exp(alpha2 * x)", RooArgSet(x, alpha1, alpha2, frac)); RooProdPdf mismatch("mismatch", "mismatch", fermi, doubleexpo, -100.); } else if (UNIFORMPLUSEXPONENTIAL_MISMATCH) { printf("USING UNIFORM-PLUS-EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n"); RooRealVar funiform("funiform", "funiform", 100., 0., 100000., ""); RooRealVar fexpo("fexpo", "fexpo", 100., 0., 100000., ""); RooAddPdf uniformplusexpo("uniformplusexpo", "uniformplusexpo", RooArgList(uniform, expo), RooArgList(funiform, fexpo), kFALSE); RooProdPdf mismatch("mismatch", "mismatch", fermi, uniformplusexpo, -100.); } else { printf("SHAPE OF MISMATCH NOT DEFINE: using EXPECTED_MISMATCH\n"); RooHistPdf mismatch("mismatch", "mismatch", x, hmismatch, 0); } printf("************************************************\n"); /*** DEFINE THE MODEL ***/ printf("***** MODEL DEFINITION *****\n"); /* variables */ Double_t integral = hdata.sumEntries(); RooRealVar nsignal("nsignal", "nsignal", 0.3 * integral, 0., integral); RooRealVar nbkg1("nbkg1", "nbkg1", 0.3 * integral, 0., integral); RooRealVar nbkg2("nbkg2", "nbkg2", 0.3 * integral, 0., integral); RooRealVar nbkg3("nbkg3", "nbkg3", 0.3 * integral, 0., integral); RooRealVar nbkg4("nbkg4", "nbkg4", 0.3 * integral, 0., integral); RooRealVar nmismatch("nmismatch", "nmismatch", 0.1 * integral, 0., integral); if (!fitBkg1) { nbkg1.setVal(0.); nbkg1.setConstant(kTRUE); } if (!fitBkg2) { nbkg2.setVal(0.); nbkg2.setConstant(kTRUE); } if (!fitBkg3) { nbkg3.setVal(0.); nbkg3.setConstant(kTRUE); } if (!fitBkg4) { nbkg4.setVal(0.); nbkg4.setConstant(kTRUE); } RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3, bkg4, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nbkg3, nbkg4, nmismatch)); #if 0 /* the model */ if (USE_ELECTRON_BACKGROUND && fitElectrons && fitPions) { printf("USING ELECTRON BACKGROUND\n"); if (NO_MISMATCH) { printf("NOT USING MISMATCH BACKGROUND\n"); nmismatch.setVal(0.); RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3/*, bkg4*/), RooArgList(nsignal, nbkg1, nbkg2, nbkg3/*, nbkg4*/)); } else { printf("USING MISMATCH BACKGROUND\n"); RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3/*, bkg4*/, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nbkg3/*, nbkg4*/, nmismatch)); } } else if (!USE_ELECTRON_BACKGROUND || !fitElectrons) { printf("NOT USING ELECTRON BACKGROUND\n"); nbkg3.setVal(0.); nbkg4.setVal(0.); if (NO_MISMATCH) { printf("NOT USING MISMATCH BACKGROUND\n"); nmismatch.setVal(0.); RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2), RooArgList(nsignal, nbkg1, nbkg2)); } else { printf("USING MISMATCH BACKGROUND\n"); RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nmismatch)); } } printf("****************************\n"); #endif /*** FIT ***/ printf("***** FIT *****\n"); printf("SIGNAL PARAMETERS BEFORE GLOBAL FIT\n"); printf("mean = %f +- %f\n", mean.getVal(), mean.getError()); printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError()); printf("tail = %f +- %f\n", tail.getVal(), tail.getError()); /* fit and draw */ if (canvas) canvas->cd(); // model.chi2FitTo(hdata, Extended(kTRUE), Verbose(kFALSE), SumW2Error(kFALSE), Range(-40., 140.), Binned(kTRUE)); // model.fitTo(hdata, Extended(kTRUE), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10), Range(-10., 10.)); model.fitTo(hdata, Range(rangeMin, rangeMax), Extended(kTRUE), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10)); printf("***************\n"); /*** DRAW ***/ #if DISPLAY RooPlot *xframe = x.frame(); hdata.plotOn(xframe, XErrorSize(0), DrawOption("PZ")); printf("plotting model...\n"); model.plotOn(xframe, LineWidth(2)/*, Precision(1.e-4)*/); printf("plotting signal...\n"); model.plotOn(xframe, Components(signal), LineWidth(2), LineColor(kRed)/*, Precision(1.e-4)*/); printf("plotting bkg1...\n"); model.plotOn(xframe, Components(bkg1), LineWidth(2), LineStyle(kDashed), LineColor(kCyan+1)); printf("plotting bkg2...\n"); model.plotOn(xframe, Components(bkg2), LineWidth(2), LineStyle(kDashed), LineColor(kGreen+1)); if (USE_ELECTRON_BACKGROUND) { model.plotOn(xframe, Components(bkg3), LineWidth(2), LineStyle(kDashed), LineColor(kYellow+1)); model.plotOn(xframe, Components(bkg4), LineWidth(2), LineStyle(kDashed), LineColor(kYellow+2)); } if (!NO_MISMATCH) model.plotOn(xframe, Components(mismatch), LineWidth(2), LineStyle(kDashed), LineColor(kGray+1)); hSignal->SetFillColor(kYellow); hSignal->SetLineWidth(0); hSignal->SetFillStyle(0); hSignal->SetMinimum(0.1); hSignal->GetXaxis()->SetRangeUser(rangeMin, rangeMax); // hSignal->Draw(); xframe->SetMinimum(0.1); xframe->Draw(); #endif if (canvas) canvas->Update(); /*** COMPUTE CHI2 ***/ Double_t datax, datapoint, datapoint_err, modelpoint; Double_t chi2 = 0.; Int_t ndf = 0; for (Int_t ibin = 0; ibin < hSignal->GetNbinsX(); ibin++) { datax = hSignal->GetBinCenter(ibin + 1); if (datax < rangeMin || datax > rangeMax) continue; datapoint = hSignal->GetBinContent(ibin + 1); datapoint_err = hSignal->GetBinError(ibin + 1); if (datapoint_err == 0.) continue; x.setVal(datax); modelpoint = model.getVal(); chi2 += (datapoint - modelpoint) * (datapoint - modelpoint) / (datapoint_err * datapoint_err); ndf++; } /*** PRINT FIT OUTPUT ***/ printf("***** CHI-SQUARE *****\n"); printf("chi-square = %f\n", chi2); printf("NDF = %d\n", ndf); printf("chi2/NDF = %f\n", ndf > 0 ? chi2 / ndf : 0.); printf("***** SIGNAL-SHAPE INFO *****\n"); printf("mean = %f +- %f\n", mean.getVal(), mean.getError()); printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError()); printf("tail = %f +- %f\n", tail.getVal(), tail.getError()); printf("pure gaus = %f +- %f\n", gaussianfrac.getVal(), gaussianfrac.getError()); printf("*****************************\n"); printf("***** COUNTS *****\n"); printf("total = %f\n", hSignal->Integral(-1, -1)); printf("integral = %f\n", hdata.sumEntries()); printf("nsignal = %f +- %f\n", nsignal.getVal(), nsignal.getError()); printf("nbkg1 = %f +- %f\n", nbkg1.getVal(), nbkg1.getError()); printf("nbkg2 = %f +- %f\n", nbkg2.getVal(), nbkg2.getError()); printf("nbkg3 = %f +- %f\n", nbkg3.getVal(), nbkg3.getError()); printf("nbkg4 = %f +- %f\n", nbkg4.getVal(), nbkg4.getError()); printf("nmismatch = %f +- %f\n", nmismatch.getVal(), nmismatch.getError()); printf("******************\n"); /*** OUTPUT FIT PARAMS ***/ param[kMean] = mean.getVal(); param_err[kMean] = mean.getError(); param[kSigma] = sigma.getVal(); param_err[kSigma] = sigma.getError(); param[kTail] = tail.getVal(); param_err[kTail] = tail.getError(); param[kTotalCounts] = hSignal->GetEntries(); param_err[kTotalCounts] = sqrt(hSignal->GetEntries()); param[kIntegralCounts] = hdata.sumEntries(); param_err[kIntegralCounts] = sqrt(hdata.sumEntries()); param[kSignalCounts] = nsignal.getVal(); param_err[kSignalCounts] = nsignal.getError(); param[kBkg1Counts] = nbkg1.getVal(); param_err[kBkg1Counts] = nbkg1.getError(); param[kBkg2Counts] = nbkg2.getVal(); param_err[kBkg2Counts] = nbkg2.getError(); param[kBkg3Counts] = nbkg3.getVal(); param_err[kBkg3Counts] = nbkg3.getError(); param[kBkg4Counts] = nbkg4.getVal(); param_err[kBkg4Counts] = nbkg4.getError(); param[kMismatchCounts] = nmismatch.getVal(); param_err[kMismatchCounts] = nmismatch.getError(); return; } //___________________________________________________________________________________ TOFpid_efficiency(Int_t ipart, Int_t icharge, Int_t icent, Int_t useTPCcut = -1, Double_t minsignalFrac = 0., Int_t nTrials = 1000) { /* prepare centrality name */ Char_t centName[1024]; if (icent < 0 || icent >= NcentralityBins) sprintf(centName, "cent0090"); else sprintf(centName, "cent%02d%02d", centralityBin[icent], centralityBin[icent + 1]); /* fraction names */ const Char_t *fractionName[AliPID::kSPECIES][AliPID::kSPECIES] = { "hSignalFraction", "hBkg4Fraction", "hBkg3Fraction", "hBkg1Fraction", "hBkg2Fraction", "hBkg4Fraction", "hSignalFraction", "hBkg3Fraction", "hBkg1Fraction", "hBkg2Fraction", "hBkg3Fraction", "hBkg4Fraction", "hSignalFraction", "hBkg1Fraction", "hBkg2Fraction", "hBkg3Fraction", "hBkg4Fraction", "hBkg1Fraction", "hSignalFraction", "hBkg2Fraction", "hBkg3Fraction", "hBkg4Fraction", "hBkg1Fraction", "hBkg2Fraction", "hSignalFraction" }; enum ETPCcut_t { kCurrentDir, k11Cut, k12Cut, k21Cut, k22Cut, k23Cut, k32Cut, k33Cut, kNTPCcuts }; const Char_t *tpccutdir[8] = { ".", "TOFpid_cutOnTPC[-1.0,1.0]", "TOFpid_cutOnTPC[1.0,2.0]", "TOFpid_cutOnTPC[-2.0,-1.0]", "TOFpid_cutOnTPC[-2.0,2.0]", "TOFpid_cutOnTPC[2.0,3.0]", "TOFpid_cutOnTPC[-3.0,-2.0]", "TOFpid_cutOnTPC[-3.0,3.0]" }; /* get data */ Char_t filename[1024]; TH1D *hAccepted[AliPID::kSPECIES][kNTPCcuts]; TH1D *hIdentified[AliPID::kSPECIES][kNTPCcuts]; TH1D *hEfficiencyIn[AliPID::kSPECIES][kNTPCcuts]; TH1D *hFraction[AliPID::kSPECIES][AliPID::kSPECIES][kNTPCcuts]; for (Int_t itpccut = 0; itpccut < kNTPCcuts; itpccut++) { /* check whether we use this cut */ if (useTPCcut >= 0 && itpccut != useTPCcut) continue; for (Int_t iipart = 0; iipart < AliPID::kSPECIES; iipart++) { /* skip electrons and muons */ if (iipart == AliPID::kMuon) continue; /* open file */ sprintf(filename, "%s/TOFspectrum_%s_%s_%s_%sID.root", tpccutdir[itpccut], centName, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)); TFile *filein = TFile::Open(filename); if (!filein || !filein->IsOpen()) { printf("cannot open %s\n", filename); return; } /* get accepted tracks */ hAccepted[iipart][itpccut] = (TH1D *)filein->Get("PostAnalysis/hAcceptedTracks"); if (!hAccepted[iipart][itpccut]) { printf("cannot find PostAnalysis/hAcceptedTracks in %s\n", filename); return; } /* get identified tracks */ hIdentified[iipart][itpccut] = (TH1D *)filein->Get("PostAnalysis/hIdentifiedCounts"); if (!hIdentified[iipart][itpccut]) { printf("cannot find PostAnalysis/hIdentifiedCounts in %s\n", filename); return; } /* compute efficiency */ hEfficiencyIn[iipart][itpccut] = new TH1D(*hIdentified[iipart][itpccut]); hEfficiencyIn[iipart][itpccut]->Divide(hEfficiencyIn[iipart][itpccut], hAccepted[iipart][itpccut], 1., 1., "B"); /* get fractions */ for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) { /* skip muons */ if (iiipart == AliPID::kMuon) continue; hFraction[iipart][iiipart][itpccut] = (TH1D *)filein->Get(Form("PostAnalysis/%s", fractionName[iipart][iiipart])); if (!hFraction[iipart][iiipart][itpccut]) { printf("cannot find PostAnalysis/%s in %s\n", fractionName[iipart][iiipart], filename); return; } } } } /* prepare output efficiency histos */ TH1D *hEfficiencyOut[AliPID::kSPECIES]; for (Int_t iipart = 0; iipart < AliPID::kSPECIES; iipart++) { hEfficiencyOut[iipart] = new TH1D(Form("hPIDEff_%d_%s_%s_%sID", icent, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NptBins, ptBin); } /*** MATRIX ***/ const Int_t nPart = 4; Int_t partIndex[4] = {AliPID::kElectron, AliPID::kPion, AliPID::kKaon, AliPID::kProton}; TMatrixD Meffin(nPart, 1); TMatrixD Meffout(nPart, 1); TMatrixD Mfrac(nPart, nPart); Double_t eff[4], effe[4], frac[4][4], frace[4][4], geneff, genfrac; Bool_t gotbestcut[4]; TH1F *hEffTemp[4]; for (Int_t iipart = 0; iipart < nPart; iipart++) hEffTemp[iipart] = new TH1F(Form("hEffTemp_%d", iipart), "", 20000, -10., 10.); /* loop over pt-bins */ for (Int_t ibin = 0; ibin < NptBins; ibin++) { /* reset temp histos */ for (Int_t iipart = 0; iipart < nPart; iipart++) hEffTemp[iipart]->Reset(); /* get measured data */ for (Int_t iipart = 0; iipart < nPart; iipart++) { /* select the best set of cuts */ Double_t signalFrac, maxsignalFrac = minsignalFrac; Int_t bestcut = -1; gotbestcut[iipart] = kFALSE; for (Int_t itpccut = 0; itpccut < kNTPCcuts; itpccut++) { /* check whether we use this cut */ if (useTPCcut >= 0 && itpccut != useTPCcut) continue; /* maximize the signal fraction */ signalFrac = hFraction[partIndex[iipart]][partIndex[iipart]][itpccut]->GetBinContent(ibin + 1); if (signalFrac > maxsignalFrac) { maxsignalFrac = signalFrac; bestcut = itpccut; gotbestcut[iipart] = kTRUE; } } if (bestcut < 0) continue; eff[iipart] = hEfficiencyIn[partIndex[iipart]][bestcut]->GetBinContent(ibin + 1); effe[iipart] = hEfficiencyIn[partIndex[iipart]][bestcut]->GetBinError(ibin + 1); for (Int_t iiipart = 0; iiipart < nPart; iiipart++) { frac[iipart][iiipart] = hFraction[partIndex[iipart]][partIndex[iiipart]][bestcut]->GetBinContent(ibin + 1); frace[iipart][iiipart] = hFraction[partIndex[iipart]][partIndex[iiipart]][bestcut]->GetBinError(ibin + 1); } } /* check best cuts */ Bool_t skip = kFALSE; for (Int_t iipart = 0; iipart < nPart; iipart++) if (!gotbestcut[iipart]) skip = kTRUE; if (skip) continue; /* loop over trials */ for (Int_t itry = 0; itry < nTrials; itry++) { /* setup matrix */ for (Int_t iipart = 0; iipart < nPart; iipart++) { geneff = gRandom->Gaus(eff[iipart], effe[iipart]); if (geneff < 0.) geneff = 0.; if (geneff > 1.) geneff = 1.; Meffin[iipart] = geneff != 0. ? 1. / geneff : 0.; for (Int_t iiipart = 0; iiipart < nPart; iiipart++) { genfrac = gRandom->Gaus(frac[iipart][iiipart], frace[iipart][iiipart]); if (genfrac < 0.) genfrac = 0.; if (genfrac > 1.) genfrac = 1.; Mfrac[iipart][iiipart] = genfrac; } } /* invert matrix */ TDecompLU lu(Mfrac); TMatrixD Minv(nPart, nPart); if (!lu.Decompose()) continue; lu.Invert(Minv); Meffout = Minv * Meffin; /* fill histos */ for (Int_t iipart = 0; iipart < nPart; iipart++) { if (Meffout.GetMatrixArray()[iipart] > 0.) hEffTemp[iipart]->Fill(1. / Meffout.GetMatrixArray()[iipart]); } } /* end of loop over trials */ /* get average and RMS */ for (Int_t iipart = 0; iipart < nPart; iipart++) { hEfficiencyOut[partIndex[iipart]]->SetBinContent(ibin + 1, hEffTemp[iipart]->GetMean()); hEfficiencyOut[partIndex[iipart]]->SetBinError(ibin + 1, hEffTemp[iipart]->GetRMS()); } } /* end of loop over pt-bins */ hEfficiencyOut[AliPID::kPion]->SetMarkerStyle(20); hEfficiencyOut[AliPID::kKaon]->SetMarkerStyle(20); hEfficiencyOut[AliPID::kProton]->SetMarkerStyle(20); hEfficiencyOut[AliPID::kElectron]->SetMarkerStyle(20); hEfficiencyOut[AliPID::kPion]->SetMarkerColor(4); hEfficiencyOut[AliPID::kKaon]->SetMarkerColor(8); hEfficiencyOut[AliPID::kProton]->SetMarkerColor(2); hEfficiencyOut[AliPID::kElectron]->SetMarkerColor(1); hEfficiencyOut[AliPID::kPion]->Draw(); hEfficiencyOut[AliPID::kKaon]->Draw("same"); hEfficiencyOut[AliPID::kProton]->Draw("same"); hEfficiencyOut[AliPID::kElectron]->Draw("same"); /* output */ TFile *fileout = TFile::Open(Form("TOFpid_efficiency_cent%d_%s_%s.root", icent, AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE"); hEfficiencyOut[AliPID::kPion]->Write(); hEfficiencyOut[AliPID::kKaon]->Write(); hEfficiencyOut[AliPID::kProton]->Write(); hEfficiencyOut[AliPID::kElectron]->Write(); fileout->Close(); } TOFpid_normalize(TH1D *h, Double_t nevents = 8.42446600000000000e+06) { Double_t counts, counts_err; for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++) { counts = h->GetBinContent(ibin + 1); counts_err = h->GetBinError(ibin + 1); counts /= h->GetBinWidth(ibin + 1); counts_err /= h->GetBinWidth(ibin + 1); h->SetBinContent(ibin + 1, counts); h->SetBinError(ibin + 1, counts_err); } h->Scale(1. / nevents); } TOFpid_normalizeAndwrite(const Char_t *fileoutname, const Char_t *corrstring = "") { TFile *fpiplus = TFile::Open("TOFpid_spectrum_pion_plus.root"); TFile *fpiminus = TFile::Open("TOFpid_spectrum_pion_minus.root"); TFile *fkaplus = TFile::Open("TOFpid_spectrum_kaon_plus.root"); TFile *fkaminus = TFile::Open("TOFpid_spectrum_kaon_minus.root"); TFile *fprplus = TFile::Open("TOFpid_spectrum_proton_plus.root"); TFile *fprminus = TFile::Open("TOFpid_spectrum_proton_minus.root"); hpiplus = (TH1F *)fpiplus->Get(Form("hSignal%sCounts", corrstring)); hpiminus = (TH1F *)fpiminus->Get(Form("hSignal%sCounts", corrstring)); hkaplus = (TH1F *)fkaplus->Get(Form("hSignal%sCounts", corrstring)); hkaminus = (TH1F *)fkaminus->Get(Form("hSignal%sCounts", corrstring)); hprplus = (TH1F *)fprplus->Get(Form("hSignal%sCounts", corrstring)); hprminus = (TH1F *)fprminus->Get(Form("hSignal%sCounts", corrstring)); hpiplus->SetName("hpiplus"); hpiminus->SetName("hpiminus"); hkaplus->SetName("hkaplus"); hkaminus->SetName("hkaminus"); hprplus->SetName("hprplus"); hprminus->SetName("hprminus"); TOFpid_normalize(hpiplus); TOFpid_normalize(hpiminus); TOFpid_normalize(hkaplus); TOFpid_normalize(hkaminus); TOFpid_normalize(hprplus); TOFpid_normalize(hprminus); TFile *fileout = TFile::Open(fileoutname, "RECREATE"); hpiplus->Write("hpiplus"); hpiminus->Write("hpiminus"); hkaplus->Write("hkaplus"); hkaminus->Write("hkaminus"); hprplus->Write("hprplus"); hprminus->Write("hprminus"); fileout->Close(); } /**************************************************************/ TOFpid_rawSpectra(const Char_t *destdir = ".") { Int_t marker[2] = {20, 21}; Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2}; Char_t *partLatex[AliPID::kSPECIES][2] = { "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}" }; Float_t deltaRapidity = rapidityMaxCut - rapidityMinCut; TFile *fileout = TFile::Open("TOF_rawSpectra.root", "RECREATE"); TH1D *hRaw; Char_t title[1024]; for (Int_t icent = -1; icent < NcentralityBins; icent++) for (Int_t icharge = 0; icharge < kNCharges; icharge++) for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) { hRaw = TOFpid_rawSpectra(destdir, ipart, icharge, icent); if (!hRaw) continue; hRaw->Scale(1. / deltaRapidity); hRaw->SetMarkerStyle(marker[icharge]); hRaw->SetMarkerColor(color[ipart]); hRaw->SetLineColor(1); hRaw->SetLineWidth(1); if (icent == -1) sprintf(title, "%s (MB);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge]); else 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]); hRaw->SetTitle(title); if (icent == -1) hRaw->SetName(Form("hRaw_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge])); else hRaw->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge])); fileout->cd(); hRaw->Write(); } fileout->Close(); } TH1D * TOFpid_rawSpectra(const Char_t *dirname, Int_t ipart, Int_t icharge, Int_t icent) { /* open data */ Char_t outfilename[1024]; if (icent < 0 || icent >= NcentralityBins) { sprintf(outfilename, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart)); } else { 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)); } TFile *filein = TFile::Open(outfilename); if (!filein || !filein->IsOpen()) { printf("cannot open %s\n", outfilename); return NULL; } /* get data */ TH1D *h = (TH1D *)filein->Get("RawSpectra/hNormalizedRawYield"); if (!h) { printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename); return NULL; } return h; } TOFpid_rawSpectra_mismatchCorrected(const Char_t *destdir = ".") { Int_t marker[2] = {20, 21}; Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2}; Char_t *partLatex[AliPID::kSPECIES][2] = { "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}" }; Float_t deltaRapidity = rapidityMaxCut - rapidityMinCut; TFile *fileout = TFile::Open("TOF_rawSpectra_mismatchCorrected.root", "RECREATE"); TH1D *hRaw; Char_t title[1024]; for (Int_t icent = -1; icent < NcentralityBins; icent++) for (Int_t icharge = 0; icharge < kNCharges; icharge++) for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) { hRaw = TOFpid_rawSpectra_mismatchCorrected(destdir, ipart, icharge, icent); if (!hRaw) continue; hRaw->Scale(1. / deltaRapidity); hRaw->SetMarkerStyle(marker[icharge]); hRaw->SetMarkerColor(color[ipart]); hRaw->SetLineColor(1); hRaw->SetLineWidth(1); if (icent == -1) sprintf(title, "%s (MB);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge]); else 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]); hRaw->SetTitle(title); if (icent == -1) hRaw->SetName(Form("hRaw_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge])); else hRaw->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge])); fileout->cd(); hRaw->Write(); } fileout->Close(); } TH1D * TOFpid_rawSpectra_mismatchCorrected(const Char_t *dirname, Int_t ipart, Int_t icharge, Int_t icent) { /* open data */ Char_t outfilename[1024]; if (icent < 0 || icent >= NcentralityBins) { sprintf(outfilename, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart)); } else { 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)); } TFile *filein = TFile::Open(outfilename); if (!filein || !filein->IsOpen()) { printf("cannot open %s\n", outfilename); return NULL; } /* get data */ TH1D *h = (TH1D *)filein->Get("RawSpectra/hNormalizedMismatchCorrectedRawYield"); if (!h) { printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename); return NULL; } return h; } /**************************************************************/ TOFpid_electronCorrection() { Int_t marker[2] = {20, 21}; Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2}; Char_t *partLatex[AliPID::kSPECIES][2] = { "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}" }; TFile *fileout = TFile::Open("TOF_electronCorrection.root", "RECREATE"); TH1D *hCorr; TProfile *pCorrAv = new TProfile("pCorrAv", "", NptBins, ptBin, "s"); Char_t title[1024]; for (Int_t icent = 0; icent < NcentralityBins; icent++) for (Int_t icharge = 0; icharge < kNCharges; icharge++) for (Int_t ipart = 2; ipart < 3; ipart++) { hCorr = TOFpid_electronCorrection(ipart, icharge, icent); hCorr->SetMarkerStyle(marker[icharge]); hCorr->SetMarkerColor(color[ipart]); hCorr->SetLineColor(1); hCorr->SetLineWidth(1); 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]); hCorr->SetTitle(title); hCorr->SetName(Form("hElectronCorr_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge])); fileout->cd(); hCorr->Write(); /* fill for average correction */ for (Int_t ipt = 0; ipt < NptBins; ipt++) { pCorrAv->Fill(hCorr->GetBinCenter(ipt + 1), hCorr->GetBinContent(ipt + 1)); } } hCorr = new TH1D("hElectronCorr_average", ";p_{T} (GeV/c);electron-subtracted pion fraction;", NptBins, ptBin); for (Int_t ipt = 0; ipt < NptBins; ipt++) { hCorr->SetBinContent(ipt + 1, pCorrAv->GetBinContent(ipt + 1)); hCorr->SetBinError(ipt + 1, pCorrAv->GetBinError(ipt + 1)); } hCorr->Write(); fileout->Close(); } TH1D * TOFpid_electronCorrection(Int_t ipart, Int_t icharge, Int_t icent) { Int_t marker[2] = {20, 21}; Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2}; Char_t *partLatex[AliPID::kSPECIES][2] = { "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}" }; TH1D *hr = TOFpid_systematics_ratio("defaultFit_electronCut", "defaultFit", ipart, icharge, icent, "electronCorrection", "", marker[icharge], color[ipart], kTRUE); TF1 fOne("fOne", "1.", fitPtMin[ipart], fitPtMax[ipart]); hr->Add(&fOne, -1.); hr->Scale(1. / 0.866); hr->Add(&fOne, 1.); return hr; } /**************************************************************/ SystematicCheck(const Char_t *defaultname = "defaultFit") { const Int_t ndata = 6; const Char_t *name[ndata] = { "signalFit_fixed_scaleSigma_09", "signalFit_fixed_scaleSigma_11", "signalFit_fixed_scaleTail_09", "signalFit_fixed_scaleTail_11", "signalFit_fixed_scaleSigma_09_scaleTail_11", "signalFit_fixed_scaleSigma_11_scaleTail_09" }; for (Int_t idata = 0; idata < ndata; idata++) SystematicCheck(name[idata], defaultname); } SystematicCheck(const Char_t *checkname, const Char_t *defaultname = "defaultFit") { gROOT->LoadMacro("HistoUtils.C"); Char_t filename1[1024]; Char_t filename2[1024]; Int_t marker[2] = {20, 25}; Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2}; Char_t *partLatex[AliPID::kSPECIES][2] = { "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}" }; TFile *fileout = TFile::Open(Form("SystematicCheck_%s.root", checkname), "RECREATE"); Char_t title[1024]; TH1 *hd; for (Int_t icent = 0; icent < NcentralityBins; icent++) for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) for (Int_t icharge = 0; icharge < kNCharges; icharge++) { if (icent < 0 || icent >= NcentralityBins) { sprintf(filename1, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", checkname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart)); sprintf(filename2, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", defaultname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart)); } else { 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)); 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)); } hd = HistoUtils_systematics(filename1, filename2, "RawSpectra/hNormalizedRawYield", NULL); hd->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge])); 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]); hd->SetTitle(title); hd->SetLineColor(1); hd->SetLineWidth(1); hd->SetMarkerStyle(marker[icharge]); hd->SetMarkerColor(color[ipart]); fileout->cd(); hd->Write(); delete hd; } fileout->Close(); } /**************************************************************/ TOFpid_systematics() { TFile *fileout = TFile::Open("TOFpid_systematics.root", "RECREATE"); TH1D *hSys; for (Int_t icent = 0; icent < NcentralityBins; icent++) for (Int_t icharge = 0; icharge < kNCharges; icharge++) for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) { hSys = TOFpid_systematics(ipart, icharge, icent); fileout->cd(); hSys->Write(Form("hRawSys_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge])); delete hSys; } fileout->Close(); } TH1D * TOFpid_systematics(Int_t ipart, Int_t icharge, Int_t icent) { TH1D *hSignalFit = TOFpid_systematics_signalFit(ipart, icharge, icent); TH1D *hBkgFit = TOFpid_systematics_bkgFit(ipart, icharge, icent); TH1D *hFitRange = TOFpid_systematics_fitRange(ipart, icharge, icent); TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin); Double_t sigsys, bkgsys, rangesys, totsys = 0.; for (Int_t ipt = 0; ipt < NptBins; ipt++) { sigsys = hSignalFit->GetBinError(ipt + 1); bkgsys = hBkgFit->GetBinError(ipt + 1); rangesys = hFitRange->GetBinError(ipt + 1); totsys = TMath::Sqrt(sigsys * sigsys + bkgsys * bkgsys + rangesys * rangesys); hSys->SetBinContent(ipt + 1, totsys); } hSys->Draw(); delete hSignalFit; delete hBkgFit; delete hFitRange; return hSys; } TH1D * TOFpid_systematics_fitRange(Int_t ipart, Int_t icharge, Int_t icent) { TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin); hArea->SetMinimum(0.8); hArea->SetMaximum(1.2); hArea->Draw(); TH1D *hr = TOFpid_systematics_ratio("default_wideRange", "default", ipart, icharge, icent, "wide", "wide-range fit;p_{T} (GeV/c);ratio wrt. default;", 20, 4); hr->Draw("same"); TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin); hSys->SetFillStyle(0); hSys->SetMarkerSize(0); for (Int_t ipt = 0; ipt < NptBins; ipt++) { Double_t val, max = 0.; if (hr->GetBinContent(ipt + 1) == 0.) continue; val = TMath::Abs(hr->GetBinContent(ipt + 1) - 1.); hSys->SetBinContent(ipt + 1, 1.); hSys->SetBinError(ipt + 1, val); } hSys->Draw("same, E2"); delete hr; return hSys; } TH1D * TOFpid_systematics_signalFit(Int_t ipart, Int_t icharge, Int_t icent) { const Int_t ndata = 6; const Char_t *name[ndata] = { "signalFit_fixed_scaleSigma_09", "signalFit_fixed_scaleSigma_11", "signalFit_fixed_scaleTail_09", "signalFit_fixed_scaleTail_11", "signalFit_fixed_scaleSigma_09_scaleTail_11", "signalFit_fixed_scaleSigma_11_scaleTail_09" }; const Char_t *title[ndata] = { "-10% #sigma;p_{T} (GeV/c); ratio", "+10% #sigma;p_{T} (GeV/c); ratio", "-10% #tau;p_{T} (GeV/c); ratio", "+10% #tau;p_{T} (GeV/c); ratio", "-10% #sigma, +10% #tau;p_{T} (GeV/c); ratio", "+10% #sigma, -10% #tau;p_{T} (GeV/c); ratio" }; Int_t marker[ndata] = {22, 28, 22, 28, 22, 28}; Int_t color[ndata] = {2, 2, 8, 8, 4, 4}; TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin); hArea->SetMinimum(0.8); hArea->SetMaximum(1.2); hArea->Draw(); TH1D *hr[ndata]; for (Int_t idata = 0; idata < ndata; idata++) { hr[idata] = TOFpid_systematics_ratio(name[idata], "defaultFit", ipart, icharge, icent, name[idata], title[idata], marker[idata], color[idata]); hr[idata]->Draw("same"); } TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin); hSys->SetFillStyle(0); hSys->SetMarkerSize(0); for (Int_t ipt = 0; ipt < NptBins; ipt++) { Double_t val, max = 0.; for (Int_t idata = 0; idata < ndata; idata++) { if (hr[idata]->GetBinContent(ipt + 1) == 0.) continue; val = TMath::Abs(hr[idata]->GetBinContent(ipt + 1) - 1.); if (val > 0.9) continue; if (val > max) max = val; } hSys->SetBinContent(ipt + 1, 1.); hSys->SetBinError(ipt + 1, max); } hSys->Draw("same, E2"); // delete hArea; // for (Int_t idata = 0; idata < ndata; idata++) // delete hr[idata]; return hSys; } TH1D * TOFpid_systematics_bkgFit(Int_t ipart, Int_t icharge, Int_t icent) { const Int_t ndata = 6; const Char_t *name[ndata] = { "bkgFit_fixed_scaleSigma_09", "bkgFit_fixed_scaleSigma_11", "bkgFit_fixed_scaleTail_09", "bkgFit_fixed_scaleTail_11", "bkgFit_fixed_scaleSigma_09_scaleTail_11", "bkgFit_fixed_scaleSigma_11_scaleTail_09" }; const Char_t *title[ndata] = { "-10% #sigma;p_{T} (GeV/c); ratio", "+10% #sigma;p_{T} (GeV/c); ratio", "-10% #tau;p_{T} (GeV/c); ratio", "+10% #tau;p_{T} (GeV/c); ratio", "-10% #sigma, +10% #tau;p_{T} (GeV/c); ratio", "+10% #sigma, -10% #tau;p_{T} (GeV/c); ratio" }; Int_t marker[ndata] = {22, 28, 22, 28, 22, 28}; Int_t color[ndata] = {2, 2, 8, 8, 4, 4}; TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin); hArea->SetMinimum(0.5); hArea->SetMaximum(1.5); hArea->Draw(); TH1D *hr[ndata]; for (Int_t idata = 0; idata < ndata; idata++) { hr[idata] = TOFpid_systematics_ratio(name[idata], "defaultFit", ipart, icharge, icent, name[idata], title[idata], marker[idata], color[idata]); hr[idata]->Draw("same"); } TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin); hSys->SetFillStyle(0); hSys->SetMarkerSize(0); for (Int_t ipt = 0; ipt < NptBins; ipt++) { Double_t val, max = 0.; for (Int_t idata = 0; idata < ndata; idata++) { if (hr[idata]->GetBinContent(ipt + 1) == 0.) continue; val = TMath::Abs(hr[idata]->GetBinContent(ipt + 1) - 1.); if (val > 0.9) continue; if (val > max) max = val; } hSys->SetBinContent(ipt + 1, 1.); hSys->SetBinError(ipt + 1, max); } hSys->Draw("same, E2"); // delete hArea; // for (Int_t idata = 0; idata < ndata; idata++) // delete hr[idata]; return hSys; } TH1D * 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) { TH1D *hr = new TH1D("hr", "", NptBins, ptBin); /* open data */ Char_t outfilename1[1024]; Char_t outfilename2[1024]; if (icent < 0 || icent >= NcentralityBins) { sprintf(outfilename1, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname1, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart)); sprintf(outfilename2, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname2, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart)); } else { 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)); 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)); } TFile *filein1 = TFile::Open(outfilename1); if (!filein1 || !filein1->IsOpen()) { printf("cannot open %s\n", outfilename1); return; } TFile *filein2 = TFile::Open(outfilename2); if (!filein2 || !filein2->IsOpen()) { printf("cannot open %s\n", outfilename2); return; } /* get data */ TH1D *h1 = (TH1D *)filein1->Get("RawSpectra/hNormalizedRawYield"); if (!h1) { printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename1); return; } TH1D *h2 = (TH1D *)filein2->Get("RawSpectra/hNormalizedRawYield"); if (!h2) { printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename2); return; } /* ratio */ if (correlated) hr->Divide(h1, h2, 1., 1., "B"); else hr->Divide(h1, h2); hr->SetNameTitle(name, title); hr->SetMarkerStyle(marker); hr->SetMarkerColor(color); filein1->Close(); filein2->Close(); return hr; }