PairAnalysisSpectrum.cxx 28.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
///////////////////////////////////////////////////////////////////////////
//                PairAnalysisSpectrum                                   //
//                                                                       //
//                                                                       //
/*

  post-processing class to extract signals, efficiency, apply corrections,
  calculate systematics and describe the spectra by some functions/models.

  Example:

  PairAnalysisSpectrum *spectrum = new PairAnalysisSpectrum("legend header","systematics");

  spectrum->SetParticleOfInterest(pdgcode);
  spectrum->SetVariable("Pt",PairAnalysisHelper::MakeLinBinning(20,0.,2.) );

  spectrum->SetSystMethod( PairAnalysisSpectrum::kSystMax );

  // add input spectra coming from PairAnalysisHistos
  spectrum->AddInput( histos->DrawSame("pM_Pt","nomc goff"),           // raw invariant mass spectrum
                      "like-sign",                                     // unique string
		      histos->DrawSame("pPt","onlymc goff sel","phi")  // optional: MC spectra for efficiency
		     );
  spectrum->AddExtractor( sig );    // signal extraction (see PairAnalysisSignalExt)

  .... add more input as much as you want


  // process all inputs
  spectrum->Process();


  // draw spectra using TTree::Draw command + some extra arguments (see Draw)
  // see Extraction for content
  spectrum->Draw("s/eff:var","","leg logY syst P");

  // fit the spectrum by some predefined function (see PairAnalysisFunction)
  spectrum->SetFitRange(0.,2.);
  spectrum->SetDefault( PairAnalysisFunction::kBoltzmann );
  spectrum->SetFitOption("RN0");
  spectrum->Fit("leg L");

 */
//                                                                       //
///////////////////////////////////////////////////////////////////////////

//#include <TObject.h>
Etienne Bechtel's avatar
Etienne Bechtel committed
48
#include "PairAnalysisSpectrum.h"
49
50

#include <TCanvas.h>
Administrator's avatar
Administrator committed
51
#include <TDatabasePDG.h>
Etienne Bechtel's avatar
Etienne Bechtel committed
52
#include <TEventList.h>
53
#include <TF1.h>
Administrator's avatar
Administrator committed
54
55
56
#include <TFile.h>
#include <TFormula.h>
#include <TGraphErrors.h>
57
#include <TH1.h>
Administrator's avatar
Administrator committed
58
#include <TList.h>
59
#include <TObjString.h>
Etienne Bechtel's avatar
Etienne Bechtel committed
60
#include <TPad.h>
Administrator's avatar
Administrator committed
61
#include <TPaveText.h>
62
63
#include <TProfile.h>
#include <TProfile2D.h>
Etienne Bechtel's avatar
Etienne Bechtel committed
64
65
#include <TROOT.h>
#include <TTree.h>
Administrator's avatar
Administrator committed
66
#include <TVectorT.h>
67
68

#include "PairAnalysisHF.h"
Administrator's avatar
Administrator committed
69
#include "PairAnalysisHistos.h"
Etienne Bechtel's avatar
Etienne Bechtel committed
70
#include "PairAnalysisSignalExt.h"
71
#include "PairAnalysisStyler.h"
Administrator's avatar
Administrator committed
72
#include "PairAnalysisVarManager.h"
73
74
75

ClassImp(PairAnalysisSpectrum)

Administrator's avatar
Administrator committed
76
77
  //______________________________________________
  PairAnalysisSpectrum::PairAnalysisSpectrum()
Etienne Bechtel's avatar
Etienne Bechtel committed
78
79
  : PairAnalysisSpectrum("spectrum", "title")
{
80
81
82
83
84
85
  ///
  /// Named Constructor
  ///
}

//______________________________________________
Administrator's avatar
Administrator committed
86
87
88
89
90
PairAnalysisSpectrum::PairAnalysisSpectrum(const char* name, const char* title)
  : PairAnalysisFunction(name, title)
  , fRawInput(0)
  , fMCInput(0)
  , fMCTruth(0)
Etienne Bechtel's avatar
Etienne Bechtel committed
91
92
  , fExtractor(0)
{
93
94
95
  ///
  /// Named Constructor
  ///
Administrator's avatar
Administrator committed
96
97
98
99
  fTree = new TTree("PAPa", "PAPa-Spectrum");
  fExt  = new Extraction;
  for (Int_t i = 0; i < 100; i++)
    fInputKeys[i] = "";
100
101
102
}

//______________________________________________
Etienne Bechtel's avatar
Etienne Bechtel committed
103
104
PairAnalysisSpectrum::~PairAnalysisSpectrum()
{
105
106
107
  ///
  /// Default Destructor
  ///
Administrator's avatar
Administrator committed
108
109
  if (fResults) delete fResults;
  if (fExtractions) delete fExtractions;
110
111
112
}

//______________________________________________
Etienne Bechtel's avatar
Etienne Bechtel committed
113
114
void PairAnalysisSpectrum::AddInput(TObjArray* raw, TString identifier, TObjArray* mc, TObjArray* truth)
{
115
116
117
118
119
  ///
  /// add input array of histograms for signals extraction "raw" and
  /// efficiency calculation "mc" and unique "idetifier" string
  ///
  fRawInput.Add(raw);
Administrator's avatar
Administrator committed
120
121
122
  if (mc) fMCInput.Add(mc);
  if (truth) fMCTruth.Add(truth);
  fInputKeys[fIdx] = identifier;
123
  fIdx++;
Administrator's avatar
Administrator committed
124
}
125
126

//______________________________________________
Etienne Bechtel's avatar
Etienne Bechtel committed
127
128
void PairAnalysisSpectrum::Init()
{
129
130
131
  ///
  /// Initialize the tree
  ///
Administrator's avatar
Administrator committed
132
  fTree->Branch("Extraction", &fExt);
133
134
135
136
  //  fTree ->Branch("hSignal","TH1",&fHistSignal,32000,0);
}

//______________________________________________
Etienne Bechtel's avatar
Etienne Bechtel committed
137
138
void PairAnalysisSpectrum::Process()
{
139
140
141
142
143
  ///
  /// process the signal extraction
  ///
  Init();

Administrator's avatar
Administrator committed
144
  Int_t i = -1;
145

Administrator's avatar
Administrator committed
146
147
148
149
  TIter nextRaw(&fRawInput);  // raw content
  TObject* obj          = NULL;
  PairAnalysisHF* hf    = NULL;
  PairAnalysisHistos* h = NULL;
150

Administrator's avatar
Administrator committed
151
152
153
  TObject* objMC          = NULL;
  PairAnalysisHF* hfMC    = NULL;
  PairAnalysisHistos* hMC = NULL;
154

Administrator's avatar
Administrator committed
155
  TObject* objMCtruth = NULL;
156
157

  /// iterate over all raw input arrays
Administrator's avatar
Administrator committed
158
159
160
  while ((obj = nextRaw())) {
    i++;
    TObject::Info("Process", "Check Extraction for %s", fInputKeys[i].Data());
Etienne Bechtel's avatar
Etienne Bechtel committed
161
162
    TObject::Info("Process", "------------------------------------------------------------"
                             "-----------");
Administrator's avatar
Administrator committed
163
164
165
    TObject::Info("Process", "Input type: %s \n", obj->ClassName());

    if (fVarBinning) {
Etienne Bechtel's avatar
Etienne Bechtel committed
166
      TObject::Info("Process", "Binning provided for %s from %.3f to %.3f", fVar.Data(), fVarBinning->Min(),
Administrator's avatar
Administrator committed
167
168
169
                    fVarBinning->Max());
      //	fVarBinning->Print();
    }
170

Administrator's avatar
Administrator committed
171
172
173
174
175
176
177
    TObjArray* histArr = dynamic_cast<TObjArray*>(obj);
    if (!(histArr)) {
      if (!(h = dynamic_cast<PairAnalysisHistos*>(obj))) {
        if (!(hf = dynamic_cast<PairAnalysisHF*>(obj))) {
          TObject::Error("Process", "No input format found");
          continue;
        }
178
      }
Administrator's avatar
Administrator committed
179
    }
180

Administrator's avatar
Administrator committed
181
    // get extractor
Etienne Bechtel's avatar
Etienne Bechtel committed
182
    PairAnalysisSignalExt* sig = dynamic_cast<PairAnalysisSignalExt*>(fExtractor.At(i));
Administrator's avatar
Administrator committed
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
    if (!sig) continue;


    /// look for MC objects
    objMC                = fMCInput.At(i);
    TObjArray* histArrMC = NULL;
    if (objMC) {
      TObject::Info("Process", "Input MC type: %s \n", objMC->ClassName());
      histArrMC = dynamic_cast<TObjArray*>(objMC);
      if (!(histArrMC)) {
        if (!(hMC = dynamic_cast<PairAnalysisHistos*>(objMC))) {
          if (!(hfMC = dynamic_cast<PairAnalysisHF*>(objMC))) {
            TObject::Error("Process", "No MC input format found");
            //continue;
          }
        }
199
      }
Administrator's avatar
Administrator committed
200
    }
201

Administrator's avatar
Administrator committed
202
203
204
205
    /// look for MC truth objects
    objMCtruth            = fMCTruth.At(i);
    TObjArray* histArrMCt = NULL;
    if (objMCtruth) {
Etienne Bechtel's avatar
Etienne Bechtel committed
206
      TObject::Info("Process", "Input MC truth type: %s \n", objMCtruth->ClassName());
Administrator's avatar
Administrator committed
207
208
      histArrMCt = dynamic_cast<TObjArray*>(objMCtruth);
    }
209
210


Administrator's avatar
Administrator committed
211
212
213
214
    /// only integrated via histos
    if (!fVarBinning) {

      // get raw input histograms
Etienne Bechtel's avatar
Etienne Bechtel committed
215
216
      if (!histArr && h) histArr = h->DrawSame("pM-wghtWeight",
                                               "nomc goff");  //NOTE w/o "can" it crashes
Administrator's avatar
Administrator committed
217
218
219
      //	if(!histArr && h) histArr = h->DrawSame("pM","can nomc goff"); //NOTE w/o "can" it crashes
      //	histArr->Print();
      // get mc input histograms TODO: think about integration
Etienne Bechtel's avatar
Etienne Bechtel committed
220
      if (!histArrMC && hMC) histArrMC = hMC->DrawSame("pM", "onlymc eff goff");  //NOTE w/o "can" it crashes
Administrator's avatar
Administrator committed
221
222
223
224
225
226
227
228
229
230


      // process raw signal extraction
      sig->Process(histArr);
      if (gErrorIgnoreLevel < kWarning) sig->Print("");

      // validate signal extraction
      if (sig->GetSignal() < 0.) continue;

      // fill the tree
Etienne Bechtel's avatar
Etienne Bechtel committed
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
      fExt->setup      = fInputKeys[i];
      fExt->setupId    = i;
      fExt->poi        = sig->GetParticleOfInterest();
      fExt->var        = 1.;  //integrated
      fExt->varE       = 0.;
      fExt->s          = sig->GetSignal();
      fExt->sE         = sig->GetSignalError();
      fExt->b          = sig->GetBackground();
      fExt->bE         = sig->GetBackgroundError();
      fExt->sb         = sig->GetSB();
      fExt->sbE        = sig->GetSBError();
      fExt->sgn        = sig->GetSignificance();
      fExt->sgnE       = sig->GetSignificanceError();
      fExt->HistSignal = NULL;  //dynamic_cast<TH1F*>(sig->GetSignalHistogram());
      fExt->eff        = 1.;    //TODO: calculate this
      fExt->effE       = 0.;    //TODO: calculate this
      fExt->signal     = sig;   //NULL;
      fExt->sref       = 1.;    //TODO: calculate this
      fExt->srefE      = 0.;    //TODO: calculate this
Administrator's avatar
Administrator committed
250
251

      fTree->Fill();
Etienne Bechtel's avatar
Etienne Bechtel committed
252
253
    }
    else {
Administrator's avatar
Administrator committed
254
255
256
257

      /// get raw 2D histogram
      if (h && !h->SetCutClass(fInputKeys[i].Data())) continue;

Etienne Bechtel's avatar
Etienne Bechtel committed
258
259
      if (!histArr && h) histArr = h->DrawSame(Form("pM_%s", fVar.Data()), "nomc goff");
      TH2* histPM = (TH2*) sig->FindObject(histArr, PairAnalysis::EPairType::kSEPM);
Administrator's avatar
Administrator committed
260
261
262
263
264
265
266
267
      if (!histPM) return;

      TObjArray tmpArr;
      tmpArr.SetOwner(kFALSE);

      /// MC
      TH1* histMC = NULL;
      if (!histArrMC && hMC) {
Etienne Bechtel's avatar
Etienne Bechtel committed
268
        histArrMC = hMC->DrawSame(Form("p%s", fVar.Data()), "goff sel",
Administrator's avatar
Administrator committed
269
                                  "phi");  // TODO: add search using fPOIpdg
270
      }
Administrator's avatar
Administrator committed
271
272
273
274
275
      if (histArrMC) {
        TH1* tmpMCnom = (TH1*) histArrMC->At(0);
        if (histArrMC->GetEntriesFast() < 2) return;
        TH1* tmpMCden = (TH1*) histArrMC->At(1);
        /// rebin and calculate efficiency
Etienne Bechtel's avatar
Etienne Bechtel committed
276
277
        histMC         = tmpMCnom->Rebin(fVarBinning->GetNrows() - 1, "effMC", fVarBinning->GetMatrixArray());
        TH1* histMCden = tmpMCden->Rebin(fVarBinning->GetNrows() - 1, "effMCden", fVarBinning->GetMatrixArray());
Administrator's avatar
Administrator committed
278
279
280
281
282
283
284
285
286
287
288
289
        histMC->Divide(histMCden);
        delete histMCden;
      }
      /// debug
      if (histMC) TObject::Info("Process", "MC histogram found and rebinned");

      /// MC
      TH1* histMCtruth = NULL;
      if (histArrMCt) {
        TH1* tmpMCtrue = (TH1*) histArrMCt->At(0);
        if (!tmpMCtrue) return;
        /// rebin and calculate efficiency
Etienne Bechtel's avatar
Etienne Bechtel committed
290
        histMCtruth = tmpMCtrue->Rebin(fVarBinning->GetNrows() - 1, "sMCtrue", fVarBinning->GetMatrixArray());
Administrator's avatar
Administrator committed
291
292
      }
      /// debug
Etienne Bechtel's avatar
Etienne Bechtel committed
293
      if (histMCtruth) TObject::Info("Process", "MCtruth reference histogram found and rebinned");
Administrator's avatar
Administrator committed
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308

      // loop over all bins
      for (Int_t bin = 0; bin < fVarBinning->GetNrows() - 1; bin++) {
        tmpArr.Clear();

        // var bin limits
        Double_t xLo = fVarBinning->GetMatrixArray()[bin];
        Double_t xHi = fVarBinning->GetMatrixArray()[bin + 1] - 0.000001;
        // axis limits
        Int_t binLo = histPM->GetYaxis()->FindBin(xLo);
        Int_t binHi = histPM->GetYaxis()->FindBin(xHi);
        // found bin limits
        Double_t fndLo = histPM->GetYaxis()->GetBinLowEdge(binLo);
        Double_t fndHi = histPM->GetYaxis()->GetBinLowEdge(binHi + 1);
        //printf("binning requested, found of %s: %.3f-%.3f, %.3f-%.3f \n", fVar.Data(), xLo,xHi, fndLo,fndHi );
Etienne Bechtel's avatar
Etienne Bechtel committed
309
        TObject::Info("Process", "Bin %d: %.3f < %s < %.3f", bin, fndLo, fVar.Data(), fndHi);
Administrator's avatar
Administrator committed
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326

        // fill array for signal extraction
        for (Int_t ih = 0; ih < histArr->GetEntriesFast(); ih++) {
          TH2* hist = (TH2*) histArr->UncheckedAt(ih);
          tmpArr.Add(hist->ProjectionX(hist->GetTitle(), binLo, binHi, "e"));
        }
        //	    tmpArr.Print();

        // process
        sig->Process(&tmpArr);
        if (gErrorIgnoreLevel < kWarning) sig->Print("");

        // validate signal extraction
        if (sig->GetSignal() < 0.) continue;
        if (TMath::IsNaN(sig->GetSignalError())) continue;

        // fill the tree
Etienne Bechtel's avatar
Etienne Bechtel committed
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
        fExt->setup      = fInputKeys[i];
        fExt->setupId    = i;
        fExt->poi        = sig->GetParticleOfInterest();
        fExt->var        = (fndHi - fndLo) / 2 + fndLo;  // center of the found bin
        fExt->varE       = (fndHi - fndLo) / 2;          // found bin width
        fExt->s          = sig->GetSignal();
        fExt->sE         = sig->GetSignalError();
        fExt->b          = sig->GetBackground();
        fExt->bE         = sig->GetBackgroundError();
        fExt->sb         = sig->GetSB();
        fExt->sbE        = sig->GetSBError();
        fExt->sgn        = sig->GetSignificance();
        fExt->sgnE       = sig->GetSignificanceError();
        fExt->HistSignal = NULL;  //dynamic_cast<TH1F*>(sig->GetSignalHistogram());
        fExt->eff        = (histMC ? histMC->GetBinContent(bin + 1) : 1.);
        fExt->effE       = (histMC ? histMC->GetBinError(bin + 1) : 0.);
        fExt->signal     = sig;
        fExt->sref       = (histMCtruth ? histMCtruth->GetBinContent(bin + 1) : 0.);
        fExt->srefE      = (histMCtruth ? histMCtruth->GetBinError(bin + 1) : 0.);
Administrator's avatar
Administrator committed
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367

        fTree->Fill();

        // set variable
        // hf->AddCutVariable( (EValueTypes)
        // 		     PairAnalysisVarManager::GetValueType(fVar.Data()),
        // 		     fVarBinning->GetMatrixArray()[bin],
        // 		     fVarBinning->GetMatrixArray()[bin+1]
        // 		     );

        // get histogram array
        //...

      }  //end binning

      /// clean up
      if (histMC) delete histMC;
      if (histMCtruth) delete histMCtruth;

    }  // end 2D

  }  // end raw content
368
369
370
371
372

  //  fTree->Print();
}

//______________________________________________
Etienne Bechtel's avatar
Etienne Bechtel committed
373
374
void PairAnalysisSpectrum::DrawSpectrum(const char* varexp, const char* selection, Option_t* option)
{
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
  //
  // TTree draw alias
  //
  /// additional plotting options:
  ///
  /// "syst":       draw and calculate systematic uncertainties according to fSystMthd (central point=mean value)
  /// "samepad":    draws spectrum into current pad
  ///
  /// "logx,y,z":   the axis are plotted in log-scale (labels are added automatically according to the range)
  /// "leg(f)":     a ("filled") legend will be created with caption=className ,
  ///               can be modified by PairAnalysisHistos::SetName("mycaption"),
  ///               change of legend position: see PairAnalysisStyler::SetLegendAlign
  ///
  /// "print":      prints the data points to stdout
  ///

  TString optString(option);
  optString.ToLower();
Etienne Bechtel's avatar
Etienne Bechtel committed
393
  printf("Plot spectrum: '%s' \t selection: '%s' \t options: '%s' \n", varexp, selection, optString.Data());
Administrator's avatar
Administrator committed
394
395
396
397
398
399
400
401
402
403
  Bool_t optLegFull = optString.Contains("legf");
  optString.ReplaceAll("legf", "");
  Bool_t optLeg = optString.Contains("leg");
  optString.ReplaceAll("leg", "");
  Bool_t optSyst = optString.Contains("syst");
  optString.ReplaceAll("syst", "");
  Bool_t optPrint = optString.Contains("print");
  optString.ReplaceAll("print", "");
  Bool_t optSamePad = optString.Contains("samepad");
  optString.ReplaceAll("samepad", "");
404
405

  /// counter
Administrator's avatar
Administrator committed
406
  Long64_t n = 1;
407
408
409
410
411
412
413
414
415

  /// load style
  PairAnalysisStyler::LoadStyle();

  /// canvas key
  TString ckey(varexp);
  Int_t ndim = ckey.CountChar(':') + 1;

  /// for variable/formula you look at, e.g. signal/eff
Administrator's avatar
Administrator committed
416
  TObjArray* carr = ckey.Tokenize(":");
417
418
419
420
  carr->SetOwner();
  //  delete carr;

  /// first variable e.g. signal in signal/eff
Administrator's avatar
Administrator committed
421
  TObjArray* oarr = ckey.Tokenize("+-*/:");
422
  oarr->SetOwner();
Administrator's avatar
Administrator committed
423
  TString fkey = ((TObjString*) oarr->At(0))->GetString();
424
425
426
  delete oarr;

  /// canvas
Administrator's avatar
Administrator committed
427
428
  ckey.ReplaceAll("/", "#");  /// canvas name does not allow '/'
  TCanvas* c = NULL;
Etienne Bechtel's avatar
Etienne Bechtel committed
429
  if (optSamePad) c = gPad->GetCanvas();
Administrator's avatar
Administrator committed
430
431
432
433
434
  else
    c = (TCanvas*) gROOT->FindObject(ckey.Data());
  if (!c) {
    TObject::Info("Draw", "create new canvas: '%s'", ckey.Data());
    c = new TCanvas(ckey.Data(), ckey.Data());
435
436
437
438
  }
  c->cd();

  /// count number of drawn objects in pad
Administrator's avatar
Administrator committed
439
440
441
  TObject* obj;
  Int_t nobj  = 0;
  TList* prim = gPad->GetListOfPrimitives();
442
  //////    if(prim->GetSize()>1) prim->RemoveLast(); // remove redraw axis histogram
Administrator's avatar
Administrator committed
443
444
445
  for (Int_t io = 0; io < prim->GetSize(); io++) {
    obj = prim->At(io);
    if (obj->InheritsFrom(TGraph::Class()) && obj != prim->At(io + 1)) nobj++;
446
447
  }

Administrator's avatar
Administrator committed
448
449
450
451
452
453
  TLegend* leg = 0;
  if ((optLeg && !nobj)) {
    //  if ( (optLeg && optTask && !nobj) || (optLeg && !optTask && !optDet) ) {
    leg = new TLegend(0. + gPad->GetLeftMargin() + gStyle->GetTickLength("Y"),
                      0. + gPad->GetBottomMargin() + gStyle->GetTickLength("X"),
                      1. - gPad->GetRightMargin() + gStyle->GetTickLength("Y"),
Etienne Bechtel's avatar
Etienne Bechtel committed
454
455
456
                      1. - gPad->GetTopMargin() + gStyle->GetTickLength("X"), GetName(), "nbNDC");
  }
  else if (optLeg && nobj) {
Administrator's avatar
Administrator committed
457
    leg = (TLegend*) prim->FindObject("TPave");
458
459
  }

Administrator's avatar
Administrator committed
460
  Info("Draw", "Basics: nobj: %d \t leg: %p", nobj, leg);
461
462

  // logaritmic style
Administrator's avatar
Administrator committed
463
464
465
466
467
468
  if (optString.Contains("logx")) gPad->SetLogx();
  if (optString.Contains("logy")) gPad->SetLogy();
  if (optString.Contains("logz")) gPad->SetLogz();
  optString.ReplaceAll("logx", "");
  optString.ReplaceAll("logy", "");
  optString.ReplaceAll("logz", "");
469

Etienne Bechtel's avatar
Etienne Bechtel committed
470
  if (ndim < 3 && !ckey.Contains("signal->")) {  // build own histogram with labels
471
472
473
474
475

    //    printf("try to get errors for %d dimensional input \n",ndim);
    TString varkey(varexp);

    // get event list
Administrator's avatar
Administrator committed
476
477
478
479
480
481
482
483
    n = fTree->Draw(">>elist", selection);
    if (!n) {
      delete carr;
      return;
    }
    TEventList* elist = (TEventList*) gDirectory->Get("elist");
    if (elist) {
      elist->SetReapplyCut(kTRUE);  // important!
484
485
486
487
      //	elist->Print("all");
      fTree->SetEventList(elist);
    }

Administrator's avatar
Administrator committed
488
    TH1D* hist = 0x0;
489
    // set up a proper histogram with labels
Administrator's avatar
Administrator committed
490
    if (varkey.Contains("setup")) {
491

Administrator's avatar
Administrator committed
492
493
      hist = new TH1D();  // activate buffer
      hist->SetNameTitle(varexp, selection);
494
495

      //read strings for all entries
Administrator's avatar
Administrator committed
496
497
498
499
      for (Long64_t i = 0; i < n; i++) {
        fTree->GetEntry(i);
        hist->Fill((fExt->setup).Data(), 1.);
        //	  printf(" read from tree: %s \n",(fExt->setup).Data());
500
501
502
      }

      hist->Draw("AXIS");
Administrator's avatar
Administrator committed
503
      hist->GetXaxis()->SetRange(1, n);
504
505
506
    }

    // get errors from tree
Administrator's avatar
Administrator committed
507
508
    TString errkey = ((TObjString*) carr->At(0))->GetString();
    errkey.ReplaceAll(fkey, fkey + "E");  /// first variables error
509
510
511
512
    //    if(ndim>1) varkey.ReplaceAll(fkey+":",fkey+":"+fkey+"E:");
    //    else       varkey.ReplaceAll(fkey,fkey+":"+fkey+"E");
    // if(ndim>1) varkey.ReplaceAll(fkey+":",fkey+":"+fkey+"E:");
    // else
Administrator's avatar
Administrator committed
513
    if (!fkey.Contains("E")) {
514
      varkey.Append(errkey.Prepend(":"));
Etienne Bechtel's avatar
Etienne Bechtel committed
515
      Info("Draw", " Appended '%s' by '%s' for error caluclation", varkey.Data(), errkey.Data());
516
517
518
519
520
521
522
    }

    /// TTree:Draw command with graphics off to get errors
    //printf("execute collect/draw command for %s \n",varkey.Data());
    fTree->Draw(varkey, selection, "goff");

    /// for one dimensional histograms
Administrator's avatar
Administrator committed
523
524
525
    Double_t* xval = new Double_t[fTree->GetSelectedRows()];
    for (Int_t ix = 0; ix < fTree->GetSelectedRows(); ix++)
      xval[ix] = 1.;
526
527

    /// setup the final tgraph WITH error bars
Administrator's avatar
Administrator committed
528
    TGraphErrors* gr = 0x0;
Etienne Bechtel's avatar
Etienne Bechtel committed
529
    if (ndim > 1) gr = new TGraphErrors(fTree->GetSelectedRows(), fTree->GetV2(), fTree->GetV1(), 0, fTree->GetV3());
Administrator's avatar
Administrator committed
530
    else
Etienne Bechtel's avatar
Etienne Bechtel committed
531
      gr = new TGraphErrors(fTree->GetSelectedRows(), fTree->GetV1(), xval, fTree->GetV2(), 0);
Administrator's avatar
Administrator committed
532
533
534
535
536
537
    delete[] xval;

    if (!gr) {
      delete carr;
      return;
    }
538
539

    /// style and legend entry
Administrator's avatar
Administrator committed
540
    PairAnalysisStyler::Style(gr, nobj);
541
    TString sel(selection);
Administrator's avatar
Administrator committed
542
543
544
545
546
547
548
549
    sel.ReplaceAll("setup.Contains", "");
    sel.ReplaceAll("(\"", "*");
    sel.ReplaceAll("\")", "*");
    sel.ReplaceAll("setup==", "");
    sel.ReplaceAll("\"", "");
    gr->SetName(Form("%s", sel.Data()));
    if (sel.Contains("setupId==") && sel.Length() < 11) {
      sel.ReplaceAll("setupId==", "");
550
      Int_t iId = sel.Atoi();
Administrator's avatar
Administrator committed
551
      gr->SetName(Form("%s", fInputKeys[iId].Data()));
552
    }
Administrator's avatar
Administrator committed
553
554
    if (optSyst) gr->SetName(GetTitle());
    if (fkey.Contains("ref")) gr->SetName("MC-truth");
555
556
557

    // sort x-values
    gr->Sort();
Etienne Bechtel's avatar
Etienne Bechtel committed
558
559
560
561
    TGraphErrors* grE = NULL;              // statistical
    TGraphErrors* grS = NULL;              // systematic
    TGraphErrors* grC = NULL;              // stat + syst
    if (optSyst && fVarBinning) {          //TODO: implement systematic calculation w/o binning
Administrator's avatar
Administrator committed
562
563
564
565
566
567
      grE           = new TGraphErrors();  // statistical graph
      grS           = new TGraphErrors();  // systematics graph
      grC           = new TGraphErrors();  // systematics graph
      Double_t* gx  = gr->GetX();
      Double_t* gy  = gr->GetY();
      Double_t* gye = gr->GetEY();
568
569

      // loop over all variable bins
Administrator's avatar
Administrator committed
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
      Int_t first = 0;  //TMath::BinarySearch(gr->GetN(),gx,xLo); //first bin
      Int_t ibin  = 0;
      while (first < gr->GetN()) {
        //      for(ibin=0; ibin<fVarBinning->GetNrows()-1; ibin++) {

        Int_t nsys    = 0;   // counter #systematics
        Double_t ysys = 0.;  // y-vlaue of systematic = mean value
        Double_t esys = 0.;  // uncertainty depends on method
        // calculate mean
        Double_t xvar = gx[first];
        for (Int_t i = first; i < gr->GetN(); i++) {
          //	  printf("gx[i]:%f xvar:%f xLo:%f\n",gx[i],xvar,xLo);
          //	  if(TMath::Abs(gx[i]-xvar)>1.e-8) break;
          if ((gx[i] - xvar) > 1.e-8) break;
          // printf("graph entry %d, found index first: %d with value %.3f \t y: %.3f+-%.3f \n",
          //        i,first,xvar,gy[i],gye[i]);
          nsys++;
          ysys += gy[i];
        }
        ysys /= (nsys ? nsys : 1);  // protect for zero division
        //	printf("bin %d, found index first: %d, %.3f \t y:%3.f \t nsys:%d\n",ibin,first,xvar,ysys,nsys);

        // y-syst. uncertainty
        //	printf("syst %f , <y> %f\n",esys,ysys);
        for (Int_t i = 0; i < nsys; i++) {
          Int_t j = first + i;
          //	  if(gx[j]!=xLo) break; // check should not be needed
          Double_t uce = 0.;
          switch (fSystMthd) {
599
            case ESystMethod::kBarlow:
Administrator's avatar
Administrator committed
600
              // I.  calc uncorr. stat. error from sub/superset w.r.t. first measurement
Etienne Bechtel's avatar
Etienne Bechtel committed
601
              uce = TMath::Sqrt(TMath::Abs(gye[j] * gye[j] - gye[first] * gye[first]));
Administrator's avatar
Administrator committed
602
603
604
605
              // II. calc max. deviation w.r.t. mean y-value incl. 1sigma* 0.9 of I.
              // NOTE: 0.9 can be change to a max value of 1->1sigma, 0.9 is more consevative
              esys = TMath::Max(esys, TMath::Abs(ysys - gy[j]) - 0.9 * uce);
              break;
Etienne Bechtel's avatar
Etienne Bechtel committed
606
            case ESystMethod::kSystMax: esys = TMath::Max(esys, TMath::Abs(ysys - gy[j])); break;
607
            case ESystMethod::kSystRMS: esys += gy[j] * gy[j]; break;
Administrator's avatar
Administrator committed
608
609
610
611
612
613
          }
          //	  printf("bin error %f \t  syst %f  from abs %f \n",gye[j],esys, TMath::Abs( gy[j] ));
        }

        // normalisation
        switch (fSystMthd) {
614
615
          case ESystMethod::kBarlow: /* nothing to be done */ break;
          case ESystMethod::kSystMax: /* nothing to be done */ break;
Etienne Bechtel's avatar
Etienne Bechtel committed
616
          case ESystMethod::kSystRMS: esys = TMath::Sqrt(TMath::Abs(esys / (nsys ? nsys : 1) - ysys * ysys)); break;
Administrator's avatar
Administrator committed
617
618
619
620
621
622
        }

        // fill statistical and systematic graph values and errors
        grE->SetPoint(ibin, xvar, ysys);            // mean
        grE->SetPointError(ibin, 0.0, gye[first]);  // stat.uncert. of first set

Etienne Bechtel's avatar
Etienne Bechtel committed
623
        Double_t boxW = (fVarBinning->Max() - fVarBinning->Min()) / (fVarBinning->GetNrows() - 1);
Administrator's avatar
Administrator committed
624
625
626
627
628
        grS->SetPoint(ibin, xvar, ysys);              // mean
        grS->SetPointError(ibin, boxW * 0.35, esys);  // systematic value

        // calculate err = sqrt(stat.**2 + syst**2)
        grC->SetPoint(ibin, xvar, ysys);  // mean
Etienne Bechtel's avatar
Etienne Bechtel committed
629
        grC->SetPointError(ibin, boxW * 0.35, TMath::Sqrt(esys * esys + gye[first] * gye[first]));
Administrator's avatar
Administrator committed
630
631
632
633
634

        // increase index counter
        first += nsys;
        ibin++;
      }  //next bin
635
      //      grS->Print();
Etienne Bechtel's avatar
Etienne Bechtel committed
636
637
    }
    else {
638
639
640
      grC = new TGraphErrors(*gr);
    }

Administrator's avatar
Administrator committed
641
    if (optPrint) grC->Print();
642

Administrator's avatar
Administrator committed
643
    Info("Draw", " Draw object with options: '%s'", optString.Data());
Etienne Bechtel's avatar
Etienne Bechtel committed
644
    if (!PairAnalysisStyler::GetFirstHistogram()) gr->Draw((optString + "A").Data());
645
    else {
Administrator's avatar
Administrator committed
646
      gr->Draw((optString + "same").Data());
647
648

      // set axis maximum
Administrator's avatar
Administrator committed
649
650
651
      Double_t* valE = grC->GetEY();
      Double_t* val  = grC->GetY();
      Int_t npnts    = grC->GetN();
652
      Int_t idx[1000];
Administrator's avatar
Administrator committed
653
      TMath::Sort(npnts, val, idx, kTRUE);  // kFALSE=increasing numbers
654

Etienne Bechtel's avatar
Etienne Bechtel committed
655
      Double_t errmin = (TMath::IsNaN(valE[idx[npnts - 1]]) ? 0. : valE[idx[npnts - 1]]);
Administrator's avatar
Administrator committed
656
      Double_t min    = (val[idx[npnts - 1]] - errmin) * 0.9;
657
      Double_t tmpmin = PairAnalysisStyler::GetFirstHistogram()->GetMinimum();
Etienne Bechtel's avatar
Etienne Bechtel committed
658
      PairAnalysisStyler::GetFirstHistogram()->SetMinimum((tmpmin < min ? tmpmin : min));
Administrator's avatar
Administrator committed
659
660
      Double_t errmax = (TMath::IsNaN(valE[idx[0]]) ? 0. : valE[idx[0]]);
      Double_t max    = (val[idx[0]] + errmax) * 1.1;
661
      Double_t tmpmax = PairAnalysisStyler::GetFirstHistogram()->GetMaximum();
Etienne Bechtel's avatar
Etienne Bechtel committed
662
      PairAnalysisStyler::GetFirstHistogram()->SetMaximum((tmpmax > max ? tmpmax : max));
663
664
665
    }

    // draw systemtaic graph ontop
Administrator's avatar
Administrator committed
666
667
668
669
    if (grS) {
      PairAnalysisStyler::Style(grE, nobj);
      grE->Draw((optString + "A").Data());
      PairAnalysisStyler::Style(grS, nobj);
670
671
672
673
674
675
      grS->SetFillColor(grS->GetLineColor());
      grS->SetFillStyle(kFEmpty);
      grS->Draw("2same");
    }

    // legend
Administrator's avatar
Administrator committed
676
677
678
679
680
681
    TString legOpt = optString + "L";
    legOpt.ReplaceAll("hist", "");
    legOpt.ReplaceAll("scat", "");
    if (legOpt.Contains("col")) legOpt = "";
    legOpt.ReplaceAll("z", "");
    legOpt.ReplaceAll("e", "");
682
683

    TString legkey = gr->GetName();
Administrator's avatar
Administrator committed
684
    if (leg) leg->AddEntry(gr, gr->GetName(), legOpt.Data());
685
686

    fSignal = grC;
Administrator's avatar
Administrator committed
687
    PairAnalysisStyler::Style(fSignal, nobj);
Etienne Bechtel's avatar
Etienne Bechtel committed
688
689
  }
  else {
690
691
    // execute tree draw command
    fTree->Draw(varexp, selection, optString.Data());
Administrator's avatar
Administrator committed
692
    fprintf(stderr, "use plain TTree::Draw command \n");
693
694
695
696
697
698
    return;
  }


  /// modify axis and titles
  //  printf("modify axis titles \n");
Administrator's avatar
Administrator committed
699
  if (!optSamePad) {
700
701
    UInt_t varx = PairAnalysisVarManager::GetValueType(fVar.Data());
    TString var(varexp);
Administrator's avatar
Administrator committed
702
    TObjArray* arr = var.Tokenize(":");
703
    arr->SetOwner();
Administrator's avatar
Administrator committed
704
705
706
    TString xt = "";
    TString yt = "Entries";
    xt         = ((TObjString*) arr->At(0))->GetString();
Etienne Bechtel's avatar
Etienne Bechtel committed
707
    if (xt.EqualTo("sb")) xt = PairAnalysisSignalExt::GetValueName(3);
Administrator's avatar
Administrator committed
708
709
710
711
712
713
714
    else if (xt.EqualTo("s"))
      xt = PairAnalysisSignalExt::GetValueName(0);
    else if (xt.EqualTo("b"))
      xt = PairAnalysisSignalExt::GetValueName(1);
    else if (xt.EqualTo("sgn"))
      xt = PairAnalysisSignalExt::GetValueName(2);
    else if (xt.EqualTo("var"))
Etienne Bechtel's avatar
Etienne Bechtel committed
715
      xt = Form("%s %s", PairAnalysisVarManager::GetValueLabel(varx), PairAnalysisVarManager::GetValueUnit(varx));
Administrator's avatar
Administrator committed
716
    else if (xt.Contains("var")) {
Etienne Bechtel's avatar
Etienne Bechtel committed
717
      xt.ReplaceAll("varE", Form("#Delta%s", PairAnalysisVarManager::GetValueLabel(varx)));
718
719
720
      xt.ReplaceAll("var", PairAnalysisVarManager::GetValueLabel(varx));
    }

Administrator's avatar
Administrator committed
721
    if (arr->GetEntriesFast() < 2) {
722
723
      PairAnalysisStyler::GetFirstHistogram()->SetXTitle(xt.Data());
      PairAnalysisStyler::GetFirstHistogram()->SetYTitle(yt.Data());
Etienne Bechtel's avatar
Etienne Bechtel committed
724
725
    }
    else {
726
      PairAnalysisStyler::GetFirstHistogram()->SetYTitle(xt.Data());
Administrator's avatar
Administrator committed
727
      xt = ((TObjString*) arr->At(1))->GetString();
Etienne Bechtel's avatar
Etienne Bechtel committed
728
      if (xt.EqualTo("sb")) xt = PairAnalysisSignalExt::GetValueName(3);
Administrator's avatar
Administrator committed
729
730
731
732
733
734
735
      else if (xt.EqualTo("s"))
        xt = PairAnalysisSignalExt::GetValueName(0);
      else if (xt.EqualTo("b"))
        xt = PairAnalysisSignalExt::GetValueName(1);
      else if (xt.EqualTo("sgn"))
        xt = PairAnalysisSignalExt::GetValueName(2);
      else if (xt.EqualTo("var"))
Etienne Bechtel's avatar
Etienne Bechtel committed
736
        xt = Form("%s %s", PairAnalysisVarManager::GetValueLabel(varx), PairAnalysisVarManager::GetValueUnit(varx));
Administrator's avatar
Administrator committed
737
      else if (xt.Contains("var")) {
Etienne Bechtel's avatar
Etienne Bechtel committed
738
        xt.ReplaceAll("varE", Form("#Delta%s", PairAnalysisVarManager::GetValueLabel(varx)));
Administrator's avatar
Administrator committed
739
        xt.ReplaceAll("var", PairAnalysisVarManager::GetValueLabel(varx));
740
741
742
743
744
      }
      PairAnalysisStyler::GetFirstHistogram()->SetXTitle(xt.Data());
    }

    /// delete array
Administrator's avatar
Administrator committed
745
    if (arr) delete arr;
746
747
748
  }

  /// clean up
Administrator's avatar
Administrator committed
749
  if (carr) delete carr;
750
751

  /// set ndivisions
Etienne Bechtel's avatar
Etienne Bechtel committed
752
  if (fVarBinning) PairAnalysisStyler::GetFirstHistogram()->SetAxisRange(fVarBinning->Min(), fVarBinning->Max(), "X");
753
754
755
756
757
  //    PairAnalysisStyler::GetFirstHistogram()->GetXaxis()->SetNdivisions(, 0, 0, kFALSE);


  // legend
  if (leg) {
Etienne Bechtel's avatar
Etienne Bechtel committed
758
759
    PairAnalysisStyler::SetLegendAttributes(leg, optLegFull);  // coordinates, margins, fillstyle, fontsize
    if (!nobj) leg->Draw();                                    // only draw the legend once
760
761
762
763
764
765
766
767
    ///    gPad->GetCanvas()->Update();
  }

  /// release eventlist
  fTree->SetEventList(0);
}

//______________________________________________
Etienne Bechtel's avatar
Etienne Bechtel committed
768
769
void PairAnalysisSpectrum::Write()
{
770
771
772
  ///
  /// write to output file
  ///
773
774
775
776
777

  /// Save old global file and folder pointer to avoid messing with FairRoot
  TFile* oldFile     = gFile;
  TDirectory* oldDir = gDirectory;

Administrator's avatar
Administrator committed
778
  TFile* fout = new TFile("test.root", "RECREATE");
779
780
781
782
  fout->cd();
  //  fTree->Print();
  fTree->Write();
  fout->Close();
783
784
785
786

  /// Restore old global file and folder pointer to avoid messing with FairRoot
  gFile      = oldFile;
  gDirectory = oldDir;
787
788
789
}

//______________________________________________
Etienne Bechtel's avatar
Etienne Bechtel committed
790
791
void PairAnalysisSpectrum::Fit(TString drawoption)
{
792
793
794
795
796
797
798
799
800
  ///
  /// Fit the spectrum according to the selected method
  ///
  /// additional plotting options:
  ///
  /// "leg":     add fit to legend if it exists
  ///

  drawoption.ToLower();
Administrator's avatar
Administrator committed
801
802
  Bool_t optLeg = drawoption.Contains("leg");
  drawoption.ReplaceAll("leg", "");
803
804
805

  //  fSignal->Print();

Administrator's avatar
Administrator committed
806
807
  Info("Fit", "Spectrum fit method: %s", fFuncSigBack->GetName());
  Int_t fitResult = fSignal->Fit(fFuncSigBack, (fFitOpt + "EX0").Data());
808
809

  // warning in case of fit issues
Administrator's avatar
Administrator committed
810
811
812
813
  if (fitResult != 0) {
    Error("Fit", "fit has error/issue (%d)", fitResult);
    return;
  }
814

Etienne Bechtel's avatar
Etienne Bechtel committed
815
  PairAnalysisStyler::Style(fFuncSigBack, static_cast<Int_t>(PairAnalysisStyler::Eidx::kFit));
Administrator's avatar
Administrator committed
816
  fFuncSigBack->SetLineColor(fSignal->GetLineColor());
817
818
819
  //  fFuncSigBack->SetLineStyle(kDashed);

  //  PairAnalysisStyler::Style(fit, PairAnalysisStyler::kFit);
Administrator's avatar
Administrator committed
820
  TF1* fit = fFuncSigBack->DrawCopy((drawoption + "same").Data());
821
822

  /// store chi2/ndf of the fit
Administrator's avatar
Administrator committed
823
824
  fDof = fFuncSigBack->GetNDF();
  if (fDof) fChi2Dof = fFuncSigBack->GetChisquare() / fFuncSigBack->GetNDF();
825
826

  /// add fit to legend
Administrator's avatar
Administrator committed
827
828
829
  if (optLeg) {
    TList* prim    = gPad->GetListOfPrimitives();
    TLegend* leg   = (TLegend*) prim->FindObject("TPave");
830
    TString legkey = fFuncSigBack->GetName();
Administrator's avatar
Administrator committed
831
832
833
834
835
836
    /// recalc legend coordinates, margins
    if (leg) {
      leg->AddEntry(fit, fFuncSigBack->GetName(), drawoption.Data());
      PairAnalysisStyler::SetLegendAttributes(leg);
      ///leg->Draw(); // was w/o !nobj
    }
837
838
  }
}