1b2ce08724c8020247b0ecda1634aa214b0b7387
[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) {
247      delete fitter;
248      return new TString("An ERROR has occured during fitting!");
249    }
250    Double_t **values = new Double_t*[dim+1] ; 
251    
252    for (Int_t i = 0; i < dim + 1; i++){
253       Int_t centries = 0;
254       if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff");
255       else  centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff");
256       
257       if (entries != centries) {
258         delete fitter;
259         delete [] values;
260         return new TString("An ERROR has occured during fitting!");
261       }
262       values[i] = new Double_t[entries];
263       memcpy(values[i],  fTree->GetV1(), entries*sizeof(Double_t)); 
264    }
265    
266    // add points to the fitter
267    for (Int_t i = 0; i < entries; i++){
268       Double_t x[1000];
269       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
270       fitter->AddPoint(x, values[dim][i], 1);
271    }
272
273    fitter->Eval();
274    fitter->GetParameters(fitParam);
275    fitter->GetCovarianceMatrix(covMatrix);
276    chi2 = fitter->GetChisquare();
277    chi2 = chi2;
278    
279    TString *preturnFormula = new TString(Form("( %e+",fitParam[0])), &returnFormula = *preturnFormula;
280    
281    for (Int_t iparam = 0; iparam < dim; iparam++) {
282      returnFormula.Append(Form("%s*(%e)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
283      if (iparam < dim-1) returnFormula.Append("+");
284    }
285    returnFormula.Append(" )");
286    delete formulaTokens;
287    delete fitter;
288    for (Int_t i = 0; i < dim + 1; i++) delete [] values[i];
289    delete[] values;
290    return preturnFormula;
291 }
292
293 //_____________________________________________________________________________
294 Double_t AliBaseCalibViewer::GetLTM(Int_t n, Double_t *array, Double_t *sigma, Double_t fraction){
295    //
296    //  returns the LTM and sigma
297    //
298    Double_t *ddata = new Double_t[n];
299    Double_t mean = 0, lsigma = 0;
300    UInt_t nPoints = 0;
301    for (UInt_t i = 0; i < (UInt_t)n; i++) {
302          ddata[nPoints]= array[nPoints];
303          nPoints++;
304    }
305    Int_t hh = TMath::Min(TMath::Nint(fraction * nPoints), Int_t(n));
306    AliMathBase::EvaluateUni(nPoints, ddata, mean, lsigma, hh);
307    if (sigma) *sigma = lsigma;
308    delete [] ddata;
309    return mean;
310 }
311
312 //_____________________________________________________________________________
313 Int_t AliBaseCalibViewer::GetBin(Float_t value, Int_t nbins, Double_t binLow, Double_t binUp){
314    // Returns the 'bin' for 'value'
315    // The interval between 'binLow' and 'binUp' is divided into 'nbins' equidistant bins
316    // avoid index out of bounds error: 'if (bin < binLow) bin = binLow' and vice versa
317    /* Begin_Latex
318          GetBin(value) = #frac{nbins - 1}{binUp - binLow} #upoint (value - binLow) +1
319       End_Latex
320    */
321    
322    Int_t bin =  TMath::Nint( (Float_t)(value - binLow) / (Float_t)(binUp - binLow) * (nbins-1) ) + 1;
323    // avoid index out of bounds:   
324    if (value < binLow) bin = 0;
325    if (value > binUp)  bin = nbins + 1;
326    return bin;
327    
328 }
329
330 //_____________________________________________________________________________
331 TH1F* AliBaseCalibViewer::SigmaCut(TH1F *histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, 
332                                    Float_t sigmaStep, Bool_t pm) {
333    //
334    // 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
335    // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'histogram'
336    // '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
337    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
338    // sigmaStep: the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
339    // pm: Decide weather Begin_Latex t > 0 End_Latex (first case) or Begin_Latex t End_Latex arbitrary (secound case)
340    // The actual work is done on the array.
341    /* Begin_Latex 
342          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    
343          or      
344          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 }
345       End_Latex  
346       begin_macro(source)
347       {
348          Float_t mean = 0;
349          Float_t sigma = 1.5;
350          Float_t sigmaMax = 4;
351          gROOT->SetStyle("Plain");
352          TH1F *distribution = new TH1F("Distribution1", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
353          TRandom rand(23);
354          for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
355          Float_t *ar = distribution->GetArray();
356          
357          TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas_SigmaCut", "", 350, 350);
358          macro_example_canvas->Divide(0,3);
359          TVirtualPad *pad1 = macro_example_canvas->cd(1);
360          pad1->SetGridy();
361          pad1->SetGridx();
362          distribution->Draw();
363          TVirtualPad *pad2 = macro_example_canvas->cd(2);
364          pad2->SetGridy();
365          pad2->SetGridx();
366          
367          TH1F *shist = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax);
368          shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
369          shist->Draw();  
370          TVirtualPad *pad3 = macro_example_canvas->cd(3);
371          pad3->SetGridy();
372          pad3->SetGridx();
373          TH1F *shistPM = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax, -1, kTRUE);
374          shistPM->Draw();   
375          return macro_example_canvas;
376       }  
377       end_macro
378    */ 
379    
380    Float_t *array = histogram->GetArray();
381    Int_t    nbins = histogram->GetXaxis()->GetNbins();
382    Float_t binLow = histogram->GetXaxis()->GetXmin();
383    Float_t binUp  = histogram->GetXaxis()->GetXmax();
384    return AliBaseCalibViewer::SigmaCut(nbins, array, mean, sigma, nbins, binLow, binUp, sigmaMax, sigmaStep, pm);
385 }   
386    
387 //_____________________________________________________________________________
388 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){
389    //
390    // 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
391    // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'array', 'n' specifies the length of the array
392    // '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
393    // 'nbins': number of bins, 'binLow': first bin, 'binUp': last bin
394    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
395    // sigmaStep: the binsize of the generated histogram
396    // Here the actual work is done.
397    
398    if (TMath::Abs(sigma) < 1.e-10) return 0;
399    Float_t binWidth = (binUp-binLow)/(nbins - 1);
400    if (sigmaStep <= 0) sigmaStep = binWidth;
401    Int_t kbins = (Int_t)(sigmaMax * sigma / sigmaStep) + 1; // + 1  due to overflow bin in histograms
402    if (pm) kbins = 2 * (Int_t)(sigmaMax * sigma / sigmaStep) + 1;
403    Float_t kbinLow = !pm ? 0 : -sigmaMax;
404    Float_t kbinUp  = sigmaMax;
405    TH1F *hist = new TH1F("sigmaCutHisto","Cumulative; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp); 
406    hist->SetDirectory(0);
407    hist->Reset();
408    
409    // calculate normalization
410    Double_t normalization = 0;
411    for (Int_t i = 0; i <= n; i++) {
412         normalization += array[i];
413    }
414    
415    // given units: units from given histogram
416    // sigma units: in units of sigma
417    // iDelta: integrate in interval (mean +- iDelta), given units
418    // x:      ofset from mean for integration, given units
419    // hist:   needs 
420    
421    // fill histogram
422    for (Float_t iDelta = 0; iDelta <= sigmaMax * sigma; iDelta += sigmaStep) {
423       // integrate array
424       Double_t valueP = array[GetBin(mean, nbins, binLow, binUp)];
425       Double_t valueM = array[GetBin(mean-binWidth, nbins, binLow, binUp)];
426       // add bin of mean value only once to the histogram
427       for (Float_t x = binWidth; x <= iDelta; x += binWidth) {
428          valueP += (mean + x <= binUp)  ? array[GetBin(mean + x, nbins, binLow, binUp)] : 0;
429          valueM += (mean-binWidth - x >= binLow) ? array[GetBin(mean-binWidth - x, nbins, binLow, binUp)] : 0; 
430       }
431
432       if (valueP / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", valueP, normalization);
433       if (valueP / normalization > 100) return hist;
434       if (valueM / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", valueM, normalization);
435       if (valueM / normalization > 100) return hist;
436       valueP = (valueP / normalization);
437       valueM = (valueM / normalization);
438       if (pm) {
439          Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
440          hist->SetBinContent(bin, valueP);
441          bin = GetBin(-iDelta/sigma, kbins, kbinLow, kbinUp);
442          hist->SetBinContent(bin, valueM);
443       }
444       else { // if (!pm)
445          Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
446          hist->SetBinContent(bin, valueP + valueM);
447       }
448    }
449    if (!pm) hist->SetMaximum(1.2);
450    return hist;
451 }
452
453 //_____________________________________________________________________________
454 TH1F* AliBaseCalibViewer::SigmaCut(Int_t n, Double_t *array, Double_t mean, Double_t sigma, Int_t nbins, Double_t *xbins, Double_t sigmaMax){
455    // 
456    // SigmaCut for variable binsize
457    // NOT YET IMPLEMENTED !!!
458    // 
459    printf("SigmaCut with variable binsize, Not yet implemented\n");
460    // avoid compiler warnings:
461    n=n;
462    mean=mean;
463    sigma=sigma;
464    nbins=nbins;
465    sigmaMax=sigmaMax;
466    array=array;
467    xbins=xbins;
468    
469    return 0;
470 }
471
472 //_____________________________________________________________________________
473 Int_t  AliBaseCalibViewer::DrawHisto1D(const Char_t* drawCommand, const Char_t* sector, const Char_t* cuts, 
474                                        const Char_t *sigmas, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM) const {
475    // 
476    // Easy drawing of data, in principle the same as EasyDraw1D
477    // Difference: A line for the mean / median / LTM is drawn 
478    // in 'sigmas' you can specify in which distance to the mean/median/LTM you want to see a line in sigma-units, separated by ';'
479    // 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.
480    // "plotMean", "plotMedian" and "plotLTM": what kind of lines do you want to see?
481    // 
482    Int_t oldOptStat = gStyle->GetOptStat();
483    gStyle->SetOptStat(0000000);
484    Double_t ltmFraction = 0.8;
485    
486    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");  
487    TVectorF nsigma(sigmasTokens->GetEntriesFast());
488    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
489       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
490       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
491       nsigma[i] = sig;
492    }
493    
494    TString drawStr(drawCommand);
495    Bool_t dangerousToDraw = drawStr.Contains(":") || drawStr.Contains(">>");
496    if (dangerousToDraw) {
497       Warning("DrawHisto1D", "The draw string must not contain ':' or '>>'.");
498       return -1;
499    }
500    drawStr += " >> tempHist";
501    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts);
502    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
503    // FIXME is this histogram deleted automatically?
504    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
505    
506    Double_t mean = TMath::Mean(entries, values);
507    Double_t median = TMath::Median(entries, values);
508    Double_t sigma = TMath::RMS(entries, values);
509    Double_t maxY = htemp->GetMaximum();
510    
511    Char_t c[500];
512    TLegend * legend = new TLegend(.7,.7, .99, .99, "Statistical information");
513
514    if (plotMean) {
515       // draw Mean
516       TLine* line = new TLine(mean, 0, mean, maxY);
517       line->SetLineColor(kRed);
518       line->SetLineWidth(2);
519       line->SetLineStyle(1);
520       line->Draw();
521       sprintf(c, "Mean: %f", mean);
522       legend->AddEntry(line, c, "l");
523       // draw sigma lines
524       for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
525          TLine* linePlusSigma = new TLine(mean + nsigma[i] * sigma, 0, mean + nsigma[i] * sigma, maxY);
526          linePlusSigma->SetLineColor(kRed);
527          linePlusSigma->SetLineStyle(2 + i);
528          linePlusSigma->Draw();
529          TLine* lineMinusSigma = new TLine(mean - nsigma[i] * sigma, 0, mean - nsigma[i] * sigma, maxY);
530          lineMinusSigma->SetLineColor(kRed);
531          lineMinusSigma->SetLineStyle(2 + i);
532          lineMinusSigma->Draw();
533          sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
534          legend->AddEntry(lineMinusSigma, c, "l");
535       }
536    }
537    if (plotMedian) {
538       // draw median
539       TLine* line = new TLine(median, 0, median, maxY);
540       line->SetLineColor(kBlue);
541       line->SetLineWidth(2);
542       line->SetLineStyle(1);
543       line->Draw();
544       sprintf(c, "Median: %f", median);
545       legend->AddEntry(line, c, "l");
546       // draw sigma lines
547       for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
548          TLine* linePlusSigma = new TLine(median + nsigma[i] * sigma, 0, median + nsigma[i]*sigma, maxY);
549          linePlusSigma->SetLineColor(kBlue);
550          linePlusSigma->SetLineStyle(2 + i);
551          linePlusSigma->Draw();
552          TLine* lineMinusSigma = new TLine(median - nsigma[i] * sigma, 0, median - nsigma[i]*sigma, maxY);
553          lineMinusSigma->SetLineColor(kBlue);
554          lineMinusSigma->SetLineStyle(2 + i);
555          lineMinusSigma->Draw();
556          sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
557          legend->AddEntry(lineMinusSigma, c, "l");
558       }
559    }
560    if (plotLTM) {
561       // draw LTM
562       Double_t ltmRms = 0;
563       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
564       TLine* line = new TLine(ltm, 0, ltm, maxY);
565       //fListOfObjectsToBeDeleted->Add(line);
566       line->SetLineColor(kGreen+2);
567       line->SetLineWidth(2);
568       line->SetLineStyle(1);
569       line->Draw();
570       sprintf(c, "LTM: %f", ltm);
571       legend->AddEntry(line, c, "l");
572       // draw sigma lines
573       for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
574          TLine* linePlusSigma = new TLine(ltm + nsigma[i] * ltmRms, 0, ltm + nsigma[i] * ltmRms, maxY);
575          //fListOfObjectsToBeDeleted->Add(linePlusSigma);
576          linePlusSigma->SetLineColor(kGreen+2);
577          linePlusSigma->SetLineStyle(2+i);
578          linePlusSigma->Draw();
579    
580          TLine* lineMinusSigma = new TLine(ltm - nsigma[i] * ltmRms, 0, ltm - nsigma[i] * ltmRms, maxY);
581          //fListOfObjectsToBeDeleted->Add(lineMinusSigma);
582          lineMinusSigma->SetLineColor(kGreen+2);
583          lineMinusSigma->SetLineStyle(2+i);
584          lineMinusSigma->Draw();
585          sprintf(c, "%i #sigma = %f", (Int_t)(nsigma[i]), (Float_t)(nsigma[i] * ltmRms));
586          legend->AddEntry(lineMinusSigma, c, "l");
587       }
588    }
589    if (!plotMean && !plotMedian && !plotLTM) return -1;
590    legend->Draw();
591    gStyle->SetOptStat(oldOptStat);
592    return 1;
593 }
594
595 //_____________________________________________________________________________
596 Int_t AliBaseCalibViewer::SigmaCut(const Char_t* drawCommand, const Char_t* sector, const Char_t* cuts, 
597                                    Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, Bool_t pm, 
598                                    const Char_t *sigmas, Float_t sigmaStep) const {
599    //
600    // Creates a histogram, where you can see, how much of the data are inside sigma-intervals 
601    // around the mean/median/LTM
602    // with drawCommand, sector and cuts you specify your input data, see EasyDraw
603    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma)
604    // sigmaStep: the binsize of the generated histogram
605    // plotMean/plotMedian/plotLTM: specifies where to put the center
606    //
607   
608    Double_t ltmFraction = 0.8;
609    
610    TString drawStr(drawCommand);
611    Bool_t dangerousToDraw = drawStr.Contains(":") || drawStr.Contains(">>");
612    if (dangerousToDraw) {
613       Warning("SigmaCut", "The draw string must not contain ':' or '>>'.");
614       return -1;
615    }
616    drawStr += " >> tempHist";
617    
618    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
619    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
620    // FIXME is this histogram deleted automatically?
621    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
622    
623    Double_t mean = TMath::Mean(entries, values);
624    Double_t median = TMath::Median(entries, values);
625    Double_t sigma = TMath::RMS(entries, values);
626    
627    TLegend * legend = new TLegend(.7,.7, .99, .99, "Cumulative");
628    //fListOfObjectsToBeDeleted->Add(legend);
629    TH1F *cutHistoMean = 0;
630    TH1F *cutHistoMedian = 0;
631    TH1F *cutHistoLTM = 0;
632    
633    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");  
634    TVectorF nsigma(sigmasTokens->GetEntriesFast());
635    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
636       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
637       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
638       nsigma[i] = sig;
639    }
640   
641    if (plotMean) {
642       cutHistoMean = SigmaCut(htemp, mean, sigma, sigmaMax, sigmaStep, pm);
643       if (cutHistoMean) {
644          //fListOfObjectsToBeDeleted->Add(cutHistoMean);
645          cutHistoMean->SetLineColor(kRed);
646          legend->AddEntry(cutHistoMean, "Mean", "l");
647          cutHistoMean->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
648          cutHistoMean->Draw();
649          DrawLines(cutHistoMean, nsigma, legend, kRed, pm);
650       } // if (cutHistoMean)
651        
652    }
653    if (plotMedian) {
654       cutHistoMedian = SigmaCut(htemp, median, sigma, sigmaMax, sigmaStep, pm);
655       if (cutHistoMedian) {
656          //fListOfObjectsToBeDeleted->Add(cutHistoMedian);
657          cutHistoMedian->SetLineColor(kBlue);
658          legend->AddEntry(cutHistoMedian, "Median", "l");
659          cutHistoMedian->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
660          if (plotMean && cutHistoMean) cutHistoMedian->Draw("same");
661             else cutHistoMedian->Draw();
662          DrawLines(cutHistoMedian, nsigma, legend, kBlue, pm);
663       }  // if (cutHistoMedian)
664    }
665    if (plotLTM) {
666       Double_t ltmRms = 0;
667       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
668       cutHistoLTM = SigmaCut(htemp, ltm, ltmRms, sigmaMax, sigmaStep, pm);
669       if (cutHistoLTM) {
670          //fListOfObjectsToBeDeleted->Add(cutHistoLTM);
671          cutHistoLTM->SetLineColor(kGreen+2);
672          legend->AddEntry(cutHistoLTM, "LTM", "l");
673          cutHistoLTM->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
674          if ((plotMean && cutHistoMean) || (plotMedian && cutHistoMedian)) cutHistoLTM->Draw("same");
675             else cutHistoLTM->Draw();
676          DrawLines(cutHistoLTM, nsigma, legend, kGreen+2, pm);
677       }
678    }
679    if (!plotMean && !plotMedian && !plotLTM) return -1;
680    legend->Draw();
681    return 1;
682 }
683
684 //_____________________________________________________________________________
685 Int_t AliBaseCalibViewer::Integrate(const Char_t* drawCommand, const Char_t* sector, const Char_t* cuts, 
686                                     Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, 
687                                     const Char_t *sigmas, Float_t sigmaStep) const {
688    //
689    // 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"   
690    // "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
691    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate 
692    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
693    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
694    // The actual work is done on the array.
695    /* Begin_Latex 
696          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 }
697       End_Latex  
698    */
699    
700    Double_t ltmFraction = 0.8;
701    // avoid compiler warnings:
702    sigmaMax = sigmaMax;
703    sigmaStep = sigmaStep;
704    
705    TString drawStr(drawCommand);
706    Bool_t dangerousToDraw = drawStr.Contains(":") || drawStr.Contains(">>");
707    if (dangerousToDraw) {
708       Warning("Integrate", "The draw string must not contain ':' or '>>'.");
709       return -1;
710    }
711    drawStr += " >> tempHist";
712    
713    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
714    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
715    TGraph *integralGraphMean   = 0;
716    TGraph *integralGraphMedian = 0;
717    TGraph *integralGraphLTM    = 0;
718    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
719    Int_t    *index  = new Int_t[entries];
720    Float_t  *xarray = new Float_t[entries];
721    Float_t  *yarray = new Float_t[entries];
722    TMath::Sort(entries, values, index, kFALSE);
723    
724    Double_t mean = TMath::Mean(entries, values);
725    Double_t median = TMath::Median(entries, values);
726    Double_t sigma = TMath::RMS(entries, values);
727    
728    // parse sigmas string
729    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");  
730    TVectorF nsigma(sigmasTokens->GetEntriesFast());
731    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
732       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
733       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
734       nsigma[i] = sig;
735    }
736   
737    TLegend * legend = new TLegend(.7,.7, .99, .99, "Integrated histogram");
738    //fListOfObjectsToBeDeleted->Add(legend);
739   
740    if (plotMean) {
741       for (Int_t i = 0; i < entries; i++) {
742          xarray[i] = (values[index[i]] - mean) / sigma; 
743          yarray[i] = float(i) / float(entries);
744       }
745       integralGraphMean = new TGraph(entries, xarray, yarray);
746       if (integralGraphMean) {
747          //fListOfObjectsToBeDeleted->Add(integralGraphMean);
748          integralGraphMean->SetLineColor(kRed);
749          legend->AddEntry(integralGraphMean, "Mean", "l");
750          integralGraphMean->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
751          integralGraphMean->Draw("alu");
752          DrawLines(integralGraphMean, nsigma, legend, kRed, kTRUE);
753       }
754    }
755    if (plotMedian) {
756       for (Int_t i = 0; i < entries; i++) {
757          xarray[i] = (values[index[i]] - median) / sigma; 
758          yarray[i] = float(i) / float(entries);
759       }
760       integralGraphMedian = new TGraph(entries, xarray, yarray);
761       if (integralGraphMedian) {
762          //fListOfObjectsToBeDeleted->Add(integralGraphMedian);
763          integralGraphMedian->SetLineColor(kBlue);
764          legend->AddEntry(integralGraphMedian, "Median", "l");
765          integralGraphMedian->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
766          if (plotMean && integralGraphMean) integralGraphMedian->Draw("samelu");
767             else integralGraphMedian->Draw("alu");
768          DrawLines(integralGraphMedian, nsigma, legend, kBlue, kTRUE);
769       }
770    }
771    if (plotLTM) {
772       Double_t ltmRms = 0;
773       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
774       for (Int_t i = 0; i < entries; i++) {
775          xarray[i] = (values[index[i]] - ltm) / ltmRms; 
776          yarray[i] = float(i) / float(entries);
777       }
778       integralGraphLTM = new TGraph(entries, xarray, yarray);
779       if (integralGraphLTM) {
780          //fListOfObjectsToBeDeleted->Add(integralGraphLTM);
781          integralGraphLTM->SetLineColor(kGreen+2);
782          legend->AddEntry(integralGraphLTM, "LTM", "l");
783          integralGraphLTM->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
784          if ((plotMean && integralGraphMean) || (plotMedian && integralGraphMedian)) integralGraphLTM->Draw("samelu");
785             else integralGraphLTM->Draw("alu");
786          DrawLines(integralGraphLTM, nsigma, legend, kGreen+2, kTRUE);
787       }
788    }
789    delete [] index;
790    delete [] xarray;
791    delete [] yarray;
792    if (!plotMean && !plotMedian && !plotLTM) return -1;
793    legend->Draw();
794    return entries;
795 }
796
797 //_____________________________________________________________________________
798 TH1F* AliBaseCalibViewer::Integrate(TH1F *histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
799    //
800    // 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"   
801    // "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
802    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate 
803    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
804    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
805    // The actual work is done on the array.
806    /* Begin_Latex 
807          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 }
808       End_Latex  
809       begin_macro(source)
810       {
811          Float_t mean = 0;
812          Float_t sigma = 1.5;
813          Float_t sigmaMax = 4;
814          gROOT->SetStyle("Plain");
815          TH1F *distribution = new TH1F("Distribution2", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
816          TRandom rand(23);
817          for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
818          Float_t *ar = distribution->GetArray();
819          
820          TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas_Integrate", "", 350, 350);
821          macro_example_canvas->Divide(0,2);
822          TVirtualPad *pad1 = macro_example_canvas->cd(1);
823          pad1->SetGridy();
824          pad1->SetGridx();
825          distribution->Draw();
826          TVirtualPad *pad2 = macro_example_canvas->cd(2);
827          pad2->SetGridy();
828          pad2->SetGridx();
829          TH1F *shist = AliTPCCalibViewer::Integrate(distribution, mean, sigma, sigmaMax);
830          shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
831          shist->Draw();  
832          
833          return macro_example_canvas_Integrate;
834       }  
835       end_macro
836    */ 
837
838    
839    Float_t *array = histogram->GetArray();
840    Int_t    nbins = histogram->GetXaxis()->GetNbins();
841    Float_t binLow = histogram->GetXaxis()->GetXmin();
842    Float_t binUp  = histogram->GetXaxis()->GetXmax();
843    return Integrate(nbins, array, nbins, binLow, binUp, mean, sigma, sigmaMax, sigmaStep);
844 }
845
846 //_____________________________________________________________________________
847 TH1F* AliBaseCalibViewer::Integrate(Int_t n, Float_t *array, Int_t nbins, Float_t binLow, Float_t binUp, 
848                                     Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
849    // 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"   
850    // "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
851    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate 
852    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
853    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
854    // Here the actual work is done.
855       
856    Bool_t givenUnits = kTRUE;
857    if (TMath::Abs(sigma) < 1.e-10 && TMath::Abs(sigmaMax) < 1.e-10) givenUnits = kFALSE;
858    if (givenUnits) {
859       sigma = 1;
860       sigmaMax = (binUp - binLow) / 2.;
861    }
862    
863    Float_t binWidth = (binUp-binLow)/(nbins - 1);
864    if (sigmaStep <= 0) sigmaStep = binWidth;
865    Int_t kbins =  (Int_t)(sigmaMax * sigma / sigmaStep) + 1;  // + 1  due to overflow bin in histograms
866    Float_t kbinLow = givenUnits ? binLow : -sigmaMax;
867    Float_t kbinUp  = givenUnits ? binUp  : sigmaMax;
868    TH1F *hist = 0; 
869    if (givenUnits)  hist = new TH1F("integratedHisto","Integrated Histogram; Given x; Fraction of included data", kbins, kbinLow, kbinUp); 
870    if (!givenUnits) hist = new TH1F("integratedHisto","Integrated Histogram; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp); 
871    hist->SetDirectory(0);
872    hist->Reset();
873    
874    // calculate normalization
875  //  printf("calculating normalization, integrating from bin 1 to %i \n", n);
876    Double_t normalization = 0;
877    for (Int_t i = 1; i <= n; i++) {
878         normalization += array[i];
879    }
880  //  printf("normalization: %f \n", normalization);
881    
882    // given units: units from given histogram
883    // sigma units: in units of sigma
884    // iDelta: integrate in interval (mean +- iDelta), given units
885    // x:      ofset from mean for integration, given units
886    // hist:   needs 
887    
888    // fill histogram
889    for (Float_t iDelta = mean - sigmaMax * sigma; iDelta <= mean + sigmaMax * sigma; iDelta += sigmaStep) {
890       // integrate array
891       Double_t value = 0;
892       for (Float_t x = mean - sigmaMax * sigma; x <= iDelta; x += binWidth) {
893          value += (x <= binUp && x >= binLow)  ? array[GetBin(x, nbins, binLow, binUp)] : 0;
894       }
895       if (value / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", value, normalization);
896       if (value / normalization > 100) return hist;
897       Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
898     //  printf("first integration bin: %i, last integration bin: %i \n", GetBin(mean - sigmaMax * sigma, nbins, binLow, binUp), GetBin(iDelta, nbins, binLow, binUp));
899     //  printf("value: %f, normalization: %f, normalized value: %f, iDelta: %f, Bin: %i \n", value, normalization, value/normalization, iDelta, bin);
900       value = (value / normalization);
901       hist->SetBinContent(bin, value);
902    }
903    return hist;
904 }
905
906 //_____________________________________________________________________________
907 void AliBaseCalibViewer::DrawLines(TH1F *histogram, TVectorF nsigma, TLegend *legend, Int_t color, Bool_t pm) const {
908    // 
909    // Private function for SigmaCut(...) and Integrate(...)
910    // Draws lines into the given histogram, specified by "nsigma", the lines are addeed to the legend
911    // 
912    
913    // start to draw the lines, loop over requested sigmas
914    Char_t c[500];
915    for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
916       if (!pm) { 
917          Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
918          TLine* lineUp = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
919          //fListOfObjectsToBeDeleted->Add(lineUp);
920          lineUp->SetLineColor(color);
921          lineUp->SetLineStyle(2 + i);
922          lineUp->Draw();
923          TLine* lineLeft = new TLine(nsigma[i], histogram->GetBinContent(bin), 0, histogram->GetBinContent(bin));
924          //fListOfObjectsToBeDeleted->Add(lineLeft);
925          lineLeft->SetLineColor(color);
926          lineLeft->SetLineStyle(2 + i);
927          lineLeft->Draw();
928          sprintf(c, "Fraction(%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
929          legend->AddEntry(lineLeft, c, "l");
930       }
931       else { // if (pm)
932          Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
933          TLine* lineUp1 = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
934          //fListOfObjectsToBeDeleted->Add(lineUp1);
935          lineUp1->SetLineColor(color);
936          lineUp1->SetLineStyle(2 + i);
937          lineUp1->Draw();
938          TLine* lineLeft1 = new TLine(nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
939          //fListOfObjectsToBeDeleted->Add(lineLeft1);
940          lineLeft1->SetLineColor(color);
941          lineLeft1->SetLineStyle(2 + i);
942          lineLeft1->Draw();
943          sprintf(c, "Fraction(+%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
944          legend->AddEntry(lineLeft1, c, "l");
945          bin = histogram->GetXaxis()->FindBin(-nsigma[i]);
946          TLine* lineUp2 = new TLine(-nsigma[i], 0, -nsigma[i], histogram->GetBinContent(bin));
947          //fListOfObjectsToBeDeleted->Add(lineUp2);
948          lineUp2->SetLineColor(color);
949          lineUp2->SetLineStyle(2 + i);
950          lineUp2->Draw();
951          TLine* lineLeft2 = new TLine(-nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
952          //fListOfObjectsToBeDeleted->Add(lineLeft2);
953          lineLeft2->SetLineColor(color);
954          lineLeft2->SetLineStyle(2 + i);
955          lineLeft2->Draw();
956          sprintf(c, "Fraction(-%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
957          legend->AddEntry(lineLeft2, c, "l");
958       }
959    }  // for (Int_t i = 0; i < nsigma.GetNoElements(); i++)   
960 }
961
962 //_____________________________________________________________________________
963 void AliBaseCalibViewer::DrawLines(TGraph *graph, TVectorF nsigma, TLegend *legend, Int_t color, Bool_t pm) const {
964    // 
965    // Private function for SigmaCut(...) and Integrate(...)
966    // Draws lines into the given histogram, specified by "nsigma", the lines are addeed to the legend
967    // 
968    
969    // start to draw the lines, loop over requested sigmas
970    Char_t c[500];
971    for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
972       if (!pm) { 
973          TLine* lineUp = new TLine(nsigma[i], 0, nsigma[i], graph->Eval(nsigma[i]));
974          //fListOfObjectsToBeDeleted->Add(lineUp);
975          lineUp->SetLineColor(color);
976          lineUp->SetLineStyle(2 + i);
977          lineUp->Draw();
978          TLine* lineLeft = new TLine(nsigma[i], graph->Eval(nsigma[i]), 0, graph->Eval(nsigma[i]));
979          //fListOfObjectsToBeDeleted->Add(lineLeft);
980          lineLeft->SetLineColor(color);
981          lineLeft->SetLineStyle(2 + i);
982          lineLeft->Draw();
983          sprintf(c, "Fraction(%f #sigma) = %f",nsigma[i], graph->Eval(nsigma[i]));
984          legend->AddEntry(lineLeft, c, "l");
985       }
986       else { // if (pm)
987          TLine* lineUp1 = new TLine(nsigma[i], 0, nsigma[i], graph->Eval(nsigma[i]));
988          //fListOfObjectsToBeDeleted->Add(lineUp1);
989          lineUp1->SetLineColor(color);
990          lineUp1->SetLineStyle(2 + i);
991          lineUp1->Draw();
992          TLine* lineLeft1 = new TLine(nsigma[i], graph->Eval(nsigma[i]), graph->GetHistogram()->GetXaxis()->GetBinLowEdge(0), graph->Eval(nsigma[i]));
993          //fListOfObjectsToBeDeleted->Add(lineLeft1);
994          lineLeft1->SetLineColor(color);
995          lineLeft1->SetLineStyle(2 + i);
996          lineLeft1->Draw();
997          sprintf(c, "Fraction(+%f #sigma) = %f",nsigma[i], graph->Eval(nsigma[i]));
998          legend->AddEntry(lineLeft1, c, "l");
999          TLine* lineUp2 = new TLine(-nsigma[i], 0, -nsigma[i], graph->Eval(-nsigma[i]));
1000          //fListOfObjectsToBeDeleted->Add(lineUp2);
1001          lineUp2->SetLineColor(color);
1002          lineUp2->SetLineStyle(2 + i);
1003          lineUp2->Draw();
1004          TLine* lineLeft2 = new TLine(-nsigma[i], graph->Eval(-nsigma[i]), graph->GetHistogram()->GetXaxis()->GetBinLowEdge(0), graph->Eval(-nsigma[i]));
1005          //fListOfObjectsToBeDeleted->Add(lineLeft2);
1006          lineLeft2->SetLineColor(color);
1007          lineLeft2->SetLineStyle(2 + i);
1008          lineLeft2->Draw();
1009          sprintf(c, "Fraction(-%f #sigma) = %f",nsigma[i], graph->Eval(-nsigma[i]));
1010          legend->AddEntry(lineLeft2, c, "l");
1011       }
1012    }  // for (Int_t i = 0; i < nsigma.GetNoElements(); i++)   
1013 }