00001 #include "NCUtils/Extrapolation/PIDSpectrum.h"
00002
00003 #include "NCUtils/Extrapolation/NCBeam.h"
00004
00005 #include "TCanvas.h"
00006 #include "TLatex.h"
00007 #include "TMath.h"
00008 #include "TPad.h"
00009 #include "TStyle.h"
00010 #include "TLegend.h"
00011
00012 #include <cmath>
00013 #include <iostream>
00014 #include <algorithm>
00015
00016 using namespace std;
00017
00018 ClassImp(PIDSpectrum)
00019
00020 PIDSpectrum::PIDSpectrum(string _name, string _title,
00021 int _nPIDBins, double _PIDmin, double _PIDmax,
00022 int _trueBinFactor, NCBeam::Info _beamInfo,
00023 double _cutValue)
00024 : TNamed(_name, _title),
00025 fNumPIDBins(_nPIDBins),
00026 fPIDmin(_PIDmin), fPIDmax(_PIDmax),
00027 fTrueBinFactor(_trueBinFactor),
00028 fBeamInfo(_beamInfo)
00029 {
00030 fEvtTypes.push_back(kNC);
00031 fEvtTypes.push_back(kNumu);
00032 fEvtTypes.push_back(kBeamNue);
00033 fEvtTypes.push_back(kOscNue);
00034 fEvtTypes.push_back(kNutau);
00035
00036 fDrawColors.resize(fEvtTypes.size());
00037 fDrawColors[kNC]=kRed;
00038 fDrawColors[kNumu]=kBlue;
00039 fDrawColors[kBeamNue]=kGreen+2;
00040 fDrawColors[kOscNue]=kOrange-3;
00041 fDrawColors[kNutau]=kMagenta+2;
00042
00043
00044 Double_t* trueEBins=MultiplyBins(fTrueBinFactor, kNumEnergyBinsFar,
00045 kEnergyBinsFar);
00046 int nTrueEBins=kNumEnergyBinsFar*fTrueBinFactor;
00047
00048
00049
00050 if (fNumPIDBins==-1) {
00051
00052 fPIDBins=new Double_t[3];
00053 fPIDBins[0]=fPIDmin; fPIDBins[2]=fPIDmax;
00054 fPIDBins[1]=_cutValue;
00055
00056
00057 fNumPIDBins=2;
00058 } else {
00059 fPIDBins=new Double_t[fNumPIDBins+1];
00060 for (int i=0; i<fNumPIDBins+1; ++i)
00061 fPIDBins[i]=fPIDmin+float(i)*(fPIDmax-fPIDmin)/fNumPIDBins;
00062 }
00063
00064
00065 fNumBinEdgesDontUseThis=fNumPIDBins+1;
00066
00067 for (vector<evtType>::iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
00068 string tmp("hTrue");
00069 tmp+=GetName();
00070 tmp+=EvtTypeAsString(*it);
00071
00072 fCmpTrue[*it]=new TH2F(tmp.c_str(), "True energy spectrum;True Energy (GeV);PID",
00073 nTrueEBins, trueEBins,
00074 fNumPIDBins, fPIDBins);
00075 tmp="t2r_";
00076 tmp+=GetName();
00077 tmp+=EvtTypeAsString(*it);
00078
00079 fTrueToReco[*it]=new TH3F(tmp.c_str(), "True-to-reco;Reco Energy;PID;True energy",
00080 kNumEnergyBinsFar, kEnergyBinsFar,
00081 fNumPIDBins, fPIDBins,
00082 nTrueEBins, trueEBins);
00083 }
00084
00085 delete[] trueEBins;
00086
00087 string tmp("hData");
00088 tmp+=GetName();
00089
00090 fData=new TH2F(tmp.c_str(), "Far Data;Reco Energy;PID",
00091 kNumEnergyBinsFar, kEnergyBinsFar,
00092 fNumPIDBins, fPIDBins);
00093
00094 tmp="hNearMC";
00095 tmp+=GetName();
00096 fNearMC=new TH2F(tmp.c_str(), "Near MC;Reco Energy;PID",
00097 kNumEnergyBinsFar, kEnergyBinsFar,
00098 fNumPIDBins, fPIDBins);
00099
00100 tmp="hNearData";
00101 tmp+=GetName();
00102 fNearData=new TH2F(tmp.c_str(), "Near Data;Reco Energy;PID",
00103 kNumEnergyBinsFar, kEnergyBinsFar,
00104 fNumPIDBins, fPIDBins);
00105
00106 }
00107
00108
00109
00110 Double_t* PIDSpectrum::MultiplyBins(int factor, int nBinsOrig,
00111 const Double_t* origBins) const
00112 {
00113 Double_t* newBins=new Double_t[(nBinsOrig+1)*factor+2];
00114
00115 for (int iOrig=0; iOrig<nBinsOrig+1; ++iOrig) {
00116 for (int j=0; j<factor; ++j) {
00117 newBins[iOrig*factor+j]=origBins[iOrig]+
00118 float(j)/factor*(origBins[iOrig+1]-origBins[iOrig]);
00119 }
00120 }
00121
00122 return newBins;
00123 }
00124
00125
00126
00127 TH2F* PIDSpectrum::GetPredicted(const NC::OscProb::OscPars* coords,
00128 const char* histname,
00129 bool doFN,
00130 TH1* nearSpectrum) const
00131 {
00132 TH2F* pred=new TH2F(histname, "Predicted total",
00133 kNumEnergyBinsFar, kEnergyBinsFar,
00134 fNumPIDBins, fPIDBins);
00135 pred->SetDirectory(0);
00136
00137 for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
00138 AddOnePredicted(pred, *it, coords);
00139 }
00140
00141
00142
00143
00144
00145 pred->SetEntries(int(pred->Integral()));
00146
00147
00148 if (doFN)
00149 FarNearCorrect(pred, fNearData, nearSpectrum ? nearSpectrum : fNearMC);
00150
00151 return pred;
00152
00153 }
00154
00155
00156
00157 TH2F* PIDSpectrum::GetOnePredicted(evtType t,
00158 const NC::OscProb::OscPars* coords,
00159 const char* histname) const
00160 {
00161
00162 string tmp = GetName();
00163 tmp+="_pred";
00164 tmp+=histname ? histname : "";
00165 tmp+=EvtTypeAsString(t);
00166 TH2F* onePred=new TH2F(tmp.c_str(),
00167 Form("Predicted %s", EvtTypeAsString(t).c_str()),
00168 kNumEnergyBinsFar, kEnergyBinsFar,
00169 fNumPIDBins, fPIDBins);
00170
00171 AddOnePredicted(onePred, t, coords);
00172
00173 return onePred;
00174 }
00175
00176
00177
00178 void PIDSpectrum::AddOnePredicted(TH2F* onePred, evtType t, const NC::OscProb::OscPars* coords) const
00179 {
00180
00181 Double_t* trueEBins=MultiplyBins(fTrueBinFactor, kNumEnergyBinsFar,
00182 kEnergyBinsFar);
00183 const int nTrueEBins=kNumEnergyBinsFar*fTrueBinFactor;
00184
00185 TH1F* oscProb=GetOscWeightHist(t, coords, nTrueEBins, trueEBins);
00186 delete[] trueEBins;
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207 TH3F* t2r=fTrueToReco.find(t)->second;
00208 TH2F* trueCmp=fCmpTrue.find(t)->second;
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252 float* pt2r = t2r->fArray;
00253
00254 for(int iTrueE = 0; iTrueE <= nTrueEBins+1; ++iTrueE){
00255 const float oscWeight = oscProb->fArray[iTrueE];
00256
00257
00258 float* pOnePred = onePred->fArray;
00259
00260 for(int iPID=0; iPID <= fNumPIDBins+1; ++iPID){
00261
00262
00263
00264 const int trueE_PIDBin = iTrueE+(nTrueEBins+2)*iPID;
00265
00266 const float pidWeight = trueCmp->fArray[trueE_PIDBin];
00267 const float oscPidWeight = oscWeight * pidWeight;
00268
00269 for(int iRecoE = 0; iRecoE <= kNumEnergyBinsFar+1; ++iRecoE){
00270
00271 const float fillValue = (*pt2r)*oscPidWeight;
00272
00273 *pOnePred += fillValue;
00274
00275
00276 ++pt2r;
00277 ++pOnePred;
00278
00279 }
00280 }
00281 }
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291 }
00292
00293
00294
00295 double PIDSpectrum::GetEventWeight(evtType t, double trueE,
00296 const NC::OscProb::OscPars* coords) const
00297 {
00298 switch (t) {
00299 case PIDSpectrum::kNC:
00300
00301
00302
00303
00304
00305
00306 return coords->TransitionProbability(NCType::kNuMuToNuMu,
00307 NCType::kNC,
00308 NCType::kBaseLineFar,
00309 trueE);
00310 break;
00311 case PIDSpectrum::kNumu:
00312 return coords->TransitionProbability(NCType::kNuMuToNuMu,
00313 NCType::kCC,
00314 NCType::kBaseLineFar,
00315 trueE);
00316 break;
00317 case PIDSpectrum::kBeamNue:
00318
00319
00320
00321
00322
00323 return 1.;
00324 break;
00325 case PIDSpectrum::kOscNue:
00326 return coords->TransitionProbability(NCType::kNuMuToNuE,
00327 NCType::kCC,
00328 NCType::kBaseLineFar,
00329 trueE);
00330 break;
00331 case PIDSpectrum::kNutau:
00332 return coords->TransitionProbability(NCType::kNuMuToNuTau,
00333 NCType::kCC,
00334 NCType::kBaseLineFar,
00335 trueE);
00336 break;
00337 default:
00338
00339 cerr << "GetOnePredicted got unknown evtType " << t << endl;
00340 assert(0);
00341 }
00342 }
00343
00344
00345
00346 TH1F* PIDSpectrum::GetOscWeightHist(evtType t,
00347 const NC::OscProb::OscPars* coords,
00348 int nBins, Double_t* binEdges) const
00349 {
00350
00351 static TH1F* oscProb=new TH1F("oscProb", "Oscillation probability", nBins, binEdges);
00352
00353
00354
00355
00356
00357
00358 for (int i=1; i<=oscProb->GetNbinsX()+1; ++i) {
00359
00360 double Emin=oscProb->GetXaxis()->GetBinLowEdge(i);
00361 double Emax=oscProb->GetXaxis()->GetBinUpEdge(i);
00362 double trueE=0.5*(Emin+Emax);
00363 oscProb->fArray[i] = GetEventWeight(t, trueE, coords);
00364 }
00365
00366 return oscProb;
00367 }
00368
00369
00370
00371 double PIDSpectrum::PIDWithinRange(double origPID) const
00372 {
00373 if(origPID>fPIDmax) return fPIDmax-1e-3;
00374 if(origPID<fPIDmin) return fPIDmin+1e-3;
00375 return origPID;
00376 }
00377
00378
00379
00380 void PIDSpectrum::FillMC(Detector::Detector_t det,
00381 const ANtpTruthInfoBeam* t,
00382 const ANtpRecoInfo* ,
00383 const ANtpAnalysisInfo* ana,
00384 double E, double scale)
00385 {
00386 evtType evt(GetEvtType(t));
00387
00388 double pidVal=PIDWithinRange(ana->separationParameter);
00389
00390
00391
00392
00393 E=std::min(E,kEnergyBinsFar[kNumEnergyBinsFar]-0.01);
00394
00395 if (det==Detector::kFar) {
00396 fCmpTrue[evt]->Fill(t->nuEnergy, pidVal, scale);
00397 fTrueToReco[evt]->Fill(E, pidVal, t->nuEnergy, scale);
00398 } else if (det==Detector::kNear) {
00399 fNearMC->Fill(E, pidVal, scale);
00400 } else {
00401 cout << "Unknown detector " << det << ". Bailing" << endl;
00402 exit(1);
00403 }
00404 }
00405
00406
00407
00408 void PIDSpectrum::FillData(Detector::Detector_t det,
00409 const ANtpTruthInfoBeam* ,
00410 const ANtpRecoInfo* ,
00411 const ANtpAnalysisInfo* ana,
00412 double E, double scale)
00413 {
00414
00415
00416
00417
00418 E=std::min(E,kEnergyBinsFar[kNumEnergyBinsFar]-0.01);
00419
00420 double pidVal=PIDWithinRange(ana->separationParameter);
00421 if (det==Detector::kFar) fData->Fill(E, pidVal, scale);
00422 else if (det==Detector::kNear) fNearData->Fill(E, pidVal, scale);
00423 else {
00424 cout << "Unknown detector " << det << ". Bailing" << endl;
00425 exit(1);
00426 }
00427 }
00428
00429
00430
00431 void PIDSpectrum::NormalizeTrueToRecos()
00432 {
00433 for (vector<evtType>::iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
00434 TH3F* t2r = fTrueToReco[*it];
00435 const int maxE = t2r->GetNbinsZ()+1;
00436 TH2F* ct = fCmpTrue[*it];
00437 TH3F* denom=(TH3F*)t2r->Clone("denom");
00438 denom->Reset();
00439 denom->SetDirectory(0);
00440 for (int iRecoE=1; iRecoE<kNumEnergyBinsFar+1; ++iRecoE) {
00441 for (int iPID=1; iPID<fNumPIDBins+1; ++iPID) {
00442
00443 for (int iTrueE=1; iTrueE<=maxE; ++iTrueE) {
00444 denom->SetBinContent(iRecoE, iPID, iTrueE,
00445 ct->GetBinContent(iTrueE, iPID));
00446 }
00447 }
00448 }
00449 t2r->Divide(denom);
00450 delete denom;
00451 }
00452 }
00453
00454
00455
00456 PIDSpectrum::evtType PIDSpectrum::GetEvtType(const ANtpTruthInfoBeam* t) const
00457 {
00458 if (t->interactionType==0) return kNC;
00459 else if (abs(t->nuFlavor)==14) return kNumu;
00460 else if (abs(t->nuFlavor)==12 && abs(t->nonOscNuFlavor)==12) return kBeamNue;
00461 else if (abs(t->nuFlavor)==12 && abs(t->nonOscNuFlavor)!=12) return kOscNue;
00462 else if (abs(t->nuFlavor)==16) return kNutau;
00463 else {
00464 cerr << "Unknown particle type " << t->nuFlavor << endl;
00465 assert(0);
00466 }
00467 }
00468
00469
00470
00471 string PIDSpectrum::EvtTypeAsString(evtType t) const
00472 {
00473
00474 switch (t) {
00475 case kNC:
00476 return "NC";
00477 case kNumu:
00478 return "NuMu";
00479 case kBeamNue:
00480 return "BeamNuE";
00481 case kOscNue:
00482 return "OscNuE";
00483 case kNutau:
00484 return "NuTau";
00485 default:
00486 cerr << "EvtTypeAsString got unknown evtType " << t << endl;
00487 return "Unknown";
00488 }
00489 }
00490
00491
00492
00493 TCanvas* PIDSpectrum::DrawSpectrum(const NC::OscProb::OscPars* coords,
00494 TH1* systratios,
00495 string suffix) const
00496 {
00497 TH2F* hPred=GetPredicted(coords, "total", false);
00498
00499 if (systratios) hPred->Multiply(systratios);
00500 string canName("can_");
00501 canName+=GetName();
00502 canName+=suffix;
00503 TCanvas *c = new TCanvas(canName.c_str(), GetTitle());
00504 c->Divide(2,2);
00505
00506 c->cd(1);
00507 hPred->DrawCopy("colz");
00508 TLatex* la=new TLatex;
00509 la->SetNDC();
00510 la->DrawLatex(0.2, 0.95, "Predicted");
00511
00512 c->cd(2);
00513 fData->DrawCopy("colz");
00514 la->DrawLatex(0.2, 0.95, "Data");
00515
00516 c->cd(3);
00517 TH2F* hRatio=(TH2F*)fData->Clone("ratio");
00518 hRatio->Divide(hPred);
00519 hRatio->Draw("colz");
00520 la->DrawLatex(0.2, 0.95, "Ratio");
00521
00522 c->cd(4);
00523 fCmpTrue.find(kNC)->second->DrawCopy("colz");
00524 la->DrawLatex(0.2, 0.95, "True NC true E");
00525
00526
00527
00528 return c;
00529 }
00530
00531
00532
00533 vector<TCanvas*> PIDSpectrum::DrawEnergySlices(const NC::OscProb::OscPars* coords,
00534 TH1* systratios,
00535 string suffix) const
00536 {
00537 return DrawSlices("energy", coords, systratios, suffix);
00538 }
00539
00540
00541
00542 vector<TCanvas*> PIDSpectrum::DrawPIDSlices(const NC::OscProb::OscPars* coords,
00543 TH1* systratios,
00544 string suffix) const
00545 {
00546 return DrawSlices("pid", coords, systratios, suffix);
00547 }
00548
00549
00550
00551 void PIDSpectrum::ResetMC()
00552 {
00553 typedef map<evtType, TH2F*>::iterator it2;
00554 for(it2 it = fCmpTrue.begin(); it != fCmpTrue.end(); ++it)
00555 it->second->Reset();
00556
00557 typedef map<evtType, TH3F*>::iterator it3;
00558 for(it3 it = fTrueToReco.begin(); it != fTrueToReco.end(); ++it)
00559 it->second->Reset();
00560 }
00561
00562
00563
00564 vector<TCanvas*> PIDSpectrum::DrawSlices(string variable,
00565 const NC::OscProb::OscPars* coords,
00566 TH1* systratios,
00567 string suffix) const
00568 {
00569
00570 assert(variable == "energy" || variable == "pid");
00571 bool energySlice=(variable=="energy");
00572
00573 vector<TCanvas*> v;
00574
00575
00576
00577
00578
00579 static int canvN=0;
00580 TCanvas* cTotal = new TCanvas((variable+"Total"+suffix+Form("%d", canvN)).c_str());
00581
00582 TCanvas* cSlices = new TCanvas((variable+"Slices"+suffix+Form("%d", canvN)).c_str());
00583 ++canvN;
00584
00585 v.push_back(cTotal); v.push_back(cSlices);
00586
00587 cSlices->cd();
00588 int nPlots = energySlice ? fNumPIDBins : kNumEnergyBinsFar;
00589 int nDiv = (int)TMath::Ceil(TMath::Sqrt(nPlots));
00590
00591
00592 if ((nDiv-1)*nDiv>=nPlots)
00593 cSlices->Divide(nDiv-1, nDiv);
00594 else if ((nDiv+1)*(nDiv-1)>=nPlots)
00595 cSlices->Divide(nDiv-1, nDiv+1);
00596 else
00597 cSlices->Divide(nDiv, nDiv);
00598
00599 TH2F* hTotal=GetPredicted(coords, "total", false);
00600 if (systratios) hTotal->Multiply(systratios);
00601
00602
00603 hTotal->SetTitle("Best fit");
00604 NC::OscProb::NoOscillations oscParsNoOsc;
00605 TH2F* hNoOsc=GetPredicted(&oscParsNoOsc, "noosc", false);
00606 hNoOsc->SetTitle("No Oscillations");
00607
00608 map<evtType, TH2F*> cmpHists;
00609
00610 for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
00611
00612 if (*it==kNC) continue;
00613 cmpHists[*it]=GetOnePredicted(*it, coords, (variable+"Slices"+suffix).c_str());
00614 cmpHists[*it]->SetLineColor(fDrawColors[*it]);
00615 }
00616
00617 hTotal->SetLineColor(kRed);
00618 hNoOsc->SetLineColor(kRed); hNoOsc->SetLineStyle(2);
00619 fData->SetLineColor(kBlack);
00620
00621
00622 gStyle->SetOptStat(0);
00623
00624
00625 cTotal->cd();
00626 TH1F* noOsc1D=(TH1F*)ProjectHist(hNoOsc, variable, "NoOsc")->DrawCopy();
00627 noOsc1D->SetLineStyle(2);
00628 if (energySlice) noOsc1D->GetXaxis()->SetRangeUser(0,19.9);
00629 ProjectHist(hTotal, variable, "Total")->DrawCopy("same");
00630
00631 for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
00632
00633 if (*it==kNC) continue;
00634 ProjectHist(cmpHists[*it], variable,
00635 string("slice")+EvtTypeAsString(*it))->DrawCopy("same");
00636 }
00637 ProjectHist(fData, variable, "Data")->DrawCopy("same e");
00638
00639 ((TPad*)gPad)->BuildLegend()->SetFillColor(kWhite);
00640 if (!energySlice && noOsc1D->GetEntries()) ((TPad*)gPad)->SetLogy();
00641
00642
00643 for (int i=1; i<nPlots+1; ++i) {
00644 cSlices->cd(i);
00645
00646 TH1F* noOsc1DSlice=(TH1F*)ProjectHist(hNoOsc, variable,
00647 "NoOsc", i, i)->DrawCopy();
00648 noOsc1DSlice->SetLineStyle(2);
00649 if (energySlice) noOsc1DSlice->GetXaxis()->SetRangeUser(0,19.9);
00650
00651 ProjectHist(hTotal, variable, "Total", i, i)->DrawCopy("same");
00652
00653 for(unsigned int n = 0; n < fEvtTypes.size(); ++n){
00654
00655 if (fEvtTypes[n]==kNC) continue;
00656
00657 ProjectHist(cmpHists[fEvtTypes[n]], variable,
00658 EvtTypeAsString(fEvtTypes[n]),
00659 i, i)->DrawCopy("same");
00660 }
00661 ProjectHist(fData, variable, "Data", i, i)->DrawCopy("same e");
00662
00663 if (!energySlice && noOsc1DSlice->GetEntries()) ((TPad*)gPad)->SetLogy();
00664 }
00665
00666 return v;
00667 }
00668
00669
00670
00671 TH1D* PIDSpectrum::ProjectHist(TH2* h, string axis, string histname,
00672 int lobin, int hibin) const
00673 {
00674 TH1D* ret;
00675 histname=axis+histname;
00676 if(lobin) histname+=TString::Format("%d", lobin);
00677 else histname+="Integ";
00678
00679 if (axis=="energy") ret = h->ProjectionX(histname.c_str(), lobin, hibin);
00680 else if (axis=="pid") ret = h->ProjectionY(histname.c_str(), lobin, hibin);
00681 else assert(0 && "unknown projection requested");
00682 ret->SetDirectory(0);
00683 TH1D* clone=(TH1D*)ret->Clone();
00684 clone->SetDirectory(0);
00685 return clone;
00686 }
00687
00688
00689
00690 vector<TCanvas*> PIDSpectrum::DrawNearSpectra(TH1* systratios, string suffix) const
00691 {
00692 vector<TCanvas*> ret;
00693 TString canvName("nearSpectra "+fBeamInfo.GetDescription());
00694 canvName+=suffix.c_str();
00695 TCanvas* c = new TCanvas(canvName, "Near detector spectra");
00696 c->Divide(2,2);
00697 c->cd(1);
00698 TH2F* nearMCShifted=(TH2F*)fNearMC->Clone("nearMCShifted");
00699 if (systratios) nearMCShifted->Multiply(systratios);
00700 nearMCShifted->Draw("colz");
00701 c->cd(2);
00702 fNearData->Draw("colz");
00703 c->cd(3);
00704 TH2F* hRatio=(TH2F*)fNearData->Clone("hRatio");
00705 hRatio->Divide(nearMCShifted);
00706 hRatio->Draw("colz");
00707 ret.push_back(c);
00708
00709 nearMCShifted->SetLineColor(kRed);
00710 fNearData->SetLineColor(kBlack);
00711
00712 TCanvas* cE=new TCanvas(canvName+"energy_1D", "Near detector 1D spectra");
00713 cE->cd();
00714 ProjectHist(nearMCShifted, "energy", "nearMC")->DrawCopy();
00715 ProjectHist(fNearData, "energy", "nearData")->DrawCopy("same e");
00716 ret.push_back(cE);
00717
00718 TCanvas* cPID=new TCanvas(canvName+"pid_1D", "Near detector 1D spectra");
00719 cPID->cd();
00720 ProjectHist(nearMCShifted, "pid", "nearMC")->DrawCopy();
00721 ProjectHist(fNearData, "pid", "nearData")->DrawCopy("same e");
00722 ret.push_back(cPID);
00723
00724 return ret;
00725 }
00726
00727
00728
00729 void PIDSpectrum::FarNearCorrect(TH2F* fdPred,
00730 const TH2F* ndData,
00731 const TH1* ndMC1D) const
00732 {
00733
00734
00735
00736
00737
00738 TH2F* ndMC = (TH2F*)ndMC1D;
00739
00740
00741
00742 for (int i=0; i<fdPred->fN; ++i) {
00743
00744
00745 if (ndMC->fArray[i]) fdPred->fArray[i]*=ndData->fArray[i]/ndMC->fArray[i];
00746 }
00747 }