TPC & TRD calibration viewer GUI, base classes (Ionut)
[u/mrichter/AliRoot.git] / STEER / AliBaseCalibViewer.cxx
1 ///////////////////////////////////////////////////////////////////////////////
2 //                                                                           //
3 //  Base class for the AliTPCCalibViewer and AliTRDCalibViewer               //
4 //  used for the calibration monitor                                         //
5 //                                                                           //
6 //  Authors:     Marian Ivanov (Marian.Ivanov@cern.ch)                       //
7 //               Jens Wiechula (Jens.Wiechula@cern.ch)                       //
8 //               Ionut Arsene  (iarsene@cern.ch)                             //
9 //                                                                           //
10 ///////////////////////////////////////////////////////////////////////////////
11
12
13 #include <iostream>
14 #include <fstream>
15 #include <TString.h>
16 #include <TRandom.h>
17 #include <TLegend.h>
18 #include <TLine.h>
19 //#include <TCanvas.h>
20 #include <TROOT.h>
21 #include <TStyle.h>
22 #include <TH1.h> 
23 #include <TH1F.h>
24 #include <TMath.h>
25 #include <THashTable.h>
26 #include <TObjString.h>
27 #include <TLinearFitter.h>
28 #include <TTreeStream.h>
29 #include <TFile.h>
30 #include <TKey.h>
31 #include <TGraph.h>
32 #include <TDirectory.h>
33 #include <TFriendElement.h>
34
35 #include "AliBaseCalibViewer.h"
36
37 ClassImp(AliBaseCalibViewer)
38
39 AliBaseCalibViewer::AliBaseCalibViewer()
40                   :TObject(),
41                    fTree(0),
42                    fFile(0),
43                    fListOfObjectsToBeDeleted(0),
44                    fTreeMustBeDeleted(0), 
45                    fAbbreviation(0), 
46                    fAppendString(0)
47 {
48   //
49   // Default constructor
50   //
51 }
52
53 //_____________________________________________________________________________
54 AliBaseCalibViewer::AliBaseCalibViewer(const AliBaseCalibViewer &c)
55                   :TObject(c),
56                    fTree(0),
57                    fFile(0),
58                    fListOfObjectsToBeDeleted(0),
59                    fTreeMustBeDeleted(0),
60                    fAbbreviation(0), 
61                    fAppendString(0)
62 {
63   //
64   // dummy AliBaseCalibViewer copy constructor
65   // not yet working!!!
66   //
67   fTree = c.fTree;
68   fTreeMustBeDeleted = c.fTreeMustBeDeleted;
69   fListOfObjectsToBeDeleted = c.fListOfObjectsToBeDeleted;
70   fAbbreviation = c.fAbbreviation;
71   fAppendString = c.fAppendString;
72 }
73
74 //_____________________________________________________________________________
75 AliBaseCalibViewer::AliBaseCalibViewer(TTree* tree)
76                   :TObject(),
77                    fTree(0),
78                    fFile(0),
79                    fListOfObjectsToBeDeleted(0),
80                    fTreeMustBeDeleted(0),
81                    fAbbreviation(0), 
82                    fAppendString(0)
83 {
84   //
85   // Constructor that initializes the calibration viewer
86   //
87   fTree = tree;
88   fTreeMustBeDeleted = kFALSE;
89   fListOfObjectsToBeDeleted = new TObjArray();
90   fAbbreviation = "~";
91   fAppendString = ".fElements";
92 }
93
94 //_____________________________________________________________________________
95 AliBaseCalibViewer::AliBaseCalibViewer(const Char_t* fileName, const Char_t* treeName)
96                   :TObject(),
97                    fTree(0),
98                    fFile(0),
99                    fListOfObjectsToBeDeleted(0),
100                    fTreeMustBeDeleted(0),
101                    fAbbreviation(0), 
102                    fAppendString(0)
103                    
104 {
105    //
106    // Constructor to initialize the calibration viewer
107    // the file 'fileName' contains the tree 'treeName'
108    //
109    fFile = new TFile(fileName, "read");
110    fTree = (TTree*) fFile->Get(treeName);
111    fTreeMustBeDeleted = kTRUE;
112    fListOfObjectsToBeDeleted = new TObjArray();
113    fAbbreviation = "~";
114    fAppendString = ".fElements";
115 }
116
117 //____________________________________________________________________________
118 AliBaseCalibViewer & AliBaseCalibViewer::operator =(const AliBaseCalibViewer & param)
119 {
120    //
121    // assignment operator - dummy
122    // not yet working!!!
123    //
124    fTree = param.fTree;
125    fTreeMustBeDeleted = param.fTreeMustBeDeleted;
126    fListOfObjectsToBeDeleted = param.fListOfObjectsToBeDeleted;
127    fAbbreviation = param.fAbbreviation;
128    fAppendString = param.fAppendString;
129    return (*this);
130 }
131
132 //_____________________________________________________________________________
133 AliBaseCalibViewer::~AliBaseCalibViewer()
134 {
135    //
136    // AliBaseCalibViewer destructor
137    // all objects will be deleted, the file will be closed, the pictures will disappear
138    //
139    if (fTree && fTreeMustBeDeleted) {
140       fTree->SetCacheSize(0);
141       fTree->Delete();
142    }
143    if (fFile) {
144       fFile->Close();
145       fFile = 0;
146    }
147
148    for (Int_t i = fListOfObjectsToBeDeleted->GetEntriesFast()-1; i >= 0; i--) {
149       delete fListOfObjectsToBeDeleted->At(i);
150    }
151    delete fListOfObjectsToBeDeleted;
152 }
153
154 //_____________________________________________________________________________
155 void AliBaseCalibViewer::Delete(Option_t* option) {
156    //
157    // Should be called from AliBaseCalibViewerGUI class only.
158    // If you use Delete() do not call the destructor.
159    // All objects (except those contained in fListOfObjectsToBeDeleted) will be deleted, the file will be closed.
160    //
161    
162    option = option;  // to avoid warnings on compiling   
163    if (fTree && fTreeMustBeDeleted) {
164       fTree->SetCacheSize(0);
165       fTree->Delete();
166    }
167    if (fFile)
168       delete fFile;
169    delete fListOfObjectsToBeDeleted;
170 }
171
172 //_____________________________________________________________________________
173 void AliBaseCalibViewer::FormatHistoLabels(TH1 *histo) const {
174    // 
175    // formats title and axis labels of histo 
176    // removes '.fElements'
177    // 
178    if (!histo) return;
179    TString replaceString(fAppendString.Data());
180    TString *str = new TString(histo->GetTitle());
181    str->ReplaceAll(replaceString, "");
182    histo->SetTitle(str->Data());
183    delete str;
184    if (histo->GetXaxis()) {
185       str = new TString(histo->GetXaxis()->GetTitle());
186       str->ReplaceAll(replaceString, "");
187       histo->GetXaxis()->SetTitle(str->Data());
188       delete str;
189    }
190    if (histo->GetYaxis()) {
191       str = new TString(histo->GetYaxis()->GetTitle());
192       str->ReplaceAll(replaceString, "");
193       histo->GetYaxis()->SetTitle(str->Data());
194       delete str;
195    }
196    if (histo->GetZaxis()) {
197       str = new TString(histo->GetZaxis()->GetTitle());
198       str->ReplaceAll(replaceString, "");
199       histo->GetZaxis()->SetTitle(str->Data());
200       delete str;
201    }
202 }
203
204 //_____________________________________________________________________________
205 TFriendElement* AliBaseCalibViewer::AddReferenceTree(const Char_t* filename, const Char_t* treename, const Char_t* refname){
206   //
207   // add a reference tree to the current tree
208   // by default the treename is 'tree' and the reference treename is 'R'
209   //
210    TFile *file = new TFile(filename);
211    fListOfObjectsToBeDeleted->Add(file);
212    TTree * tree = (TTree*)file->Get(treename);
213    return AddFriend(tree, refname);
214 }
215
216 //_____________________________________________________________________________
217 TString* AliBaseCalibViewer::Fit(const Char_t* drawCommand, const Char_t* formula, const Char_t* cuts, 
218                                  Double_t & chi2, TVectorD &fitParam, TMatrixD &covMatrix){
219    //
220    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
221    // returns chi2, fitParam and covMatrix
222    // returns TString with fitted formula
223    //
224    
225    TString formulaStr(formula); 
226    TString drawStr(drawCommand);
227    TString cutStr(cuts);
228    
229    // abbreviations:
230    drawStr.ReplaceAll(fAbbreviation, fAppendString);
231    cutStr.ReplaceAll(fAbbreviation, fAppendString);
232    formulaStr.ReplaceAll(fAbbreviation, fAppendString);
233    
234    formulaStr.ReplaceAll("++", fAbbreviation);
235    TObjArray* formulaTokens = formulaStr.Tokenize(fAbbreviation.Data()); 
236    Int_t dim = formulaTokens->GetEntriesFast();
237    
238    fitParam.ResizeTo(dim);
239    covMatrix.ResizeTo(dim,dim);
240    
241    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
242    fitter->StoreData(kTRUE);   
243    fitter->ClearPoints();
244    
245    Int_t entries = Draw(drawStr.Data(), cutStr.Data(), "goff");
246    if (entries == -1) return new TString("An ERROR has occured during fitting!");
247    Double_t **values = new Double_t*[dim+1] ; 
248    
249    for (Int_t i = 0; i < dim + 1; i++){
250       Int_t centries = 0;
251       if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff");
252       else  centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff");
253       
254       if (entries != centries) return new TString("An ERROR has occured during fitting!");
255       values[i] = new Double_t[entries];
256       memcpy(values[i],  fTree->GetV1(), entries*sizeof(Double_t)); 
257    }
258    
259    // add points to the fitter
260    for (Int_t i = 0; i < entries; i++){
261       Double_t x[1000];
262       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
263       fitter->AddPoint(x, values[dim][i], 1);
264    }
265
266    fitter->Eval();
267    fitter->GetParameters(fitParam);
268    fitter->GetCovarianceMatrix(covMatrix);
269    chi2 = fitter->GetChisquare();
270    chi2 = chi2;
271    
272    TString *preturnFormula = new TString(Form("( %e+",fitParam[0])), &returnFormula = *preturnFormula;
273    
274    for (Int_t iparam = 0; iparam < dim; iparam++) {
275      returnFormula.Append(Form("%s*(%e)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
276      if (iparam < dim-1) returnFormula.Append("+");
277    }
278    returnFormula.Append(" )");
279    delete formulaTokens;
280    delete fitter;
281    delete[] values;
282    return preturnFormula;
283 }
284
285 //_____________________________________________________________________________
286 Double_t AliBaseCalibViewer::GetLTM(Int_t n, Double_t *array, Double_t *sigma, Double_t fraction){
287    //
288    //  returns the LTM and sigma
289    //
290    Double_t *ddata = new Double_t[n];
291    Double_t mean = 0, lsigma = 0;
292    UInt_t nPoints = 0;
293    for (UInt_t i = 0; i < (UInt_t)n; i++) {
294          ddata[nPoints]= array[nPoints];
295          nPoints++;
296    }
297    Int_t hh = TMath::Min(TMath::Nint(fraction * nPoints), Int_t(n));
298    AliMathBase::EvaluateUni(nPoints, ddata, mean, lsigma, hh);
299    if (sigma) *sigma = lsigma;
300    delete [] ddata;
301    return mean;
302 }
303
304 //_____________________________________________________________________________
305 Int_t AliBaseCalibViewer::GetBin(Float_t value, Int_t nbins, Double_t binLow, Double_t binUp){
306    // Returns the 'bin' for 'value'
307    // The interval between 'binLow' and 'binUp' is divided into 'nbins' equidistant bins
308    // avoid index out of bounds error: 'if (bin < binLow) bin = binLow' and vice versa
309    /* Begin_Latex
310          GetBin(value) = #frac{nbins - 1}{binUp - binLow} #upoint (value - binLow) +1
311       End_Latex
312    */
313    
314    Int_t bin =  TMath::Nint( (Float_t)(value - binLow) / (Float_t)(binUp - binLow) * (nbins-1) ) + 1;
315    // avoid index out of bounds:   
316    if (value < binLow) bin = 0;
317    if (value > binUp)  bin = nbins + 1;
318    return bin;
319    
320 }
321
322 //_____________________________________________________________________________
323 TH1F* AliBaseCalibViewer::SigmaCut(TH1F *histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, 
324                                    Float_t sigmaStep, Bool_t pm) {
325    //
326    // Creates a cumulative histogram Begin_Latex S(t, #mu, #sigma) End_Latex, where you can see, how much of the data are inside sigma-intervals around the mean value
327    // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'histogram'
328    // 'mean' and 'sigma' are Begin_Latex #mu End_Latex and  Begin_Latex #sigma End_Latex of the distribution in 'histogram', to be specified by the user
329    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
330    // sigmaStep: the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
331    // pm: Decide weather Begin_Latex t > 0 End_Latex (first case) or Begin_Latex t End_Latex arbitrary (secound case)
332    // The actual work is done on the array.
333    /* Begin_Latex 
334          f(x, #mu, #sigma)     #Rightarrow       S(t, #mu, #sigma) = #frac{#int_{#mu}^{#mu + t #sigma} f(x, #mu, #sigma) dx + #int_{#mu}^{#mu - t #sigma} f(x, #mu, #sigma) dx }{ #int_{-#infty}^{+#infty} f(x, #mu, #sigma) dx } ,    for  t > 0    
335          or      
336          f(x, #mu, #sigma)     #Rightarrow       S(t, #mu, #sigma) = #frac{#int_{#mu}^{#mu + t #sigma} f(x, #mu, #sigma) dx}{ #int_{-#infty}^{+#infty} f(x, #mu, #sigma) dx }
337       End_Latex  
338       begin_macro(source)
339       {
340          Float_t mean = 0;
341          Float_t sigma = 1.5;
342          Float_t sigmaMax = 4;
343          gROOT->SetStyle("Plain");
344          TH1F *distribution = new TH1F("Distribution1", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
345          TRandom rand(23);
346          for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
347          Float_t *ar = distribution->GetArray();
348          
349          TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas_SigmaCut", "", 350, 350);
350          macro_example_canvas->Divide(0,3);
351          TVirtualPad *pad1 = macro_example_canvas->cd(1);
352          pad1->SetGridy();
353          pad1->SetGridx();
354          distribution->Draw();
355          TVirtualPad *pad2 = macro_example_canvas->cd(2);
356          pad2->SetGridy();
357          pad2->SetGridx();
358          
359          TH1F *shist = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax);
360          shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
361          shist->Draw();  
362          TVirtualPad *pad3 = macro_example_canvas->cd(3);
363          pad3->SetGridy();
364          pad3->SetGridx();
365          TH1F *shistPM = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax, -1, kTRUE);
366          shistPM->Draw();   
367          return macro_example_canvas;
368       }  
369       end_macro
370    */ 
371    
372    Float_t *array = histogram->GetArray();
373    Int_t    nbins = histogram->GetXaxis()->GetNbins();
374    Float_t binLow = histogram->GetXaxis()->GetXmin();
375    Float_t binUp  = histogram->GetXaxis()->GetXmax();
376    return AliBaseCalibViewer::SigmaCut(nbins, array, mean, sigma, nbins, binLow, binUp, sigmaMax, sigmaStep, pm);
377 }   
378    
379 //_____________________________________________________________________________
380 TH1F* AliBaseCalibViewer::SigmaCut(Int_t n, Float_t *array, Float_t mean, Float_t sigma, Int_t nbins, Float_t binLow, Float_t binUp, Float_t sigmaMax, Float_t sigmaStep, Bool_t pm){
381    //
382    // Creates a histogram Begin_Latex S(t, #mu, #sigma) End_Latex, where you can see, how much of the data are inside sigma-intervals around the mean value
383    // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'array', 'n' specifies the length of the array
384    // 'mean' and 'sigma' are Begin_Latex #mu End_Latex and  Begin_Latex #sigma End_Latex of the distribution in 'array', to be specified by the user
385    // 'nbins': number of bins, 'binLow': first bin, 'binUp': last bin
386    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
387    // sigmaStep: the binsize of the generated histogram
388    // Here the actual work is done.
389    
390    if (TMath::Abs(sigma) < 1.e-10) return 0;
391    Float_t binWidth = (binUp-binLow)/(nbins - 1);
392    if (sigmaStep <= 0) sigmaStep = binWidth;
393    Int_t kbins = (Int_t)(sigmaMax * sigma / sigmaStep) + 1; // + 1  due to overflow bin in histograms
394    if (pm) kbins = 2 * (Int_t)(sigmaMax * sigma / sigmaStep) + 1;
395    Float_t kbinLow = !pm ? 0 : -sigmaMax;
396    Float_t kbinUp  = sigmaMax;
397    TH1F *hist = new TH1F("sigmaCutHisto","Cumulative; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp); 
398    hist->SetDirectory(0);
399    hist->Reset();
400    
401    // calculate normalization
402    Double_t normalization = 0;
403    for (Int_t i = 0; i <= n; i++) {
404         normalization += array[i];
405    }
406    
407    // given units: units from given histogram
408    // sigma units: in units of sigma
409    // iDelta: integrate in interval (mean +- iDelta), given units
410    // x:      ofset from mean for integration, given units
411    // hist:   needs 
412    
413    // fill histogram
414    for (Float_t iDelta = 0; iDelta <= sigmaMax * sigma; iDelta += sigmaStep) {
415       // integrate array
416       Double_t valueP = array[GetBin(mean, nbins, binLow, binUp)];
417       Double_t valueM = array[GetBin(mean-binWidth, nbins, binLow, binUp)];
418       // add bin of mean value only once to the histogram
419       for (Float_t x = binWidth; x <= iDelta; x += binWidth) {
420          valueP += (mean + x <= binUp)  ? array[GetBin(mean + x, nbins, binLow, binUp)] : 0;
421          valueM += (mean-binWidth - x >= binLow) ? array[GetBin(mean-binWidth - x, nbins, binLow, binUp)] : 0; 
422       }
423
424       if (valueP / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", valueP, normalization);
425       if (valueP / normalization > 100) return hist;
426       if (valueM / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", valueM, normalization);
427       if (valueM / normalization > 100) return hist;
428       valueP = (valueP / normalization);
429       valueM = (valueM / normalization);
430       if (pm) {
431          Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
432          hist->SetBinContent(bin, valueP);
433          bin = GetBin(-iDelta/sigma, kbins, kbinLow, kbinUp);
434          hist->SetBinContent(bin, valueM);
435       }
436       else { // if (!pm)
437          Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
438          hist->SetBinContent(bin, valueP + valueM);
439       }
440    }
441    if (!pm) hist->SetMaximum(1.2);
442    return hist;
443 }
444
445 //_____________________________________________________________________________
446 TH1F* AliBaseCalibViewer::SigmaCut(Int_t n, Double_t *array, Double_t mean, Double_t sigma, Int_t nbins, Double_t *xbins, Double_t sigmaMax){
447    // 
448    // SigmaCut for variable binsize
449    // NOT YET IMPLEMENTED !!!
450    // 
451    printf("SigmaCut with variable binsize, Not yet implemented\n");
452    // avoid compiler warnings:
453    n=n;
454    mean=mean;
455    sigma=sigma;
456    nbins=nbins;
457    sigmaMax=sigmaMax;
458    array=array;
459    xbins=xbins;
460    
461    return 0;
462 }
463
464 //_____________________________________________________________________________
465 Int_t  AliBaseCalibViewer::DrawHisto1D(const Char_t* drawCommand, const Char_t* sector, const Char_t* cuts, 
466                                        const Char_t *sigmas, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM) const {
467    // 
468    // Easy drawing of data, in principle the same as EasyDraw1D
469    // Difference: A line for the mean / median / LTM is drawn 
470    // in 'sigmas' you can specify in which distance to the mean/median/LTM you want to see a line in sigma-units, separated by ';'
471    // example: sigmas = "2; 4; 6;"  at Begin_Latex 2 #sigma End_Latex, Begin_Latex 4 #sigma End_Latex and Begin_Latex 6 #sigma End_Latex  a line is drawn.
472    // "plotMean", "plotMedian" and "plotLTM": what kind of lines do you want to see?
473    // 
474    Int_t oldOptStat = gStyle->GetOptStat();
475    gStyle->SetOptStat(0000000);
476    Double_t ltmFraction = 0.8;
477    
478    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");  
479    TVectorF nsigma(sigmasTokens->GetEntriesFast());
480    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
481       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
482       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
483       nsigma[i] = sig;
484    }
485    
486    TString drawStr(drawCommand);
487    Bool_t dangerousToDraw = drawStr.Contains(":") || drawStr.Contains(">>");
488    if (dangerousToDraw) {
489       Warning("DrawHisto1D", "The draw string must not contain ':' or '>>'.");
490       return -1;
491    }
492    drawStr += " >> tempHist";
493    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts);
494    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
495    // FIXME is this histogram deleted automatically?
496    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
497    
498    Double_t mean = TMath::Mean(entries, values);
499    Double_t median = TMath::Median(entries, values);
500    Double_t sigma = TMath::RMS(entries, values);
501    Double_t maxY = htemp->GetMaximum();
502    
503    Char_t c[500];
504    TLegend * legend = new TLegend(.7,.7, .99, .99, "Statistical information");
505
506    if (plotMean) {
507       // draw Mean
508       TLine* line = new TLine(mean, 0, mean, maxY);
509       line->SetLineColor(kRed);
510       line->SetLineWidth(2);
511       line->SetLineStyle(1);
512       line->Draw();
513       sprintf(c, "Mean: %f", mean);
514       legend->AddEntry(line, c, "l");
515       // draw sigma lines
516       for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
517          TLine* linePlusSigma = new TLine(mean + nsigma[i] * sigma, 0, mean + nsigma[i] * sigma, maxY);
518          linePlusSigma->SetLineColor(kRed);
519          linePlusSigma->SetLineStyle(2 + i);
520          linePlusSigma->Draw();
521          TLine* lineMinusSigma = new TLine(mean - nsigma[i] * sigma, 0, mean - nsigma[i] * sigma, maxY);
522          lineMinusSigma->SetLineColor(kRed);
523          lineMinusSigma->SetLineStyle(2 + i);
524          lineMinusSigma->Draw();
525          sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
526          legend->AddEntry(lineMinusSigma, c, "l");
527       }
528    }
529    if (plotMedian) {
530       // draw median
531       TLine* line = new TLine(median, 0, median, maxY);
532       line->SetLineColor(kBlue);
533       line->SetLineWidth(2);
534       line->SetLineStyle(1);
535       line->Draw();
536       sprintf(c, "Median: %f", median);
537       legend->AddEntry(line, c, "l");
538       // draw sigma lines
539       for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
540          TLine* linePlusSigma = new TLine(median + nsigma[i] * sigma, 0, median + nsigma[i]*sigma, maxY);
541          linePlusSigma->SetLineColor(kBlue);
542          linePlusSigma->SetLineStyle(2 + i);
543          linePlusSigma->Draw();
544          TLine* lineMinusSigma = new TLine(median - nsigma[i] * sigma, 0, median - nsigma[i]*sigma, maxY);
545          lineMinusSigma->SetLineColor(kBlue);
546          lineMinusSigma->SetLineStyle(2 + i);
547          lineMinusSigma->Draw();
548          sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
549          legend->AddEntry(lineMinusSigma, c, "l");
550       }
551    }
552    if (plotLTM) {
553       // draw LTM
554       Double_t ltmRms = 0;
555       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
556       TLine* line = new TLine(ltm, 0, ltm, maxY);
557       //fListOfObjectsToBeDeleted->Add(line);
558       line->SetLineColor(kGreen+2);
559       line->SetLineWidth(2);
560       line->SetLineStyle(1);
561       line->Draw();
562       sprintf(c, "LTM: %f", ltm);
563       legend->AddEntry(line, c, "l");
564       // draw sigma lines
565       for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
566          TLine* linePlusSigma = new TLine(ltm + nsigma[i] * ltmRms, 0, ltm + nsigma[i] * ltmRms, maxY);
567          //fListOfObjectsToBeDeleted->Add(linePlusSigma);
568          linePlusSigma->SetLineColor(kGreen+2);
569          linePlusSigma->SetLineStyle(2+i);
570          linePlusSigma->Draw();
571    
572          TLine* lineMinusSigma = new TLine(ltm - nsigma[i] * ltmRms, 0, ltm - nsigma[i] * ltmRms, maxY);
573          //fListOfObjectsToBeDeleted->Add(lineMinusSigma);
574          lineMinusSigma->SetLineColor(kGreen+2);
575          lineMinusSigma->SetLineStyle(2+i);
576          lineMinusSigma->Draw();
577          sprintf(c, "%i #sigma = %f", (Int_t)(nsigma[i]), (Float_t)(nsigma[i] * ltmRms));
578          legend->AddEntry(lineMinusSigma, c, "l");
579       }
580    }
581    if (!plotMean && !plotMedian && !plotLTM) return -1;
582    legend->Draw();
583    gStyle->SetOptStat(oldOptStat);
584    return 1;
585 }
586
587 //_____________________________________________________________________________
588 Int_t AliBaseCalibViewer::SigmaCut(const Char_t* drawCommand, const Char_t* sector, const Char_t* cuts, 
589                                    Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, Bool_t pm, 
590                                    const Char_t *sigmas, Float_t sigmaStep) const {
591    //
592    // Creates a histogram, where you can see, how much of the data are inside sigma-intervals 
593    // around the mean/median/LTM
594    // with drawCommand, sector and cuts you specify your input data, see EasyDraw
595    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma)
596    // sigmaStep: the binsize of the generated histogram
597    // plotMean/plotMedian/plotLTM: specifies where to put the center
598    //
599   
600    Double_t ltmFraction = 0.8;
601    
602    TString drawStr(drawCommand);
603    Bool_t dangerousToDraw = drawStr.Contains(":") || drawStr.Contains(">>");
604    if (dangerousToDraw) {
605       Warning("SigmaCut", "The draw string must not contain ':' or '>>'.");
606       return -1;
607    }
608    drawStr += " >> tempHist";
609    
610    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
611    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
612    // FIXME is this histogram deleted automatically?
613    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
614    
615    Double_t mean = TMath::Mean(entries, values);
616    Double_t median = TMath::Median(entries, values);
617    Double_t sigma = TMath::RMS(entries, values);
618    
619    TLegend * legend = new TLegend(.7,.7, .99, .99, "Cumulative");
620    //fListOfObjectsToBeDeleted->Add(legend);
621    TH1F *cutHistoMean = 0;
622    TH1F *cutHistoMedian = 0;
623    TH1F *cutHistoLTM = 0;
624    
625    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");  
626    TVectorF nsigma(sigmasTokens->GetEntriesFast());
627    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
628       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
629       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
630       nsigma[i] = sig;
631    }
632   
633    if (plotMean) {
634       cutHistoMean = SigmaCut(htemp, mean, sigma, sigmaMax, sigmaStep, pm);
635       if (cutHistoMean) {
636          //fListOfObjectsToBeDeleted->Add(cutHistoMean);
637          cutHistoMean->SetLineColor(kRed);
638          legend->AddEntry(cutHistoMean, "Mean", "l");
639          cutHistoMean->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
640          cutHistoMean->Draw();
641          DrawLines(cutHistoMean, nsigma, legend, kRed, pm);
642       } // if (cutHistoMean)
643        
644    }
645    if (plotMedian) {
646       cutHistoMedian = SigmaCut(htemp, median, sigma, sigmaMax, sigmaStep, pm);
647       if (cutHistoMedian) {
648          //fListOfObjectsToBeDeleted->Add(cutHistoMedian);
649          cutHistoMedian->SetLineColor(kBlue);
650          legend->AddEntry(cutHistoMedian, "Median", "l");
651          cutHistoMedian->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
652          if (plotMean && cutHistoMean) cutHistoMedian->Draw("same");
653             else cutHistoMedian->Draw();
654          DrawLines(cutHistoMedian, nsigma, legend, kBlue, pm);
655       }  // if (cutHistoMedian)
656    }
657    if (plotLTM) {
658       Double_t ltmRms = 0;
659       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
660       cutHistoLTM = SigmaCut(htemp, ltm, ltmRms, sigmaMax, sigmaStep, pm);
661       if (cutHistoLTM) {
662          //fListOfObjectsToBeDeleted->Add(cutHistoLTM);
663          cutHistoLTM->SetLineColor(kGreen+2);
664          legend->AddEntry(cutHistoLTM, "LTM", "l");
665          cutHistoLTM->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
666          if ((plotMean && cutHistoMean) || (plotMedian && cutHistoMedian)) cutHistoLTM->Draw("same");
667             else cutHistoLTM->Draw();
668          DrawLines(cutHistoLTM, nsigma, legend, kGreen+2, pm);
669       }
670    }
671    if (!plotMean && !plotMedian && !plotLTM) return -1;
672    legend->Draw();
673    return 1;
674 }
675
676 //_____________________________________________________________________________
677 Int_t AliBaseCalibViewer::Integrate(const Char_t* drawCommand, const Char_t* sector, const Char_t* cuts, 
678                                     Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, 
679                                     const Char_t *sigmas, Float_t sigmaStep) const {
680    //
681    // Creates an integrated histogram Begin_Latex S(t, #mu, #sigma) End_Latex, out of the input distribution distribution Begin_Latex f(x, #mu, #sigma) End_Latex, given in "histogram"   
682    // "mean" and "sigma" are Begin_Latex #mu End_Latex and  Begin_Latex #sigma End_Latex of the distribution in "histogram", to be specified by the user
683    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate 
684    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
685    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
686    // The actual work is done on the array.
687    /* Begin_Latex 
688          f(x, #mu, #sigma)     #Rightarrow       S(t, #mu, #sigma) = #frac{#int_{-#infty}^{#mu + t #sigma} f(x, #mu, #sigma) dx}{ #int_{-#infty}^{+#infty} f(x, #mu, #sigma) dx }
689       End_Latex  
690    */
691    
692    Double_t ltmFraction = 0.8;
693    // avoid compiler warnings:
694    sigmaMax = sigmaMax;
695    sigmaStep = sigmaStep;
696    
697    TString drawStr(drawCommand);
698    Bool_t dangerousToDraw = drawStr.Contains(":") || drawStr.Contains(">>");
699    if (dangerousToDraw) {
700       Warning("Integrate", "The draw string must not contain ':' or '>>'.");
701       return -1;
702    }
703    drawStr += " >> tempHist";
704    
705    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
706    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
707    TGraph *integralGraphMean   = 0;
708    TGraph *integralGraphMedian = 0;
709    TGraph *integralGraphLTM    = 0;
710    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
711    Int_t    *index  = new Int_t[entries];
712    Float_t  *xarray = new Float_t[entries];
713    Float_t  *yarray = new Float_t[entries];
714    TMath::Sort(entries, values, index, kFALSE);
715    
716    Double_t mean = TMath::Mean(entries, values);
717    Double_t median = TMath::Median(entries, values);
718    Double_t sigma = TMath::RMS(entries, values);
719    
720    // parse sigmas string
721    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");  
722    TVectorF nsigma(sigmasTokens->GetEntriesFast());
723    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
724       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
725       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
726       nsigma[i] = sig;
727    }
728   
729    TLegend * legend = new TLegend(.7,.7, .99, .99, "Integrated histogram");
730    //fListOfObjectsToBeDeleted->Add(legend);
731   
732    if (plotMean) {
733       for (Int_t i = 0; i < entries; i++) {
734          xarray[i] = (values[index[i]] - mean) / sigma; 
735          yarray[i] = float(i) / float(entries);
736       }
737       integralGraphMean = new TGraph(entries, xarray, yarray);
738       if (integralGraphMean) {
739          //fListOfObjectsToBeDeleted->Add(integralGraphMean);
740          integralGraphMean->SetLineColor(kRed);
741          legend->AddEntry(integralGraphMean, "Mean", "l");
742          integralGraphMean->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
743          integralGraphMean->Draw("alu");
744          DrawLines(integralGraphMean, nsigma, legend, kRed, kTRUE);
745       }
746    }
747    if (plotMedian) {
748       for (Int_t i = 0; i < entries; i++) {
749          xarray[i] = (values[index[i]] - median) / sigma; 
750          yarray[i] = float(i) / float(entries);
751       }
752       integralGraphMedian = new TGraph(entries, xarray, yarray);
753       if (integralGraphMedian) {
754          //fListOfObjectsToBeDeleted->Add(integralGraphMedian);
755          integralGraphMedian->SetLineColor(kBlue);
756          legend->AddEntry(integralGraphMedian, "Median", "l");
757          integralGraphMedian->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
758          if (plotMean && integralGraphMean) integralGraphMedian->Draw("samelu");
759             else integralGraphMedian->Draw("alu");
760          DrawLines(integralGraphMedian, nsigma, legend, kBlue, kTRUE);
761       }
762    }
763    if (plotLTM) {
764       Double_t ltmRms = 0;
765       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
766       for (Int_t i = 0; i < entries; i++) {
767          xarray[i] = (values[index[i]] - ltm) / ltmRms; 
768          yarray[i] = float(i) / float(entries);
769       }
770       integralGraphLTM = new TGraph(entries, xarray, yarray);
771       if (integralGraphLTM) {
772          //fListOfObjectsToBeDeleted->Add(integralGraphLTM);
773          integralGraphLTM->SetLineColor(kGreen+2);
774          legend->AddEntry(integralGraphLTM, "LTM", "l");
775          integralGraphLTM->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
776          if ((plotMean && integralGraphMean) || (plotMedian && integralGraphMedian)) integralGraphLTM->Draw("samelu");
777             else integralGraphLTM->Draw("alu");
778          DrawLines(integralGraphLTM, nsigma, legend, kGreen+2, kTRUE);
779       }
780    }
781    if (!plotMean && !plotMedian && !plotLTM) return -1;
782    legend->Draw();
783    return entries;
784 }
785
786 //_____________________________________________________________________________
787 TH1F* AliBaseCalibViewer::Integrate(TH1F *histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
788    //
789    // Creates an integrated histogram Begin_Latex S(t, #mu, #sigma) End_Latex, out of the input distribution distribution Begin_Latex f(x, #mu, #sigma) End_Latex, given in "histogram"   
790    // "mean" and "sigma" are Begin_Latex #mu End_Latex and  Begin_Latex #sigma End_Latex of the distribution in "histogram", to be specified by the user
791    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate 
792    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
793    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
794    // The actual work is done on the array.
795    /* Begin_Latex 
796          f(x, #mu, #sigma)     #Rightarrow       S(t, #mu, #sigma) = #frac{#int_{-#infty}^{#mu + t #sigma} f(x, #mu, #sigma) dx}{ #int_{-#infty}^{+#infty} f(x, #mu, #sigma) dx }
797       End_Latex  
798       begin_macro(source)
799       {
800          Float_t mean = 0;
801          Float_t sigma = 1.5;
802          Float_t sigmaMax = 4;
803          gROOT->SetStyle("Plain");
804          TH1F *distribution = new TH1F("Distribution2", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
805          TRandom rand(23);
806          for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
807          Float_t *ar = distribution->GetArray();
808          
809          TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas_Integrate", "", 350, 350);
810          macro_example_canvas->Divide(0,2);
811          TVirtualPad *pad1 = macro_example_canvas->cd(1);
812          pad1->SetGridy();
813          pad1->SetGridx();
814          distribution->Draw();
815          TVirtualPad *pad2 = macro_example_canvas->cd(2);
816          pad2->SetGridy();
817          pad2->SetGridx();
818          TH1F *shist = AliTPCCalibViewer::Integrate(distribution, mean, sigma, sigmaMax);
819          shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
820          shist->Draw();  
821          
822          return macro_example_canvas_Integrate;
823       }  
824       end_macro
825    */ 
826
827    
828    Float_t *array = histogram->GetArray();
829    Int_t    nbins = histogram->GetXaxis()->GetNbins();
830    Float_t binLow = histogram->GetXaxis()->GetXmin();
831    Float_t binUp  = histogram->GetXaxis()->GetXmax();
832    return Integrate(nbins, array, nbins, binLow, binUp, mean, sigma, sigmaMax, sigmaStep);
833 }
834
835 //_____________________________________________________________________________
836 TH1F* AliBaseCalibViewer::Integrate(Int_t n, Float_t *array, Int_t nbins, Float_t binLow, Float_t binUp, 
837                                     Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
838    // Creates an integrated histogram Begin_Latex S(t, #mu, #sigma) End_Latex, out of the input distribution distribution Begin_Latex f(x, #mu, #sigma) End_Latex, given in "histogram"   
839    // "mean" and "sigma" are Begin_Latex #mu End_Latex and  Begin_Latex #sigma End_Latex of the distribution in "histogram", to be specified by the user
840    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate 
841    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
842    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
843    // Here the actual work is done.
844       
845    Bool_t givenUnits = kTRUE;
846    if (TMath::Abs(sigma) < 1.e-10 && TMath::Abs(sigmaMax) < 1.e-10) givenUnits = kFALSE;
847    if (givenUnits) {
848       sigma = 1;
849       sigmaMax = (binUp - binLow) / 2.;
850    }
851    
852    Float_t binWidth = (binUp-binLow)/(nbins - 1);
853    if (sigmaStep <= 0) sigmaStep = binWidth;
854    Int_t kbins =  (Int_t)(sigmaMax * sigma / sigmaStep) + 1;  // + 1  due to overflow bin in histograms
855    Float_t kbinLow = givenUnits ? binLow : -sigmaMax;
856    Float_t kbinUp  = givenUnits ? binUp  : sigmaMax;
857    TH1F *hist = 0; 
858    if (givenUnits)  hist = new TH1F("integratedHisto","Integrated Histogram; Given x; Fraction of included data", kbins, kbinLow, kbinUp); 
859    if (!givenUnits) hist = new TH1F("integratedHisto","Integrated Histogram; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp); 
860    hist->SetDirectory(0);
861    hist->Reset();
862    
863    // calculate normalization
864  //  printf("calculating normalization, integrating from bin 1 to %i \n", n);
865    Double_t normalization = 0;
866    for (Int_t i = 1; i <= n; i++) {
867         normalization += array[i];
868    }
869  //  printf("normalization: %f \n", normalization);
870    
871    // given units: units from given histogram
872    // sigma units: in units of sigma
873    // iDelta: integrate in interval (mean +- iDelta), given units
874    // x:      ofset from mean for integration, given units
875    // hist:   needs 
876    
877    // fill histogram
878    for (Float_t iDelta = mean - sigmaMax * sigma; iDelta <= mean + sigmaMax * sigma; iDelta += sigmaStep) {
879       // integrate array
880       Double_t value = 0;
881       for (Float_t x = mean - sigmaMax * sigma; x <= iDelta; x += binWidth) {
882          value += (x <= binUp && x >= binLow)  ? array[GetBin(x, nbins, binLow, binUp)] : 0;
883       }
884       if (value / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", value, normalization);
885       if (value / normalization > 100) return hist;
886       Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
887     //  printf("first integration bin: %i, last integration bin: %i \n", GetBin(mean - sigmaMax * sigma, nbins, binLow, binUp), GetBin(iDelta, nbins, binLow, binUp));
888     //  printf("value: %f, normalization: %f, normalized value: %f, iDelta: %f, Bin: %i \n", value, normalization, value/normalization, iDelta, bin);
889       value = (value / normalization);
890       hist->SetBinContent(bin, value);
891    }
892    return hist;
893 }
894
895 //_____________________________________________________________________________
896 void AliBaseCalibViewer::DrawLines(TH1F *histogram, TVectorF nsigma, TLegend *legend, Int_t color, Bool_t pm) const {
897    // 
898    // Private function for SigmaCut(...) and Integrate(...)
899    // Draws lines into the given histogram, specified by "nsigma", the lines are addeed to the legend
900    // 
901    
902    // start to draw the lines, loop over requested sigmas
903    Char_t c[500];
904    for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
905       if (!pm) { 
906          Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
907          TLine* lineUp = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
908          //fListOfObjectsToBeDeleted->Add(lineUp);
909          lineUp->SetLineColor(color);
910          lineUp->SetLineStyle(2 + i);
911          lineUp->Draw();
912          TLine* lineLeft = new TLine(nsigma[i], histogram->GetBinContent(bin), 0, histogram->GetBinContent(bin));
913          //fListOfObjectsToBeDeleted->Add(lineLeft);
914          lineLeft->SetLineColor(color);
915          lineLeft->SetLineStyle(2 + i);
916          lineLeft->Draw();
917          sprintf(c, "Fraction(%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
918          legend->AddEntry(lineLeft, c, "l");
919       }
920       else { // if (pm)
921          Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
922          TLine* lineUp1 = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
923          //fListOfObjectsToBeDeleted->Add(lineUp1);
924          lineUp1->SetLineColor(color);
925          lineUp1->SetLineStyle(2 + i);
926          lineUp1->Draw();
927          TLine* lineLeft1 = new TLine(nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
928          //fListOfObjectsToBeDeleted->Add(lineLeft1);
929          lineLeft1->SetLineColor(color);
930          lineLeft1->SetLineStyle(2 + i);
931          lineLeft1->Draw();
932          sprintf(c, "Fraction(+%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
933          legend->AddEntry(lineLeft1, c, "l");
934          bin = histogram->GetXaxis()->FindBin(-nsigma[i]);
935          TLine* lineUp2 = new TLine(-nsigma[i], 0, -nsigma[i], histogram->GetBinContent(bin));
936          //fListOfObjectsToBeDeleted->Add(lineUp2);
937          lineUp2->SetLineColor(color);
938          lineUp2->SetLineStyle(2 + i);
939          lineUp2->Draw();
940          TLine* lineLeft2 = new TLine(-nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
941          //fListOfObjectsToBeDeleted->Add(lineLeft2);
942          lineLeft2->SetLineColor(color);
943          lineLeft2->SetLineStyle(2 + i);
944          lineLeft2->Draw();
945          sprintf(c, "Fraction(-%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
946          legend->AddEntry(lineLeft2, c, "l");
947       }
948    }  // for (Int_t i = 0; i < nsigma.GetNoElements(); i++)   
949 }
950
951 //_____________________________________________________________________________
952 void AliBaseCalibViewer::DrawLines(TGraph *graph, TVectorF nsigma, TLegend *legend, Int_t color, Bool_t pm) const {
953    // 
954    // Private function for SigmaCut(...) and Integrate(...)
955    // Draws lines into the given histogram, specified by "nsigma", the lines are addeed to the legend
956    // 
957    
958    // start to draw the lines, loop over requested sigmas
959    Char_t c[500];
960    for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
961       if (!pm) { 
962          TLine* lineUp = new TLine(nsigma[i], 0, nsigma[i], graph->Eval(nsigma[i]));
963          //fListOfObjectsToBeDeleted->Add(lineUp);
964          lineUp->SetLineColor(color);
965          lineUp->SetLineStyle(2 + i);
966          lineUp->Draw();
967          TLine* lineLeft = new TLine(nsigma[i], graph->Eval(nsigma[i]), 0, graph->Eval(nsigma[i]));
968          //fListOfObjectsToBeDeleted->Add(lineLeft);
969          lineLeft->SetLineColor(color);
970          lineLeft->SetLineStyle(2 + i);
971          lineLeft->Draw();
972          sprintf(c, "Fraction(%f #sigma) = %f",nsigma[i], graph->Eval(nsigma[i]));
973          legend->AddEntry(lineLeft, c, "l");
974       }
975       else { // if (pm)
976          TLine* lineUp1 = new TLine(nsigma[i], 0, nsigma[i], graph->Eval(nsigma[i]));
977          //fListOfObjectsToBeDeleted->Add(lineUp1);
978          lineUp1->SetLineColor(color);
979          lineUp1->SetLineStyle(2 + i);
980          lineUp1->Draw();
981          TLine* lineLeft1 = new TLine(nsigma[i], graph->Eval(nsigma[i]), graph->GetHistogram()->GetXaxis()->GetBinLowEdge(0), graph->Eval(nsigma[i]));
982          //fListOfObjectsToBeDeleted->Add(lineLeft1);
983          lineLeft1->SetLineColor(color);
984          lineLeft1->SetLineStyle(2 + i);
985          lineLeft1->Draw();
986          sprintf(c, "Fraction(+%f #sigma) = %f",nsigma[i], graph->Eval(nsigma[i]));
987          legend->AddEntry(lineLeft1, c, "l");
988          TLine* lineUp2 = new TLine(-nsigma[i], 0, -nsigma[i], graph->Eval(-nsigma[i]));
989          //fListOfObjectsToBeDeleted->Add(lineUp2);
990          lineUp2->SetLineColor(color);
991          lineUp2->SetLineStyle(2 + i);
992          lineUp2->Draw();
993          TLine* lineLeft2 = new TLine(-nsigma[i], graph->Eval(-nsigma[i]), graph->GetHistogram()->GetXaxis()->GetBinLowEdge(0), graph->Eval(-nsigma[i]));
994          //fListOfObjectsToBeDeleted->Add(lineLeft2);
995          lineLeft2->SetLineColor(color);
996          lineLeft2->SetLineStyle(2 + i);
997          lineLeft2->Draw();
998          sprintf(c, "Fraction(-%f #sigma) = %f",nsigma[i], graph->Eval(-nsigma[i]));
999          legend->AddEntry(lineLeft2, c, "l");
1000       }
1001    }  // for (Int_t i = 0; i < nsigma.GetNoElements(); i++)   
1002 }