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    delete [] index;
 782    delete [] xarray;
 783    delete [] yarray;
 784    if (!plotMean && !plotMedian && !plotLTM) return -1;
 785    legend->Draw();
 786    return entries;
 787 }
 788
 789 //_____________________________________________________________________________
 790 TH1F* AliBaseCalibViewer::Integrate(TH1F *histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
 791    //
 792    // 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"
 793    // "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
 794    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
 795    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
 796    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
 797    // The actual work is done on the array.
 798    /* Begin_Latex
 799          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 }
 800       End_Latex
 801       begin_macro(source)
 802       {
 803          Float_t mean = 0;
 804          Float_t sigma = 1.5;
 805          Float_t sigmaMax = 4;
 806          gROOT->SetStyle("Plain");
 807          TH1F *distribution = new TH1F("Distribution2", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
 808          TRandom rand(23);
 809          for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
 810          Float_t *ar = distribution->GetArray();
 811
 812          TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas_Integrate", "", 350, 350);
 813          macro_example_canvas->Divide(0,2);
 814          TVirtualPad *pad1 = macro_example_canvas->cd(1);
 815          pad1->SetGridy();
 816          pad1->SetGridx();
 817          distribution->Draw();
 818          TVirtualPad *pad2 = macro_example_canvas->cd(2);
 819          pad2->SetGridy();
 820          pad2->SetGridx();
 821          TH1F *shist = AliTPCCalibViewer::Integrate(distribution, mean, sigma, sigmaMax);
 822          shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
 823          shist->Draw();
 824
 825          return macro_example_canvas_Integrate;
 826       }
 827       end_macro
 828    */
 829
 830
 831    Float_t *array = histogram->GetArray();
 832    Int_t    nbins = histogram->GetXaxis()->GetNbins();
 833    Float_t binLow = histogram->GetXaxis()->GetXmin();
 834    Float_t binUp  = histogram->GetXaxis()->GetXmax();
 835    return Integrate(nbins, array, nbins, binLow, binUp, mean, sigma, sigmaMax, sigmaStep);
 836 }
 837
 838 //_____________________________________________________________________________
 839 TH1F* AliBaseCalibViewer::Integrate(Int_t n, Float_t *array, Int_t nbins, Float_t binLow, Float_t binUp,
 840                                     Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
 841    // 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"
 842    // "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
 843    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
 844    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
 845    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
 846    // Here the actual work is done.
 847
 848    Bool_t givenUnits = kTRUE;
 849    if (TMath::Abs(sigma) < 1.e-10 && TMath::Abs(sigmaMax) < 1.e-10) givenUnits = kFALSE;
 850    if (givenUnits) {
 851       sigma = 1;
 852       sigmaMax = (binUp - binLow) / 2.;
 853    }
 854
 855    Float_t binWidth = (binUp-binLow)/(nbins - 1);
 856    if (sigmaStep <= 0) sigmaStep = binWidth;
 857    Int_t kbins =  (Int_t)(sigmaMax * sigma / sigmaStep) + 1;  // + 1  due to overflow bin in histograms
 858    Float_t kbinLow = givenUnits ? binLow : -sigmaMax;
 859    Float_t kbinUp  = givenUnits ? binUp  : sigmaMax;
 860    TH1F *hist = 0;
 861    if (givenUnits)  hist = new TH1F("integratedHisto","Integrated Histogram; Given x; Fraction of included data", kbins, kbinLow, kbinUp);
 862    if (!givenUnits) hist = new TH1F("integratedHisto","Integrated Histogram; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp);
 863    hist->SetDirectory(0);
 864    hist->Reset();
 865
 866    // calculate normalization
 867  //  printf("calculating normalization, integrating from bin 1 to %i \n", n);
 868    Double_t normalization = 0;
 869    for (Int_t i = 1; i <= n; i++) {
 870         normalization += array[i];
 871    }
 872  //  printf("normalization: %f \n", normalization);
 873
 874    // given units: units from given histogram
 875    // sigma units: in units of sigma
 876    // iDelta: integrate in interval (mean +- iDelta), given units
 877    // x:      ofset from mean for integration, given units
 878    // hist:   needs
 879
 880    // fill histogram
 881    for (Float_t iDelta = mean - sigmaMax * sigma; iDelta <= mean + sigmaMax * sigma; iDelta += sigmaStep) {
 882       // integrate array
 883       Double_t value = 0;
 884       for (Float_t x = mean - sigmaMax * sigma; x <= iDelta; x += binWidth) {
 885          value += (x <= binUp && x >= binLow)  ? array[GetBin(x, nbins, binLow, binUp)] : 0;
 886       }
 887       if (value / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", value, normalization);
 888       if (value / normalization > 100) return hist;
 889       Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
 890     //  printf("first integration bin: %i, last integration bin: %i \n", GetBin(mean - sigmaMax * sigma, nbins, binLow, binUp), GetBin(iDelta, nbins, binLow, binUp));
 891     //  printf("value: %f, normalization: %f, normalized value: %f, iDelta: %f, Bin: %i \n", value, normalization, value/normalization, iDelta, bin);
 892       value = (value / normalization);
 893       hist->SetBinContent(bin, value);
 894    }
 895    return hist;
 896 }
 897
 898 //_____________________________________________________________________________
 899 void AliBaseCalibViewer::DrawLines(TH1F *histogram, TVectorF nsigma, TLegend *legend, Int_t color, Bool_t pm) const {
 900    //
 901    // Private function for SigmaCut(...) and Integrate(...)
 902    // Draws lines into the given histogram, specified by "nsigma", the lines are addeed to the legend
 903    //
 904
 905    // start to draw the lines, loop over requested sigmas
 906    Char_t c[500];
 907    for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
 908       if (!pm) {
 909          Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
 910          TLine* lineUp = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
 911          //fListOfObjectsToBeDeleted->Add(lineUp);
 912          lineUp->SetLineColor(color);
 913          lineUp->SetLineStyle(2 + i);
 914          lineUp->Draw();
 915          TLine* lineLeft = new TLine(nsigma[i], histogram->GetBinContent(bin), 0, histogram->GetBinContent(bin));
 916          //fListOfObjectsToBeDeleted->Add(lineLeft);
 917          lineLeft->SetLineColor(color);
 918          lineLeft->SetLineStyle(2 + i);
 919          lineLeft->Draw();
 920          sprintf(c, "Fraction(%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
 921          legend->AddEntry(lineLeft, c, "l");
 922       }
 923       else { // if (pm)
 924          Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
 925          TLine* lineUp1 = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
 926          //fListOfObjectsToBeDeleted->Add(lineUp1);
 927          lineUp1->SetLineColor(color);
 928          lineUp1->SetLineStyle(2 + i);
 929          lineUp1->Draw();
 930          TLine* lineLeft1 = new TLine(nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
 931          //fListOfObjectsToBeDeleted->Add(lineLeft1);
 932          lineLeft1->SetLineColor(color);
 933          lineLeft1->SetLineStyle(2 + i);
 934          lineLeft1->Draw();
 935          sprintf(c, "Fraction(+%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
 936          legend->AddEntry(lineLeft1, c, "l");
 937          bin = histogram->GetXaxis()->FindBin(-nsigma[i]);
 938          TLine* lineUp2 = new TLine(-nsigma[i], 0, -nsigma[i], histogram->GetBinContent(bin));
 939          //fListOfObjectsToBeDeleted->Add(lineUp2);
 940          lineUp2->SetLineColor(color);
 941          lineUp2->SetLineStyle(2 + i);
 942          lineUp2->Draw();
 943          TLine* lineLeft2 = new TLine(-nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
 944          //fListOfObjectsToBeDeleted->Add(lineLeft2);
 945          lineLeft2->SetLineColor(color);
 946          lineLeft2->SetLineStyle(2 + i);
 947          lineLeft2->Draw();
 948          sprintf(c, "Fraction(-%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
 949          legend->AddEntry(lineLeft2, c, "l");
 950       }
 951    }  // for (Int_t i = 0; i < nsigma.GetNoElements(); i++)
 952 }
 953
 954 //_____________________________________________________________________________
 955 void AliBaseCalibViewer::DrawLines(TGraph *graph, TVectorF nsigma, TLegend *legend, Int_t color, Bool_t pm) const {
 956    //
 957    // Private function for SigmaCut(...) and Integrate(...)
 958    // Draws lines into the given histogram, specified by "nsigma", the lines are addeed to the legend
 959    //
 960
 961    // start to draw the lines, loop over requested sigmas
 962    Char_t c[500];
 963    for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
 964       if (!pm) {
 965          TLine* lineUp = new TLine(nsigma[i], 0, nsigma[i], graph->Eval(nsigma[i]));
 966          //fListOfObjectsToBeDeleted->Add(lineUp);
 967          lineUp->SetLineColor(color);
 968          lineUp->SetLineStyle(2 + i);
 969          lineUp->Draw();
 970          TLine* lineLeft = new TLine(nsigma[i], graph->Eval(nsigma[i]), 0, graph->Eval(nsigma[i]));
 971          //fListOfObjectsToBeDeleted->Add(lineLeft);
 972          lineLeft->SetLineColor(color);
 973          lineLeft->SetLineStyle(2 + i);
 974          lineLeft->Draw();
 975          sprintf(c, "Fraction(%f #sigma) = %f",nsigma[i], graph->Eval(nsigma[i]));
 976          legend->AddEntry(lineLeft, c, "l");
 977       }
 978       else { // if (pm)
 979          TLine* lineUp1 = new TLine(nsigma[i], 0, nsigma[i], graph->Eval(nsigma[i]));
 980          //fListOfObjectsToBeDeleted->Add(lineUp1);
 981          lineUp1->SetLineColor(color);
 982          lineUp1->SetLineStyle(2 + i);
 983          lineUp1->Draw();
 984          TLine* lineLeft1 = new TLine(nsigma[i], graph->Eval(nsigma[i]), graph->GetHistogram()->GetXaxis()->GetBinLowEdge(0), graph->Eval(nsigma[i]));
 985          //fListOfObjectsToBeDeleted->Add(lineLeft1);
 986          lineLeft1->SetLineColor(color);
 987          lineLeft1->SetLineStyle(2 + i);
 988          lineLeft1->Draw();
 989          sprintf(c, "Fraction(+%f #sigma) = %f",nsigma[i], graph->Eval(nsigma[i]));
 990          legend->AddEntry(lineLeft1, c, "l");
 991          TLine* lineUp2 = new TLine(-nsigma[i], 0, -nsigma[i], graph->Eval(-nsigma[i]));
 992          //fListOfObjectsToBeDeleted->Add(lineUp2);
 993          lineUp2->SetLineColor(color);
 994          lineUp2->SetLineStyle(2 + i);
 995          lineUp2->Draw();
 996          TLine* lineLeft2 = new TLine(-nsigma[i], graph->Eval(-nsigma[i]), graph->GetHistogram()->GetXaxis()->GetBinLowEdge(0), graph->Eval(-nsigma[i]));
 997          //fListOfObjectsToBeDeleted->Add(lineLeft2);
 998          lineLeft2->SetLineColor(color);
 999          lineLeft2->SetLineStyle(2 + i);
1000          lineLeft2->Draw();
1001          sprintf(c, "Fraction(-%f #sigma) = %f",nsigma[i], graph->Eval(-nsigma[i]));
1002          legend->AddEntry(lineLeft2, c, "l");
1003       }
1004    }  // for (Int_t i = 0; i < nsigma.GetNoElements(); i++)
1005 }