TPC/AliTPCCalibViewer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //  Class for viewing/visualizing TPC calibration data                       //
  20 //  base on  TTree functionality for visualization                           //
  21 //                                                                           //
  22 //  Create a list of AliTPCCalPads, arrange them in an TObjArray.            //
  23 //  Pass this TObjArray to MakeTree and create the calibration Tree          //
  24 //  While craating this tree some statistical information are calculated     //
  25 //  Open the viewer with this Tree: AliTPCCalibViewer v("CalibTree.root")    //
  26 //  Have fun!                                                                //
  27 //  EasyDraw("CETmean~-CETmean_mean", "A", "(CETmean~-CETmean_mean)>0")      //
  28 //                                                                           //
  29 //  If you like to click, we recommand you the                               //
  30 //    AliTPCCalibViewerGUI                                                   //
  31 //                                                                           //
  32 //    THE DOCUMENTATION IS STILL NOT COMPLETED !!!!                          //
  33 //                                                                           //
  34 ///////////////////////////////////////////////////////////////////////////////
  35
  36 //
  37 // ROOT includes
  38 //
  39 #include <iostream>
  40 #include <fstream>
  41 #include <TString.h>
  42 #include <TRandom.h>
  43 #include <TLegend.h>
  44 #include <TLine.h>
  45 #include <TCanvas.h>
  46 #include <TROOT.h>
  47 #include <TStyle.h>
  48 #include <TH1.h>
  49 #include <TH1F.h>
  50 #include <THashTable.h>
  51 #include <TObjString.h>
  52 #include "TTreeStream.h"
  53 #include "TFile.h"
  54 #include "TKey.h"
  55 #include "TGraph.h"
  56 #include "AliTPCCalibPulser.h"
  57 #include "AliTPCCalibPedestal.h"
  58 #include "AliTPCCalibCE.h"
  59 // #include "TObjArray.h"
  60 // #include "TObjString.h"
  61 // #include "TString.h"
  62 // #include "AliTPCCalPad.h"
  63
  64
  65 //
  66 // AliRoot includes
  67 //
  68 #include "AliTPCCalibViewer.h"
  69
  70 ClassImp(AliTPCCalibViewer)
  71
  72 AliTPCCalibViewer::AliTPCCalibViewer()
  73                   :TObject(),
  74                    fTree(0),
  75                    fFile(0),
  76                    fListOfObjectsToBeDeleted(0),
  77                    fTreeMustBeDeleted(0)
  78 {
  79   //
  80   // Default constructor
  81   //
  82
  83 }
  84
  85 //_____________________________________________________________________________
  86 AliTPCCalibViewer::AliTPCCalibViewer(const AliTPCCalibViewer &c)
  87                   :TObject(c),
  88                    fTree(0),
  89                    fFile(0),
  90                    fListOfObjectsToBeDeleted(0),
  91                    fTreeMustBeDeleted(0)
  92 {
  93   //
  94   // dummy AliTPCCalibViewer copy constructor
  95   // not yet working!!!
  96   //
  97   fTree = c.fTree;
  98   fTreeMustBeDeleted = c.fTreeMustBeDeleted;
  99   //fFile = new TFile(*(c.fFile));
 100   fListOfObjectsToBeDeleted = c.fListOfObjectsToBeDeleted;
 101 }
 102
 103 //_____________________________________________________________________________
 104 AliTPCCalibViewer::AliTPCCalibViewer(TTree* tree)
 105                   :TObject(),
 106                    fTree(0),
 107                    fFile(0),
 108                    fListOfObjectsToBeDeleted(0),
 109                    fTreeMustBeDeleted(0)
 110 {
 111   //
 112   // Constructor that initializes the calibration viewer
 113   //
 114   fTree = tree;
 115   fTreeMustBeDeleted = kFALSE;
 116   fListOfObjectsToBeDeleted = new TObjArray();
 117 }
 118
 119 //_____________________________________________________________________________
 120 AliTPCCalibViewer::AliTPCCalibViewer(char* fileName, char* treeName)
 121                   :TObject(),
 122                    fTree(0),
 123                    fFile(0),
 124                    fListOfObjectsToBeDeleted(0),
 125                    fTreeMustBeDeleted(0)
 126 {
 127    //
 128    // Constructor to initialize the calibration viewer
 129    // the file 'fileName' contains the tree 'treeName'
 130    //
 131    fFile = new TFile(fileName, "read");
 132    fTree = (TTree*) fFile->Get(treeName);
 133    fTreeMustBeDeleted = kTRUE;
 134    fListOfObjectsToBeDeleted = new TObjArray();
 135 }
 136
 137 //____________________________________________________________________________
 138 AliTPCCalibViewer & AliTPCCalibViewer::operator =(const AliTPCCalibViewer & param)
 139 {
 140    //
 141    // assignment operator - dummy
 142    // not yet working!!!
 143    //
 144    fTree = param.fTree;
 145    fTreeMustBeDeleted = param.fTreeMustBeDeleted;
 146    //fFile = new TFile(*(param.fFile));
 147    fListOfObjectsToBeDeleted = param.fListOfObjectsToBeDeleted;
 148    return (*this);
 149 }
 150
 151 //_____________________________________________________________________________
 152 AliTPCCalibViewer::~AliTPCCalibViewer()
 153 {
 154    //
 155    // AliTPCCalibViewer destructor
 156    // all objects will be deleted, the file will be closed, the pictures will disappear
 157    //
 158    if (fTree && fTreeMustBeDeleted) {
 159       fTree->SetCacheSize(0);
 160       fTree->Delete();
 161       //delete fTree;
 162    }
 163    if (fFile) {
 164       fFile->Close();
 165       fFile = 0;
 166    }
 167
 168    for (Int_t i = fListOfObjectsToBeDeleted->GetEntriesFast()-1; i >= 0; i--) {
 169       //cout << "Index " << i << " trying to delete the following object: " << fListOfObjectsToBeDeleted->At(i)->GetName() << "..."<< endl;
 170       delete fListOfObjectsToBeDeleted->At(i);
 171    }
 172    delete fListOfObjectsToBeDeleted;
 173 }
 174
 175 //_____________________________________________________________________________
 176 void AliTPCCalibViewer::Delete(Option_t* option) {
 177    //
 178    // Should be called from AliTPCCalibViewerGUI class only.
 179    // If you use Delete() do not call the destructor.
 180    // All objects (except those contained in fListOfObjectsToBeDeleted) will be deleted, the file will be closed.
 181    //
 182
 183    option = option;  // to avoid warnings on compiling
 184    if (fTree && fTreeMustBeDeleted) {
 185       fTree->SetCacheSize(0);
 186       fTree->Delete();
 187    }
 188    if (fFile)
 189       delete fFile;
 190    delete fListOfObjectsToBeDeleted;
 191 }
 192
 193 //_____________________________________________________________________________
 194 Int_t AliTPCCalibViewer::EasyDraw(const char* drawCommand, const char* sector, const char* cuts, const char* drawOptions, Bool_t writeDrawCommand) const {
 195   //
 196   // easy drawing of data, use '~' for abbreviation of '.fElements'
 197   // example: EasyDraw("CETmean~-CETmean_mean", "A", "(CETmean~-CETmean_mean)>0")
 198  // sector: sector-number - only the specified sector will be drwawn
 199   //         'A'/'C' or 'a'/'c' - side A/C will be drawn
 200   //         'ALL' - whole TPC will be drawn, projected on one side
 201   // cuts: specifies cuts
 202   // drawOptions: draw options like 'same'
 203   // writeDrawCommand: write the command, that is passed to TTree::Draw
 204   //
 205
 206    TString drawStr(drawCommand);
 207    TString sectorStr(sector);
 208    sectorStr.ToUpper();
 209    TString cutStr("");
 210    //TString drawOptionsStr("profcolz ");
 211    TString drawOptionsStr("");
 212    TRandom rnd(0);
 213    Int_t rndNumber = rnd.Integer(10000);
 214
 215    if (drawOptions && strcmp(drawOptions, "") != 0)
 216       drawOptionsStr += drawOptions;
 217    else
 218       drawOptionsStr += "profcolz";
 219
 220    if (sectorStr == "A") {
 221       drawStr += ":gy.fElements:gx.fElements>>prof";
 222       drawStr += rndNumber;
 223       drawStr += "(330,-250,250,330,-250,250)";
 224       cutStr += "(sector/18)%2==0 ";
 225    }
 226    else if  (sectorStr == "C") {
 227       drawStr += ":gy.fElements:gx.fElements>>prof";
 228       drawStr += rndNumber;
 229       drawStr += "(330,-250,250,330,-250,250)";
 230       cutStr += "(sector/18)%2==1 ";
 231    }
 232    else if  (sectorStr == "ALL") {
 233       drawStr += ":gy.fElements:gx.fElements>>prof";
 234       drawStr += rndNumber;
 235       drawStr += "(330,-250,250,330,-250,250)";
 236    }
 237    else if (sectorStr.IsDigit()) {
 238       Int_t isec = sectorStr.Atoi();
 239       drawStr += ":rpad.fElements:row.fElements>>prof";
 240       drawStr += rndNumber;
 241       if (isec < 36 && isec >= 0)
 242          drawStr += "(63,0,63,108,-54,54)";
 243       else if (isec < 72 && isec >= 36)
 244          drawStr += "(96,0,96,140,-70,70)";
 245       else {
 246          Error("EasyDraw","The TPC contains only sectors between 0 and 71.");
 247          return -1;
 248       }
 249       cutStr += "(sector==";
 250       cutStr += isec;
 251       cutStr += ") ";
 252    }
 253
 254    if (cuts && cuts[0] != 0) {
 255       if (cutStr.Length() != 0) cutStr += "&& ";
 256       cutStr += "(";
 257       cutStr += cuts;
 258       cutStr += ")";
 259    }
 260    drawStr.ReplaceAll("~", ".fElements");
 261    cutStr.ReplaceAll("~", ".fElements");
 262    if (writeDrawCommand) std::cout << "fTree->Draw(\"" << drawStr << "\", \"" <<  cutStr << "\", \"" << drawOptionsStr << "\");" << std::endl;
 263    return fTree->Draw(drawStr.Data(), cutStr.Data(), drawOptionsStr.Data());
 264 }
 265
 266
 267 Int_t AliTPCCalibViewer::EasyDraw(const char* drawCommand, Int_t sector, const char* cuts, const char* drawOptions, Bool_t writeDrawCommand) const {
 268   //
 269   // easy drawing of data, use '~' for abbreviation of '.fElements'
 270   // example: EasyDraw("CETmean~-CETmean_mean", 34, "(CETmean~-CETmean_mean)>0")
 271   // sector: sector-number - only the specified sector will be drwawn
 272   // cuts: specifies cuts
 273   // drawOptions: draw options like 'same'
 274   // writeDrawCommand: write the command, that is passed to TTree::Draw
 275   //
 276    if (sector >= 0 && sector < 72) {
 277       char sectorChr[3];
 278       sprintf(sectorChr, "%i", sector);
 279       return EasyDraw(drawCommand, sectorChr, cuts, drawOptions, writeDrawCommand);
 280    }
 281    Error("EasyDraw","The TPC contains only sectors between 0 and 71.");
 282    return -1;
 283 }
 284
 285
 286 //_____________________________________________________________________________
 287 Int_t AliTPCCalibViewer::EasyDraw1D(const char* drawCommand, const char* sector, const char* cuts, const char* drawOptions, Bool_t writeDrawCommand) const {
 288   //
 289   // easy drawing of data, use '~' for abbreviation of '.fElements'
 290   // example: EasyDraw("CETmean~-CETmean_mean", "A", "(CETmean~-CETmean_mean)>0")
 291   // sector: sector-number - the specified sector will be drwawn
 292   //         'A'/'C' or 'a'/'c' - side A/C will be drawn
 293   //         'ALL' - whole TPC will be drawn, projected on one side
 294   // cuts: specifies cuts
 295   // drawOptions: draw options like 'same'
 296   // writeDrawCommand: write the command, that is passed to TTree::Draw
 297   //
 298
 299    TString drawStr(drawCommand);
 300    TString sectorStr(sector);
 301    TString drawOptionsStr(drawOptions);
 302    sectorStr.ToUpper();
 303    TString cutStr("");
 304
 305    if (sectorStr == "A")
 306       cutStr += "(sector/18)%2==0 ";
 307    else if  (sectorStr == "C")
 308       cutStr += "(sector/18)%2==1 ";
 309    else if (sectorStr.IsDigit()) {
 310       Int_t isec = sectorStr.Atoi();
 311       if (isec < 0 || isec > 71) {
 312          Error("EasyDraw","The TPC contains only sectors between 0 and 71.");
 313          return -1;
 314       }
 315       cutStr += "(sector==";
 316       cutStr += isec;
 317       cutStr += ") ";
 318    }
 319
 320    if (cuts && cuts[0] != 0) {
 321       if (cutStr.Length() != 0) cutStr += "&& ";
 322       cutStr += "(";
 323       cutStr += cuts;
 324       cutStr += ")";
 325    }
 326
 327    drawStr.ReplaceAll("~", ".fElements");
 328    cutStr.ReplaceAll("~", ".fElements");
 329    if (writeDrawCommand) std::cout << "fTree->Draw(\"" << drawStr << "\", \"" <<  cutStr << "\", \"" << drawOptionsStr << "\");" << std::endl;
 330    return fTree->Draw(drawStr.Data(), cutStr.Data(), drawOptionsStr.Data());
 331 }
 332
 333
 334 Int_t AliTPCCalibViewer::EasyDraw1D(const char* drawCommand, Int_t sector, const char* cuts, const char* drawOptions, Bool_t writeDrawCommand) const {
 335   //
 336   // easy drawing of data, use '~' for abbreviation of '.fElements'
 337   // example: EasyDraw("CETmean~-CETmean_mean", 34, "(CETmean~-CETmean_mean)>0")
 338   // sector: sector-number - the specified sector will be drwawn
 339   // cuts: specifies cuts
 340   // drawOptions: draw options like 'same'
 341   // writeDrawCommand: write the command, that is passed to TTree::Draw
 342   //
 343
 344    if (sector >= 0 && sector < 72) {
 345       char sectorChr[3];
 346       sprintf(sectorChr, "%i", sector);
 347       return EasyDraw1D(drawCommand, sectorChr, cuts, drawOptions, writeDrawCommand);
 348    }
 349   Error("EasyDraw","The TPC contains only sectors between 0 and 71.");
 350   return -1;
 351 }
 352
 353
 354 Int_t  AliTPCCalibViewer::DrawHisto1D(const char* drawCommand, Int_t sector, const char* cuts, const char *sigmas, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM) const {
 355    //
 356    // Easy drawing of data, in principle the same as EasyDraw1D
 357    // Difference: A line for the mean / median / LTM is drawn
 358    // in 'sigmas' you can specify in which distance to the mean/median/LTM you want to see a line in sigma-units, separated by ';'
 359    // 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.
 360    // "plotMean", "plotMedian" and "plotLTM": what kind of lines do you want to see?
 361    //
 362    if (sector >= 0 && sector < 72) {
 363       char sectorChr[3];
 364       sprintf(sectorChr, "%i", sector);
 365       return DrawHisto1D(drawCommand, sectorChr, cuts, sigmas, plotMean, plotMedian, plotLTM);
 366    }
 367    Error("DrawHisto1D","The TPC contains only sectors between 0 and 71.");
 368    return -1;
 369 }
 370
 371
 372 Int_t  AliTPCCalibViewer::DrawHisto1D(const char* drawCommand, const char* sector, const char* cuts, const char *sigmas, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM) const {
 373    //
 374    // Easy drawing of data, in principle the same as EasyDraw1D
 375    // Difference: A line for the mean / median / LTM is drawn
 376    // in 'sigmas' you can specify in which distance to the mean/median/LTM you want to see a line in sigma-units, separated by ';'
 377    // 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.
 378    // "plotMean", "plotMedian" and "plotLTM": what kind of lines do you want to see?
 379    //
 380    Int_t oldOptStat = gStyle->GetOptStat();
 381    gStyle->SetOptStat(0000000);
 382    Double_t ltmFraction = 0.8;
 383
 384    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");
 385    TVectorF nsigma(sigmasTokens->GetEntriesFast());
 386    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
 387       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
 388       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
 389       nsigma[i] = sig;
 390    }
 391
 392    TString drawStr(drawCommand);
 393    drawStr += " >> tempHist";
 394    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts);
 395    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
 396    // FIXME is this histogram deleted automatically?
 397    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
 398
 399    Double_t mean = TMath::Mean(entries, values);
 400    Double_t median = TMath::Median(entries, values);
 401    Double_t sigma = TMath::RMS(entries, values);
 402    Double_t maxY = htemp->GetMaximum();
 403
 404    char c[500];
 405    TLegend * legend = new TLegend(.7,.7, .99, .99, "Statistical information");
 406 //    sprintf(c, "%s, sector: %i", type, sector);
 407    fListOfObjectsToBeDeleted->Add(legend);
 408
 409    if (plotMean) {
 410       // draw Mean
 411       TLine* line = new TLine(mean, 0, mean, maxY);
 412       fListOfObjectsToBeDeleted->Add(line);
 413       line->SetLineColor(kRed);
 414       line->SetLineWidth(2);
 415       line->SetLineStyle(1);
 416       line->Draw();
 417       sprintf(c, "Mean: %f", mean);
 418       legend->AddEntry(line, c, "l");
 419       // draw sigma lines
 420       for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
 421          TLine* linePlusSigma = new TLine(mean + nsigma[i] * sigma, 0, mean + nsigma[i] * sigma, maxY);
 422          fListOfObjectsToBeDeleted->Add(linePlusSigma);
 423          linePlusSigma->SetLineColor(kRed);
 424          linePlusSigma->SetLineStyle(2 + i);
 425          linePlusSigma->Draw();
 426          TLine* lineMinusSigma = new TLine(mean - nsigma[i] * sigma, 0, mean - nsigma[i] * sigma, maxY);
 427          fListOfObjectsToBeDeleted->Add(lineMinusSigma);
 428          lineMinusSigma->SetLineColor(kRed);
 429          lineMinusSigma->SetLineStyle(2 + i);
 430          lineMinusSigma->Draw();
 431          sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
 432          legend->AddEntry(lineMinusSigma, c, "l");
 433       }
 434    }
 435    if (plotMedian) {
 436       // draw median
 437       TLine* line = new TLine(median, 0, median, maxY);
 438       fListOfObjectsToBeDeleted->Add(line);
 439       line->SetLineColor(kBlue);
 440       line->SetLineWidth(2);
 441       line->SetLineStyle(1);
 442       line->Draw();
 443       sprintf(c, "Median: %f", median);
 444       legend->AddEntry(line, c, "l");
 445       // draw sigma lines
 446       for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
 447          TLine* linePlusSigma = new TLine(median + nsigma[i] * sigma, 0, median + nsigma[i]*sigma, maxY);
 448          fListOfObjectsToBeDeleted->Add(linePlusSigma);
 449          linePlusSigma->SetLineColor(kBlue);
 450          linePlusSigma->SetLineStyle(2 + i);
 451          linePlusSigma->Draw();
 452          TLine* lineMinusSigma = new TLine(median - nsigma[i] * sigma, 0, median - nsigma[i]*sigma, maxY);
 453          fListOfObjectsToBeDeleted->Add(lineMinusSigma);
 454          lineMinusSigma->SetLineColor(kBlue);
 455          lineMinusSigma->SetLineStyle(2 + i);
 456          lineMinusSigma->Draw();
 457          sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
 458          legend->AddEntry(lineMinusSigma, c, "l");
 459       }
 460    }
 461    if (plotLTM) {
 462       // draw LTM
 463       Double_t ltmRms = 0;
 464       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
 465       TLine* line = new TLine(ltm, 0, ltm, maxY);
 466       fListOfObjectsToBeDeleted->Add(line);
 467       line->SetLineColor(kGreen+2);
 468       line->SetLineWidth(2);
 469       line->SetLineStyle(1);
 470       line->Draw();
 471       sprintf(c, "LTM: %f", ltm);
 472       legend->AddEntry(line, c, "l");
 473       // draw sigma lines
 474       for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
 475          TLine* linePlusSigma = new TLine(ltm + nsigma[i] * ltmRms, 0, ltm + nsigma[i] * ltmRms, maxY);
 476          fListOfObjectsToBeDeleted->Add(linePlusSigma);
 477          linePlusSigma->SetLineColor(kGreen+2);
 478          linePlusSigma->SetLineStyle(2+i);
 479          linePlusSigma->Draw();
 480
 481          TLine* lineMinusSigma = new TLine(ltm - nsigma[i] * ltmRms, 0, ltm - nsigma[i] * ltmRms, maxY);
 482          fListOfObjectsToBeDeleted->Add(lineMinusSigma);
 483          lineMinusSigma->SetLineColor(kGreen+2);
 484          lineMinusSigma->SetLineStyle(2+i);
 485          lineMinusSigma->Draw();
 486          sprintf(c, "%i #sigma = %f", (Int_t)(nsigma[i]), (Float_t)(nsigma[i] * ltmRms));
 487          legend->AddEntry(lineMinusSigma, c, "l");
 488       }
 489    }
 490    if (!plotMean && !plotMedian && !plotLTM) return -1;
 491    legend->Draw();
 492    gStyle->SetOptStat(oldOptStat);
 493    return 1;
 494 }
 495
 496
 497 Int_t AliTPCCalibViewer::SigmaCut(const char* drawCommand, Int_t sector, const char* cuts, Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, Bool_t pm, const char *sigmas, Float_t sigmaStep) const {
 498    //
 499    // 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
 500    // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'array', 'n' specifies the length of the array
 501    // '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
 502    // 'nbins': number of bins, 'binLow': first bin, 'binUp': last bin
 503    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
 504    // sigmaStep: the binsize of the generated histogram
 505    // Begin_Latex
 506    // 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 }
 507    // End_Latex
 508    //
 509    //
 510    // Creates a histogram, where you can see, how much of the data are inside sigma-intervals
 511    // around the mean/median/LTM
 512    // with drawCommand, sector and cuts you specify your input data, see EasyDraw
 513    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma)
 514    // sigmaStep: the binsize of the generated histogram
 515    // plotMean/plotMedian/plotLTM: specifies where to put the center
 516    //
 517    if (sector >= 0 && sector < 72) {
 518       char sectorChr[3];
 519       sprintf(sectorChr, "%i", sector);
 520       return SigmaCut(drawCommand, sectorChr, cuts, sigmaMax, plotMean, plotMedian, plotLTM, pm, sigmas, sigmaStep);
 521    }
 522    Error("SigmaCut","The TPC contains only sectors between 0 and 71.");
 523    return -1;
 524 }
 525
 526
 527 Int_t AliTPCCalibViewer::SigmaCut(const char* drawCommand, const char* sector, const char* cuts, Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, Bool_t pm, const char *sigmas, Float_t sigmaStep) const {
 528    //
 529    // Creates a histogram, where you can see, how much of the data are inside sigma-intervals
 530    // around the mean/median/LTM
 531    // with drawCommand, sector and cuts you specify your input data, see EasyDraw
 532    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma)
 533    // sigmaStep: the binsize of the generated histogram
 534    // plotMean/plotMedian/plotLTM: specifies where to put the center
 535    //
 536
 537    Double_t ltmFraction = 0.8;
 538
 539    TString drawStr(drawCommand);
 540    drawStr += " >> tempHist";
 541
 542    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
 543    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
 544    // FIXME is this histogram deleted automatically?
 545    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
 546
 547    Double_t mean = TMath::Mean(entries, values);
 548    Double_t median = TMath::Median(entries, values);
 549    Double_t sigma = TMath::RMS(entries, values);
 550
 551    TLegend * legend = new TLegend(.7,.7, .99, .99, "Cumulative");
 552    fListOfObjectsToBeDeleted->Add(legend);
 553    TH1F *cutHistoMean = 0;
 554    TH1F *cutHistoMedian = 0;
 555    TH1F *cutHistoLTM = 0;
 556
 557    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");
 558    TVectorF nsigma(sigmasTokens->GetEntriesFast());
 559    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
 560       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
 561       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
 562       nsigma[i] = sig;
 563    }
 564
 565    if (plotMean) {
 566       cutHistoMean = AliTPCCalibViewer::SigmaCut(htemp, mean, sigma, sigmaMax, sigmaStep, pm);
 567       if (cutHistoMean) {
 568          fListOfObjectsToBeDeleted->Add(cutHistoMean);
 569          cutHistoMean->SetLineColor(kRed);
 570          legend->AddEntry(cutHistoMean, "Mean", "l");
 571          cutHistoMean->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 572          cutHistoMean->Draw();
 573          DrawLines(cutHistoMean, nsigma, legend, kRed, pm);
 574       } // if (cutHistoMean)
 575
 576    }
 577    if (plotMedian) {
 578       cutHistoMedian = AliTPCCalibViewer::SigmaCut(htemp, median, sigma, sigmaMax, sigmaStep, pm);
 579       if (cutHistoMedian) {
 580          fListOfObjectsToBeDeleted->Add(cutHistoMedian);
 581          cutHistoMedian->SetLineColor(kBlue);
 582          legend->AddEntry(cutHistoMedian, "Median", "l");
 583          cutHistoMedian->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 584          if (plotMean && cutHistoMean) cutHistoMedian->Draw("same");
 585             else cutHistoMedian->Draw();
 586          DrawLines(cutHistoMedian, nsigma, legend, kBlue, pm);
 587       }  // if (cutHistoMedian)
 588    }
 589    if (plotLTM) {
 590       Double_t ltmRms = 0;
 591       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
 592       cutHistoLTM = AliTPCCalibViewer::SigmaCut(htemp, ltm, ltmRms, sigmaMax, sigmaStep, pm);
 593       if (cutHistoLTM) {
 594          fListOfObjectsToBeDeleted->Add(cutHistoLTM);
 595          cutHistoLTM->SetLineColor(kGreen+2);
 596          legend->AddEntry(cutHistoLTM, "LTM", "l");
 597          cutHistoLTM->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 598          if (plotMean && cutHistoMean || plotMedian && cutHistoMedian) cutHistoLTM->Draw("same");
 599             else cutHistoLTM->Draw();
 600          DrawLines(cutHistoLTM, nsigma, legend, kGreen+2, pm);
 601       }
 602    }
 603    if (!plotMean && !plotMedian && !plotLTM) return -1;
 604    legend->Draw();
 605    return 1;
 606 }
 607
 608 Int_t AliTPCCalibViewer::SigmaCutNew(const char* drawCommand, const char* sector, const char* cuts, Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, Bool_t pm, const char *sigmas, Float_t sigmaStep) const {
 609    //
 610    // Creates a histogram, where you can see, how much of the data are inside sigma-intervals
 611    // around the mean/median/LTM
 612    // with drawCommand, sector and cuts you specify your input data, see EasyDraw
 613    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma)
 614    // sigmaStep: the binsize of the generated histogram
 615    // plotMean/plotMedian/plotLTM: specifies where to put the center
 616    //
 617
 618    // Double_t ltmFraction = 0.8;  //unused
 619    // avoid compiler warnings:
 620    sigmaMax = sigmaMax;
 621    pm = pm;
 622    sigmaStep = sigmaStep;
 623
 624    TString drawStr(drawCommand);
 625    drawStr += " >> tempHist";
 626
 627    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
 628    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
 629    TGraph *cutGraphMean   = 0;
 630    // TGraph *cutGraphMedian = 0;
 631    // TGraph *cutGraphLTM    = 0;
 632    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
 633    Int_t    *index  = new Int_t[entries];
 634    Float_t  *xarray = new Float_t[entries];
 635    Float_t  *yarray = new Float_t[entries];
 636    TMath::Sort(entries, values, index, kFALSE);
 637
 638    Double_t mean = TMath::Mean(entries, values);
 639    // Double_t median = TMath::Median(entries, values);
 640    Double_t sigma = TMath::RMS(entries, values);
 641
 642    TLegend * legend = new TLegend(.7,.7, .99, .99, "Cumulative");
 643    fListOfObjectsToBeDeleted->Add(legend);
 644
 645    // parse sigmas string
 646    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");
 647    TVectorF nsigma(sigmasTokens->GetEntriesFast());
 648    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
 649       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
 650       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
 651       nsigma[i] = sig;
 652    }
 653
 654    if (plotMean) {
 655       for (Int_t i = 0; i < entries; i++) {
 656          xarray[i] = TMath::Abs(values[index[i]] - mean) / sigma;
 657          yarray[i] = float(i) / float(entries);
 658       }
 659       cutGraphMean = new TGraph(entries, xarray, yarray);
 660       if (cutGraphMean) {
 661          fListOfObjectsToBeDeleted->Add(cutGraphMean);
 662          cutGraphMean->SetLineColor(kRed);
 663          legend->AddEntry(cutGraphMean, "Mean", "l");
 664          cutGraphMean->SetTitle(Form("%s, Cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 665          cutGraphMean->Draw("alu");
 666          DrawLines(cutGraphMean, nsigma, legend, kRed, kTRUE);
 667       }
 668    }
 669    /*
 670    if (plotMedian) {
 671       cutHistoMedian = AliTPCCalibViewer::SigmaCut(htemp, median, sigma, sigmaMax, sigmaStep, pm);
 672       if (cutHistoMedian) {
 673          fListOfObjectsToBeDeleted->Add(cutHistoMedian);
 674          cutHistoMedian->SetLineColor(kBlue);
 675          legend->AddEntry(cutHistoMedian, "Median", "l");
 676          cutHistoMedian->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 677          if (plotMean && cutHistoMean) cutHistoMedian->Draw("same");
 678             else cutHistoMedian->Draw();
 679          DrawLines(cutHistoMedian, nsigma, legend, kBlue, pm);
 680       }  // if (cutHistoMedian)
 681    }
 682    if (plotLTM) {
 683       Double_t ltmRms = 0;
 684       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
 685       cutHistoLTM = AliTPCCalibViewer::SigmaCut(htemp, ltm, ltmRms, sigmaMax, sigmaStep, pm);
 686       if (cutHistoLTM) {
 687          fListOfObjectsToBeDeleted->Add(cutHistoLTM);
 688          cutHistoLTM->SetLineColor(kGreen+2);
 689          legend->AddEntry(cutHistoLTM, "LTM", "l");
 690          cutHistoLTM->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 691          if (plotMean && cutHistoMean || plotMedian && cutHistoMedian) cutHistoLTM->Draw("same");
 692             else cutHistoLTM->Draw();
 693          DrawLines(cutHistoLTM, nsigma, legend, kGreen+2, pm);
 694       }
 695    }*/
 696    if (!plotMean && !plotMedian && !plotLTM) return -1;
 697    legend->Draw();
 698    return 1;
 699 }
 700
 701
 702
 703
 704 Int_t AliTPCCalibViewer::Integrate(const char* drawCommand,       Int_t sector, const char* cuts, Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, const char *sigmas, Float_t sigmaStep) const {
 705    //
 706    // 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"
 707    // "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
 708    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
 709    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
 710    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
 711    // The actual work is done on the array.
 712    /* Begin_Latex
 713          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 }
 714       End_Latex
 715    */
 716    if (sector >= 0 && sector < 72) {
 717       char sectorChr[3];
 718       sprintf(sectorChr, "%i", sector);
 719       return Integrate(drawCommand, sectorChr, cuts, sigmaMax, plotMean, plotMedian, plotLTM, sigmas, sigmaStep);
 720    }
 721    Error("Integrate","The TPC contains only sectors between 0 and 71.");
 722    return -1;
 723
 724 }
 725
 726
 727 Int_t AliTPCCalibViewer::IntegrateOld(const char* drawCommand, const char* sector, const char* cuts, Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, const char *sigmas, Float_t sigmaStep) const {
 728    //
 729    // 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"
 730    // "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
 731    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
 732    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
 733    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
 734    // The actual work is done on the array.
 735    /* Begin_Latex
 736          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 }
 737       End_Latex
 738    */
 739
 740    Double_t ltmFraction = 0.8;
 741
 742    TString drawStr(drawCommand);
 743    drawStr += " >> tempHist";
 744
 745    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
 746    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
 747    // FIXME is this histogram deleted automatically?
 748    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
 749
 750    Double_t mean = TMath::Mean(entries, values);
 751    Double_t median = TMath::Median(entries, values);
 752    Double_t sigma = TMath::RMS(entries, values);
 753
 754    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");
 755    TVectorF nsigma(sigmasTokens->GetEntriesFast());
 756    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
 757       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
 758       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
 759       nsigma[i] = sig;
 760    }
 761
 762    TLegend * legend = new TLegend(.7,.7, .99, .99, "Integrated histogram");
 763    fListOfObjectsToBeDeleted->Add(legend);
 764    TH1F *integralHistoMean = 0;
 765    TH1F *integralHistoMedian = 0;
 766    TH1F *integralHistoLTM = 0;
 767
 768    if (plotMean) {
 769       integralHistoMean = AliTPCCalibViewer::Integrate(htemp, mean, sigma, sigmaMax, sigmaStep);
 770       if (integralHistoMean) {
 771          fListOfObjectsToBeDeleted->Add(integralHistoMean);
 772          integralHistoMean->SetLineColor(kRed);
 773          legend->AddEntry(integralHistoMean, "Mean", "l");
 774          integralHistoMean->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 775          integralHistoMean->Draw();
 776          DrawLines(integralHistoMean, nsigma, legend, kRed, kTRUE);
 777       }
 778    }
 779    if (plotMedian) {
 780       integralHistoMedian = AliTPCCalibViewer::Integrate(htemp, median, sigma, sigmaMax, sigmaStep);
 781       if (integralHistoMedian) {
 782          fListOfObjectsToBeDeleted->Add(integralHistoMedian);
 783          integralHistoMedian->SetLineColor(kBlue);
 784          legend->AddEntry(integralHistoMedian, "Median", "l");
 785          integralHistoMedian->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 786          if (plotMean && integralHistoMean) integralHistoMedian->Draw("same");
 787             else integralHistoMedian->Draw();
 788          DrawLines(integralHistoMedian, nsigma, legend, kBlue, kTRUE);
 789       }
 790    }
 791    if (plotLTM) {
 792       Double_t ltmRms = 0;
 793       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
 794       integralHistoLTM = AliTPCCalibViewer::Integrate(htemp, ltm, ltmRms, sigmaMax, sigmaStep);
 795       if (integralHistoLTM) {
 796          fListOfObjectsToBeDeleted->Add(integralHistoLTM);
 797          integralHistoLTM->SetLineColor(kGreen+2);
 798          legend->AddEntry(integralHistoLTM, "LTM", "l");
 799          integralHistoLTM->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 800          if (plotMean && integralHistoMean || plotMedian && integralHistoMedian) integralHistoLTM->Draw("same");
 801             else integralHistoLTM->Draw();
 802          DrawLines(integralHistoLTM, nsigma, legend, kGreen+2, kTRUE);
 803       }
 804    }
 805    if (!plotMean && !plotMedian && !plotLTM) return -1;
 806    legend->Draw();
 807    return 1;
 808 }
 809
 810
 811 Int_t AliTPCCalibViewer::Integrate(const char* drawCommand, const char* sector, const char* cuts, Float_t sigmaMax, Bool_t plotMean, Bool_t plotMedian, Bool_t plotLTM, const char *sigmas, Float_t sigmaStep) const {
 812    //
 813    // 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"
 814    // "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
 815    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
 816    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
 817    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
 818    // The actual work is done on the array.
 819    /* Begin_Latex
 820          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 }
 821       End_Latex
 822    */
 823
 824    Double_t ltmFraction = 0.8;
 825    // avoid compiler warnings:
 826    sigmaMax = sigmaMax;
 827    sigmaStep = sigmaStep;
 828
 829    TString drawStr(drawCommand);
 830    drawStr += " >> tempHist";
 831
 832    Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
 833    TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
 834    TGraph *integralGraphMean   = 0;
 835    TGraph *integralGraphMedian = 0;
 836    TGraph *integralGraphLTM    = 0;
 837    Double_t *values = fTree->GetV1();  // value is the array containing 'entries' numbers
 838    Int_t    *index  = new Int_t[entries];
 839    Float_t  *xarray = new Float_t[entries];
 840    Float_t  *yarray = new Float_t[entries];
 841    TMath::Sort(entries, values, index, kFALSE);
 842
 843    Double_t mean = TMath::Mean(entries, values);
 844    Double_t median = TMath::Median(entries, values);
 845    Double_t sigma = TMath::RMS(entries, values);
 846
 847    // parse sigmas string
 848    TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");
 849    TVectorF nsigma(sigmasTokens->GetEntriesFast());
 850    for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
 851       TString str(((TObjString*)sigmasTokens->At(i))->GetString());
 852       Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
 853       nsigma[i] = sig;
 854    }
 855
 856    TLegend * legend = new TLegend(.7,.7, .99, .99, "Integrated histogram");
 857    fListOfObjectsToBeDeleted->Add(legend);
 858
 859    if (plotMean) {
 860       for (Int_t i = 0; i < entries; i++) {
 861          xarray[i] = (values[index[i]] - mean) / sigma;
 862          yarray[i] = float(i) / float(entries);
 863       }
 864       integralGraphMean = new TGraph(entries, xarray, yarray);
 865       if (integralGraphMean) {
 866          fListOfObjectsToBeDeleted->Add(integralGraphMean);
 867          integralGraphMean->SetLineColor(kRed);
 868          legend->AddEntry(integralGraphMean, "Mean", "l");
 869          integralGraphMean->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 870          integralGraphMean->Draw("alu");
 871          DrawLines(integralGraphMean, nsigma, legend, kRed, kTRUE);
 872       }
 873    }
 874    if (plotMedian) {
 875       for (Int_t i = 0; i < entries; i++) {
 876          xarray[i] = (values[index[i]] - median) / sigma;
 877          yarray[i] = float(i) / float(entries);
 878       }
 879       integralGraphMedian = new TGraph(entries, xarray, yarray);
 880       if (integralGraphMedian) {
 881          fListOfObjectsToBeDeleted->Add(integralGraphMedian);
 882          integralGraphMedian->SetLineColor(kBlue);
 883          legend->AddEntry(integralGraphMedian, "Median", "l");
 884          integralGraphMedian->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 885          if (plotMean && integralGraphMean) integralGraphMedian->Draw("samelu");
 886             else integralGraphMedian->Draw("alu");
 887          DrawLines(integralGraphMedian, nsigma, legend, kBlue, kTRUE);
 888       }
 889    }
 890    if (plotLTM) {
 891       Double_t ltmRms = 0;
 892       Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
 893       for (Int_t i = 0; i < entries; i++) {
 894          xarray[i] = (values[index[i]] - ltm) / ltmRms;
 895          yarray[i] = float(i) / float(entries);
 896       }
 897       integralGraphLTM = new TGraph(entries, xarray, yarray);
 898       if (integralGraphLTM) {
 899          fListOfObjectsToBeDeleted->Add(integralGraphLTM);
 900          integralGraphLTM->SetLineColor(kGreen+2);
 901          legend->AddEntry(integralGraphLTM, "LTM", "l");
 902          integralGraphLTM->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
 903          if (plotMean && integralGraphMean || plotMedian && integralGraphMedian) integralGraphLTM->Draw("samelu");
 904             else integralGraphLTM->Draw("alu");
 905          DrawLines(integralGraphLTM, nsigma, legend, kGreen+2, kTRUE);
 906       }
 907    }
 908    if (!plotMean && !plotMedian && !plotLTM) return -1;
 909    legend->Draw();
 910    return entries;
 911 }
 912
 913
 914 void AliTPCCalibViewer::DrawLines(TH1F *histogram, TVectorF nsigma, TLegend *legend, Int_t color, Bool_t pm) const {
 915    //
 916    // Private function for SigmaCut(...) and Integrate(...)
 917    // Draws lines into the given histogram, specified by "nsigma", the lines are addeed to the legend
 918    //
 919
 920    // start to draw the lines, loop over requested sigmas
 921    char c[500];
 922    for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
 923       if (!pm) {
 924          Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
 925          TLine* lineUp = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
 926          fListOfObjectsToBeDeleted->Add(lineUp);
 927          lineUp->SetLineColor(color);
 928          lineUp->SetLineStyle(2 + i);
 929          lineUp->Draw();
 930          TLine* lineLeft = new TLine(nsigma[i], histogram->GetBinContent(bin), 0, histogram->GetBinContent(bin));
 931          fListOfObjectsToBeDeleted->Add(lineLeft);
 932          lineLeft->SetLineColor(color);
 933          lineLeft->SetLineStyle(2 + i);
 934          lineLeft->Draw();
 935          sprintf(c, "Fraction(%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
 936          legend->AddEntry(lineLeft, c, "l");
 937       }
 938       else { // if (pm)
 939          Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
 940          TLine* lineUp1 = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
 941          fListOfObjectsToBeDeleted->Add(lineUp1);
 942          lineUp1->SetLineColor(color);
 943          lineUp1->SetLineStyle(2 + i);
 944          lineUp1->Draw();
 945          TLine* lineLeft1 = new TLine(nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
 946          fListOfObjectsToBeDeleted->Add(lineLeft1);
 947          lineLeft1->SetLineColor(color);
 948          lineLeft1->SetLineStyle(2 + i);
 949          lineLeft1->Draw();
 950          sprintf(c, "Fraction(+%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
 951          legend->AddEntry(lineLeft1, c, "l");
 952          bin = histogram->GetXaxis()->FindBin(-nsigma[i]);
 953          TLine* lineUp2 = new TLine(-nsigma[i], 0, -nsigma[i], histogram->GetBinContent(bin));
 954          fListOfObjectsToBeDeleted->Add(lineUp2);
 955          lineUp2->SetLineColor(color);
 956          lineUp2->SetLineStyle(2 + i);
 957          lineUp2->Draw();
 958          TLine* lineLeft2 = new TLine(-nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
 959          fListOfObjectsToBeDeleted->Add(lineLeft2);
 960          lineLeft2->SetLineColor(color);
 961          lineLeft2->SetLineStyle(2 + i);
 962          lineLeft2->Draw();
 963          sprintf(c, "Fraction(-%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
 964          legend->AddEntry(lineLeft2, c, "l");
 965       }
 966    }  // for (Int_t i = 0; i < nsigma.GetNoElements(); i++)
 967 }
 968
 969
 970 void AliTPCCalibViewer::DrawLines(TGraph *graph, TVectorF nsigma, TLegend *legend, Int_t color, Bool_t pm) const {
 971    //
 972    // Private function for SigmaCut(...) and Integrate(...)
 973    // Draws lines into the given histogram, specified by "nsigma", the lines are addeed to the legend
 974    //
 975
 976    // start to draw the lines, loop over requested sigmas
 977    char c[500];
 978    for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
 979       if (!pm) {
 980          TLine* lineUp = new TLine(nsigma[i], 0, nsigma[i], graph->Eval(nsigma[i]));
 981          fListOfObjectsToBeDeleted->Add(lineUp);
 982          lineUp->SetLineColor(color);
 983          lineUp->SetLineStyle(2 + i);
 984          lineUp->Draw();
 985          TLine* lineLeft = new TLine(nsigma[i], graph->Eval(nsigma[i]), 0, graph->Eval(nsigma[i]));
 986          fListOfObjectsToBeDeleted->Add(lineLeft);
 987          lineLeft->SetLineColor(color);
 988          lineLeft->SetLineStyle(2 + i);
 989          lineLeft->Draw();
 990          sprintf(c, "Fraction(%f #sigma) = %f",nsigma[i], graph->Eval(nsigma[i]));
 991          legend->AddEntry(lineLeft, c, "l");
 992       }
 993       else { // if (pm)
 994          TLine* lineUp1 = new TLine(nsigma[i], 0, nsigma[i], graph->Eval(nsigma[i]));
 995          fListOfObjectsToBeDeleted->Add(lineUp1);
 996          lineUp1->SetLineColor(color);
 997          lineUp1->SetLineStyle(2 + i);
 998          lineUp1->Draw();
 999          TLine* lineLeft1 = new TLine(nsigma[i], graph->Eval(nsigma[i]), graph->GetHistogram()->GetXaxis()->GetBinLowEdge(0), graph->Eval(nsigma[i]));
1000          fListOfObjectsToBeDeleted->Add(lineLeft1);
1001          lineLeft1->SetLineColor(color);
1002          lineLeft1->SetLineStyle(2 + i);
1003          lineLeft1->Draw();
1004          sprintf(c, "Fraction(+%f #sigma) = %f",nsigma[i], graph->Eval(nsigma[i]));
1005          legend->AddEntry(lineLeft1, c, "l");
1006          TLine* lineUp2 = new TLine(-nsigma[i], 0, -nsigma[i], graph->Eval(-nsigma[i]));
1007          fListOfObjectsToBeDeleted->Add(lineUp2);
1008          lineUp2->SetLineColor(color);
1009          lineUp2->SetLineStyle(2 + i);
1010          lineUp2->Draw();
1011          TLine* lineLeft2 = new TLine(-nsigma[i], graph->Eval(-nsigma[i]), graph->GetHistogram()->GetXaxis()->GetBinLowEdge(0), graph->Eval(-nsigma[i]));
1012          fListOfObjectsToBeDeleted->Add(lineLeft2);
1013          lineLeft2->SetLineColor(color);
1014          lineLeft2->SetLineStyle(2 + i);
1015          lineLeft2->Draw();
1016          sprintf(c, "Fraction(-%f #sigma) = %f",nsigma[i], graph->Eval(-nsigma[i]));
1017          legend->AddEntry(lineLeft2, c, "l");
1018       }
1019    }  // for (Int_t i = 0; i < nsigma.GetNoElements(); i++)
1020 }
1021
1022
1023
1024
1025
1026 /////////////////
1027 // Array tools //
1028 /////////////////
1029
1030
1031 Int_t AliTPCCalibViewer::GetBin(Float_t value, Int_t nbins, Double_t binLow, Double_t binUp){
1032    // Returns the 'bin' for 'value'
1033    // The interval between 'binLow' and 'binUp' is divided into 'nbins' equidistant bins
1034    // avoid index out of bounds error: 'if (bin < binLow) bin = binLow' and vice versa
1035    /* Begin_Latex
1036          GetBin(value) = #frac{nbins - 1}{binUp - binLow} #upoint (value - binLow) +1
1037       End_Latex
1038    */
1039
1040    Int_t bin =  TMath::Nint( (Float_t)(value - binLow) / (Float_t)(binUp - binLow) * (nbins-1) ) + 1;
1041    // avoid index out of bounds:
1042    if (value < binLow) bin = 0;
1043    if (value > binUp)  bin = nbins + 1;
1044    return bin;
1045
1046 }
1047
1048
1049 Double_t AliTPCCalibViewer::GetLTM(Int_t n, Double_t *array, Double_t *sigma, Double_t fraction){
1050    //
1051    //  returns the LTM and sigma
1052    //
1053    Double_t *ddata = new Double_t[n];
1054    Double_t mean = 0, lsigma = 0;
1055    UInt_t nPoints = 0;
1056    for (UInt_t i = 0; i < (UInt_t)n; i++) {
1057          ddata[nPoints]= array[nPoints];
1058          nPoints++;
1059    }
1060    Int_t hh = TMath::Min(TMath::Nint(fraction * nPoints), Int_t(n));
1061    AliMathBase::EvaluateUni(nPoints, ddata, mean, lsigma, hh);
1062    if (sigma) *sigma = lsigma;
1063    delete [] ddata;
1064    return mean;
1065 }
1066
1067
1068 TH1F* AliTPCCalibViewer::SigmaCut(TH1F *histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep, Bool_t pm) {
1069    //
1070    // 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
1071    // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'histogram'
1072    // '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
1073    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
1074    // sigmaStep: the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
1075    // pm: Decide weather Begin_Latex t > 0 End_Latex (first case) or Begin_Latex t End_Latex arbitrary (secound case)
1076    // The actual work is done on the array.
1077    /* Begin_Latex
1078          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
1079          or
1080          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 }
1081       End_Latex
1082       begin_macro(source)
1083       {
1084          Float_t mean = 0;
1085          Float_t sigma = 1.5;
1086          Float_t sigmaMax = 4;
1087          gROOT->SetStyle("Plain");
1088          TH1F *distribution = new TH1F("Distribution1", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
1089          TRandom rand(23);
1090          for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
1091          Float_t *ar = distribution->GetArray();
1092
1093          TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas_SigmaCut", "", 350, 350);
1094          macro_example_canvas->Divide(0,3);
1095          TVirtualPad *pad1 = macro_example_canvas->cd(1);
1096          pad1->SetGridy();
1097          pad1->SetGridx();
1098          distribution->Draw();
1099          TVirtualPad *pad2 = macro_example_canvas->cd(2);
1100          pad2->SetGridy();
1101          pad2->SetGridx();
1102
1103          TH1F *shist = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax);
1104          shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
1105          shist->Draw();
1106          TVirtualPad *pad3 = macro_example_canvas->cd(3);
1107          pad3->SetGridy();
1108          pad3->SetGridx();
1109          TH1F *shistPM = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax, -1, kTRUE);
1110          shistPM->Draw();
1111          return macro_example_canvas;
1112       }
1113       end_macro
1114    */
1115
1116    Float_t *array = histogram->GetArray();
1117    Int_t    nbins = histogram->GetXaxis()->GetNbins();
1118    Float_t binLow = histogram->GetXaxis()->GetXmin();
1119    Float_t binUp  = histogram->GetXaxis()->GetXmax();
1120    return AliTPCCalibViewer::SigmaCut(nbins, array, mean, sigma, nbins, binLow, binUp, sigmaMax, sigmaStep, pm);
1121 }
1122
1123
1124
1125 TH1F* AliTPCCalibViewer::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){
1126    //
1127    // 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
1128    // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'array', 'n' specifies the length of the array
1129    // '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
1130    // 'nbins': number of bins, 'binLow': first bin, 'binUp': last bin
1131    // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
1132    // sigmaStep: the binsize of the generated histogram
1133    // Here the actual work is done.
1134
1135    if (sigma == 0) return 0;
1136    Float_t binWidth = (binUp-binLow)/(nbins - 1);
1137    if (sigmaStep <= 0) sigmaStep = binWidth;
1138    Int_t kbins = (Int_t)(sigmaMax * sigma / sigmaStep) + 1; // + 1  due to overflow bin in histograms
1139    if (pm) kbins = 2 * (Int_t)(sigmaMax * sigma / sigmaStep) + 1;
1140    Float_t kbinLow = !pm ? 0 : -sigmaMax;
1141    Float_t kbinUp  = sigmaMax;
1142    TH1F *hist = new TH1F("sigmaCutHisto","Cumulative; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp);
1143    hist->SetDirectory(0);
1144    hist->Reset();
1145
1146    // calculate normalization
1147    Double_t normalization = 0;
1148    for (Int_t i = 0; i <= n; i++) {
1149         normalization += array[i];
1150    }
1151
1152    // given units: units from given histogram
1153    // sigma units: in units of sigma
1154    // iDelta: integrate in interval (mean +- iDelta), given units
1155    // x:      ofset from mean for integration, given units
1156    // hist:   needs
1157
1158 //    printf("nbins: %i, binLow: %f, binUp: %f \n", nbins, binLow, binUp);
1159    // fill histogram
1160    for (Float_t iDelta = 0; iDelta <= sigmaMax * sigma; iDelta += sigmaStep) {
1161       // integrate array
1162       Double_t valueP = array[GetBin(mean, nbins, binLow, binUp)];
1163       Double_t valueM = array[GetBin(mean-binWidth, nbins, binLow, binUp)];
1164       // add bin of mean value only once to the histogram
1165 //       printf("++ adding bins: ");
1166       for (Float_t x = binWidth; x <= iDelta; x += binWidth) {
1167          valueP += (mean + x <= binUp)  ? array[GetBin(mean + x, nbins, binLow, binUp)] : 0;
1168          valueM += (mean-binWidth - x >= binLow) ? array[GetBin(mean-binWidth - x, nbins, binLow, binUp)] : 0;
1169 //          printf("%i, ", GetBin(mean + x, nbins, binLow, binUp));
1170       }
1171 //       printf("\n");
1172       if (valueP / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", valueP, normalization);
1173       if (valueP / normalization > 100) return hist;
1174       if (valueM / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", valueM, normalization);
1175       if (valueM / normalization > 100) return hist;
1176       valueP = (valueP / normalization);
1177       valueM = (valueM / normalization);
1178       if (pm) {
1179          Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
1180          hist->SetBinContent(bin, valueP);
1181          bin = GetBin(-iDelta/sigma, kbins, kbinLow, kbinUp);
1182          hist->SetBinContent(bin, valueM);
1183       }
1184       else { // if (!pm)
1185          Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
1186          hist->SetBinContent(bin, valueP + valueM);
1187 //          printf("  first integration bin: %i, last integration bin in + direction: %i \n", GetBin(mean+binWidth, nbins, binLow, binUp), GetBin(iDelta, nbins, binLow, binUp));
1188 //          printf("  first integration bin: %i, last integration bin in - direction: %i \n", GetBin(mean+binWidth, nbins, binLow, binUp), GetBin(-iDelta, nbins, binLow, binUp));
1189 //          printf("  value: %f, normalization: %f, iDelta: %f, Bin: %i \n", valueP+valueM, normalization, iDelta, bin);
1190       }
1191    }
1192    //hist->SetMaximum(0.7);
1193    if (!pm) hist->SetMaximum(1.2);
1194    return hist;
1195 }
1196
1197
1198 TH1F* AliTPCCalibViewer::SigmaCut(Int_t n, Double_t *array, Double_t mean, Double_t sigma, Int_t nbins, Double_t *xbins, Double_t sigmaMax){
1199    //
1200    // SigmaCut for variable binsize
1201    // NOT YET IMPLEMENTED !!!
1202    //
1203    printf("SigmaCut with variable binsize, Not yet implemented\n");
1204    // avoid compiler warnings:
1205    n=n;
1206    mean=mean;
1207    sigma=sigma;
1208    nbins=nbins;
1209    sigmaMax=sigmaMax;
1210    array=array;
1211    xbins=xbins;
1212
1213    return 0;
1214 }
1215
1216
1217 TH1F* AliTPCCalibViewer::Integrate(TH1F *histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
1218    //
1219    // 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"
1220    // "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
1221    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
1222    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
1223    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
1224    // The actual work is done on the array.
1225    /* Begin_Latex
1226          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 }
1227       End_Latex
1228       begin_macro(source)
1229       {
1230          Float_t mean = 0;
1231          Float_t sigma = 1.5;
1232          Float_t sigmaMax = 4;
1233          gROOT->SetStyle("Plain");
1234          TH1F *distribution = new TH1F("Distribution2", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
1235          TRandom rand(23);
1236          for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
1237          Float_t *ar = distribution->GetArray();
1238
1239          TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas_Integrate", "", 350, 350);
1240          macro_example_canvas->Divide(0,2);
1241          TVirtualPad *pad1 = macro_example_canvas->cd(1);
1242          pad1->SetGridy();
1243          pad1->SetGridx();
1244          distribution->Draw();
1245          TVirtualPad *pad2 = macro_example_canvas->cd(2);
1246          pad2->SetGridy();
1247          pad2->SetGridx();
1248          TH1F *shist = AliTPCCalibViewer::Integrate(distribution, mean, sigma, sigmaMax);
1249          shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
1250          shist->Draw();
1251
1252          return macro_example_canvas_Integrate;
1253       }
1254       end_macro
1255    */
1256
1257
1258    Float_t *array = histogram->GetArray();
1259    Int_t    nbins = histogram->GetXaxis()->GetNbins();
1260    Float_t binLow = histogram->GetXaxis()->GetXmin();
1261    Float_t binUp  = histogram->GetXaxis()->GetXmax();
1262    return AliTPCCalibViewer::Integrate(nbins, array, nbins, binLow, binUp, mean, sigma, sigmaMax, sigmaStep);
1263 }
1264
1265
1266 TH1F* AliTPCCalibViewer::Integrate(Int_t n, Float_t *array, Int_t nbins, Float_t binLow, Float_t binUp, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
1267    // 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"
1268    // "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
1269    // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
1270    // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
1271    // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
1272    // Here the actual work is done.
1273
1274    Bool_t givenUnits = kTRUE;
1275    if (sigma != 0 && sigmaMax != 0) givenUnits = kFALSE;
1276    if (givenUnits) {
1277       sigma = 1;
1278       sigmaMax = (binUp - binLow) / 2.;
1279    }
1280
1281    Float_t binWidth = (binUp-binLow)/(nbins - 1);
1282    if (sigmaStep <= 0) sigmaStep = binWidth;
1283    Int_t kbins =  (Int_t)(sigmaMax * sigma / sigmaStep) + 1;  // + 1  due to overflow bin in histograms
1284    Float_t kbinLow = givenUnits ? binLow : -sigmaMax;
1285    Float_t kbinUp  = givenUnits ? binUp  : sigmaMax;
1286    TH1F *hist = 0;
1287    if (givenUnits)  hist = new TH1F("integratedHisto","Integrated Histogram; Given x; Fraction of included data", kbins, kbinLow, kbinUp);
1288    if (!givenUnits) hist = new TH1F("integratedHisto","Integrated Histogram; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp);
1289    hist->SetDirectory(0);
1290    hist->Reset();
1291
1292    // calculate normalization
1293  //  printf("calculating normalization, integrating from bin 1 to %i \n", n);
1294    Double_t normalization = 0;
1295    for (Int_t i = 1; i <= n; i++) {
1296         normalization += array[i];
1297    }
1298  //  printf("normalization: %f \n", normalization);
1299
1300    // given units: units from given histogram
1301    // sigma units: in units of sigma
1302    // iDelta: integrate in interval (mean +- iDelta), given units
1303    // x:      ofset from mean for integration, given units
1304    // hist:   needs
1305
1306    // fill histogram
1307    for (Float_t iDelta = mean - sigmaMax * sigma; iDelta <= mean + sigmaMax * sigma; iDelta += sigmaStep) {
1308       // integrate array
1309       Double_t value = 0;
1310       for (Float_t x = mean - sigmaMax * sigma; x <= iDelta; x += binWidth) {
1311          value += (x <= binUp && x >= binLow)  ? array[GetBin(x, nbins, binLow, binUp)] : 0;
1312       }
1313       if (value / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail  +++ \n", value, normalization);
1314       if (value / normalization > 100) return hist;
1315       Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
1316     //  printf("first integration bin: %i, last integration bin: %i \n", GetBin(mean - sigmaMax * sigma, nbins, binLow, binUp), GetBin(iDelta, nbins, binLow, binUp));
1317     //  printf("value: %f, normalization: %f, normalized value: %f, iDelta: %f, Bin: %i \n", value, normalization, value/normalization, iDelta, bin);
1318       value = (value / normalization);
1319       hist->SetBinContent(bin, value);
1320    }
1321    return hist;
1322 }
1323
1324
1325
1326
1327
1328 ////////////////////////
1329 // end of Array tools //
1330 ////////////////////////
1331
1332
1333
1334 //_____________________________________________________________________________
1335 AliTPCCalPad* AliTPCCalibViewer::GetCalPad(const char* desiredData, char* cuts, char* calPadName) const {
1336   //
1337   // creates a AliTPCCalPad out of the 'desiredData'
1338   // the functionality of EasyDraw1D is used
1339   // calPadName specifies the name of the created AliTPCCalPad
1340   //  - this takes a while -
1341   //
1342    TString drawStr(desiredData);
1343    drawStr.Append(":channel~");
1344    AliTPCCalPad * createdCalPad = new AliTPCCalPad(calPadName, calPadName);
1345    Int_t entries = 0;
1346    for (Int_t sec = 0; sec < 72; sec++) {
1347       entries = EasyDraw1D(drawStr.Data(), (Int_t)sec, cuts, "goff");
1348       if (entries == -1) return 0;
1349       for (Int_t i = 0; i < entries; i++)
1350          createdCalPad->GetCalROC(sec)->SetValue((UInt_t)(fTree->GetV2()[i]), (Float_t)(fTree->GetV1()[i]));
1351    }
1352    return createdCalPad;
1353 }
1354
1355 //_____________________________________________________________________________
1356 AliTPCCalROC* AliTPCCalibViewer::GetCalROC(const char* desiredData, UInt_t sector, char* cuts) const {
1357   //
1358   // creates a AliTPCCalROC out of the desiredData
1359   // the functionality of EasyDraw1D is used
1360   // sector specifies the sector of the created AliTPCCalROC
1361   //
1362    TString drawStr(desiredData);
1363    drawStr.Append(":channel~");
1364    Int_t entries = EasyDraw1D(drawStr.Data(), (Int_t)sector, cuts, "goff");
1365    if (entries == -1) return 0;
1366    AliTPCCalROC * createdROC = new AliTPCCalROC(sector);
1367    for (Int_t i = 0; i < entries; i++)
1368       createdROC->SetValue((UInt_t)(fTree->GetV2()[i]), fTree->GetV1()[i]);
1369    return createdROC;
1370 }
1371
1372
1373 TObjArray* AliTPCCalibViewer::GetListOfVariables(Bool_t printList) {
1374   //
1375   // scan the tree  - produces a list of available variables in the tree
1376   // printList: print the list to the screen, after the scan is done
1377   //
1378    TObjArray* arr = new TObjArray();
1379    TObjString* str = 0;
1380    Int_t nentries = fTree->GetListOfBranches()->GetEntries();
1381    for (Int_t i = 0; i < nentries; i++) {
1382       str = new TObjString(fTree->GetListOfBranches()->At(i)->GetName());
1383       str->String().ReplaceAll("_Median", "");
1384       str->String().ReplaceAll("_Mean", "");
1385       str->String().ReplaceAll("_RMS", "");
1386       str->String().ReplaceAll("_LTM", "");
1387       str->String().ReplaceAll("_OutlierCutted", "");
1388       str->String().ReplaceAll(".", "");
1389       if (!arr->FindObject(str) &&
1390           !(str->String() == "channel" || str->String() == "gx" || str->String() == "gy" ||
1391             str->String() == "lx" || str->String() == "ly" || str->String() == "pad" ||
1392             str->String() == "row" || str->String() == "rpad" || str->String() == "sector"  ))
1393          arr->Add(str);
1394    }
1395    arr->Sort();
1396
1397    if (printList) {
1398       TIterator* iter = arr->MakeIterator();
1399       iter->Reset();
1400       TObjString* currentStr = 0;
1401       while ( (currentStr = (TObjString*)(iter->Next())) ) {
1402          std::cout << currentStr->GetString().Data() << std::endl;
1403       }
1404       delete iter;
1405    }
1406    return arr;
1407 }
1408
1409
1410 TObjArray* AliTPCCalibViewer::GetListOfNormalizationVariables(Bool_t printList) const{
1411   //
1412   // produces a list of available variables for normalization in the tree
1413   // printList: print the list to the screen, after the scan is done
1414   //
1415    TObjArray* arr = new TObjArray();
1416    arr->Add(new TObjString("_Mean"));
1417    arr->Add(new TObjString("_Mean_OutlierCutted"));
1418    arr->Add(new TObjString("_Median"));
1419    arr->Add(new TObjString("_Median_OutlierCutted"));
1420    arr->Add(new TObjString("_LTM"));
1421    arr->Add(new TObjString("_LTM_OutlierCutted"));
1422    arr->Add(new TObjString("LFitIntern_4_8.fElements"));
1423    arr->Add(new TObjString("GFitIntern_Lin.fElements"));
1424    arr->Add(new TObjString("GFitIntern_Par.fElements"));
1425    arr->Add(new TObjString("FitLinLocal"));
1426    arr->Add(new TObjString("FitLinGlobal"));
1427    arr->Add(new TObjString("FitParLocal"));
1428    arr->Add(new TObjString("FitParGlobal"));
1429
1430    if (printList) {
1431       TIterator* iter = arr->MakeIterator();
1432       iter->Reset();
1433       TObjString* currentStr = 0;
1434       while ((currentStr = (TObjString*)(iter->Next()))) {
1435          std::cout << currentStr->GetString().Data() << std::endl;
1436       }
1437       delete iter;
1438    }
1439    return arr;
1440 }
1441
1442
1443 TFriendElement* AliTPCCalibViewer::AddReferenceTree(const char* filename, const char* treename, const char* refname){
1444   //
1445   // add a reference tree to the current tree
1446   // by default the treename is 'calPads' and the reference treename is 'R'
1447   //
1448    TFile *file = new TFile(filename);
1449    fListOfObjectsToBeDeleted->Add(file);
1450    TTree * tree = (TTree*)file->Get(treename);
1451    return AddFriend(tree, refname);
1452 }
1453
1454
1455 TObjArray* AliTPCCalibViewer::GetArrayOfCalPads(){
1456   //
1457   // Returns a TObjArray with all AliTPCCalPads that are stored in the tree
1458   //  - this takes a while -
1459   //
1460    TObjArray *listOfCalPads = GetListOfVariables();
1461    TObjArray *calPadsArray = new TObjArray();
1462    Int_t numberOfCalPads = listOfCalPads->GetEntries();
1463    for (Int_t i = 0; i < numberOfCalPads; i++) {
1464      std::cout << "Creating calPad " << (i+1) << " of " << numberOfCalPads << "\r" << std::flush;
1465       char* calPadName = (char*)((TObjString*)(listOfCalPads->At(i)))->GetString().Data();
1466       TString drawCommand = ((TObjString*)(listOfCalPads->At(i)))->GetString();
1467       drawCommand.Append("~");
1468       AliTPCCalPad* calPad = GetCalPad(drawCommand.Data(), "", calPadName);
1469       calPadsArray->Add(calPad);
1470    }
1471    std::cout << std::endl;
1472    listOfCalPads->Delete();
1473    delete listOfCalPads;
1474    return calPadsArray;
1475 }
1476
1477
1478 TString* AliTPCCalibViewer::Fit(const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, TVectorD &fitParam, TMatrixD &covMatrix){
1479    //
1480    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
1481    // returns chi2, fitParam and covMatrix
1482    // returns TString with fitted formula
1483    //
1484
1485    TString formulaStr(formula);
1486    TString drawStr(drawCommand);
1487    TString cutStr(cuts);
1488
1489    // abbreviations:
1490    drawStr.ReplaceAll("~",".fElements");
1491    cutStr.ReplaceAll("~",".fElements");
1492    formulaStr.ReplaceAll("~", ".fElements");
1493
1494    formulaStr.ReplaceAll("++", "~");
1495    TObjArray* formulaTokens = formulaStr.Tokenize("~");
1496    Int_t dim = formulaTokens->GetEntriesFast();
1497
1498    fitParam.ResizeTo(dim);
1499    covMatrix.ResizeTo(dim,dim);
1500
1501    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
1502    fitter->StoreData(kTRUE);
1503    fitter->ClearPoints();
1504
1505    Int_t entries = Draw(drawStr.Data(), cutStr.Data(), "goff");
1506    if (entries == -1) return new TString("An ERROR has occured during fitting!");
1507    Double_t **values = new Double_t*[dim+1] ;
1508
1509    for (Int_t i = 0; i < dim + 1; i++){
1510       Int_t centries = 0;
1511       if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff");
1512       else  centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff");
1513
1514       if (entries != centries) return new TString("An ERROR has occured during fitting!");
1515       values[i] = new Double_t[entries];
1516       memcpy(values[i],  fTree->GetV1(), entries*sizeof(Double_t));
1517    }
1518
1519    // add points to the fitter
1520    for (Int_t i = 0; i < entries; i++){
1521       Double_t x[1000];
1522       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
1523       fitter->AddPoint(x, values[dim][i], 1);
1524    }
1525
1526    fitter->Eval();
1527    fitter->GetParameters(fitParam);
1528    fitter->GetCovarianceMatrix(covMatrix);
1529    chi2 = fitter->GetChisquare();
1530    chi2 = chi2;
1531
1532    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
1533
1534    for (Int_t iparam = 0; iparam < dim; iparam++) {
1535      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
1536      if (iparam < dim-1) returnFormula.Append("+");
1537    }
1538    returnFormula.Append(" )");
1539    delete formulaTokens;
1540    delete fitter;
1541    delete[] values;
1542    return preturnFormula;
1543 }
1544
1545
1546 void AliTPCCalibViewer::MakeTreeWithObjects(const char * fileName, TObjArray * array, const char * mapFileName) {
1547   //
1548   // Write tree with all available information
1549   // im mapFileName is speciefied, the Map information are also written to the tree
1550   // AliTPCCalPad-Objects are written directly to the tree, so that they can be accessd later on
1551   // (does not work!!!)
1552   //
1553    AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
1554
1555    TObjArray* mapIROCs = 0;
1556    TObjArray* mapOROCs = 0;
1557    TVectorF *mapIROCArray = 0;
1558    TVectorF *mapOROCArray = 0;
1559    Int_t mapEntries = 0;
1560    TString* mapNames = 0;
1561
1562    if (mapFileName) {
1563       TFile mapFile(mapFileName, "read");
1564
1565       TList* listOfROCs = mapFile.GetListOfKeys();
1566       mapEntries = listOfROCs->GetEntries()/2;
1567       mapIROCs = new TObjArray(mapEntries*2);
1568       mapOROCs = new TObjArray(mapEntries*2);
1569       mapIROCArray = new TVectorF[mapEntries];
1570       mapOROCArray = new TVectorF[mapEntries];
1571
1572       mapNames = new TString[mapEntries];
1573       for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1574          TString rocName(((TKey*)(listOfROCs->At(ivalue*2)))->GetName());
1575          rocName.Remove(rocName.Length()-4, 4);
1576          mapIROCs->AddAt((AliTPCCalROC*)mapFile.Get((rocName + "IROC").Data()), ivalue);
1577          mapOROCs->AddAt((AliTPCCalROC*)mapFile.Get((rocName + "OROC").Data()), ivalue);
1578          mapNames[ivalue].Append(rocName);
1579       }
1580
1581       for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1582          mapIROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(0));
1583          mapOROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(36));
1584
1585          for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(0); ichannel++)
1586             (mapIROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapIROCs->At(ivalue)))->GetValue(ichannel);
1587          for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(36); ichannel++)
1588             (mapOROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapOROCs->At(ivalue)))->GetValue(ichannel);
1589       }
1590
1591    } //  if (mapFileName)
1592
1593    TTreeSRedirector cstream(fileName);
1594    Int_t arrayEntries = array->GetEntries();
1595
1596    // Read names of AliTPCCalPads and save them in names[]
1597    TString* names = new TString[arrayEntries];
1598    for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
1599       names[ivalue].Append(((AliTPCCalPad*)array->At(ivalue))->GetName());
1600
1601    for (UInt_t isector = 0; isector < tpcROCinstance->GetNSectors(); isector++) {
1602
1603       TVectorF *vectorArray = new TVectorF[arrayEntries];
1604       for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
1605          vectorArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
1606
1607
1608       //
1609       // fill vectors of variable per pad
1610       //
1611       TVectorF *posArray = new TVectorF[8];
1612       for (Int_t ivalue = 0; ivalue < 8; ivalue++)
1613          posArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
1614
1615       Float_t posG[3] = {0};
1616       Float_t posL[3] = {0};
1617       Int_t ichannel = 0;
1618       for (UInt_t irow = 0; irow < tpcROCinstance->GetNRows(isector); irow++) {
1619          for (UInt_t ipad = 0; ipad < tpcROCinstance->GetNPads(isector, irow); ipad++) {
1620             tpcROCinstance->GetPositionLocal(isector, irow, ipad, posL);
1621             tpcROCinstance->GetPositionGlobal(isector, irow, ipad, posG);
1622             posArray[0][ichannel] = irow;
1623             posArray[1][ichannel] = ipad;
1624             posArray[2][ichannel] = posL[0];
1625             posArray[3][ichannel] = posL[1];
1626             posArray[4][ichannel] = posG[0];
1627             posArray[5][ichannel] = posG[1];
1628             posArray[6][ichannel] = (Int_t)(ipad - (Double_t)(tpcROCinstance->GetNPads(isector, irow))/2);
1629             posArray[7][ichannel] = ichannel;
1630
1631             // loop over array containing AliTPCCalPads
1632             for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1633                AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue);
1634                AliTPCCalROC* calROC = calPad->GetCalROC(isector);
1635                if (calROC)
1636                   (vectorArray[ivalue])[ichannel] = calROC->GetValue(irow, ipad);
1637                else
1638                   (vectorArray[ivalue])[ichannel] = 0;
1639             }
1640             ichannel++;
1641          }
1642       }
1643       AliTPCCalROC dummyROC(0);
1644       for  (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1645          AliTPCCalROC *roc = ((AliTPCCalPad*)array->At(ivalue))->GetCalROC(isector);
1646          if (!roc) roc = &dummyROC;
1647          cstream << "calPads" <<
1648             (Char_t*)((names[ivalue] + ".=").Data()) << &vectorArray[ivalue];
1649          cstream << "calPads" <<
1650             (Char_t*)((names[ivalue] + "Pad.=").Data()) << roc;
1651       }
1652
1653       if (mapFileName) {
1654          for  (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1655             if (isector < 36)
1656                cstream << "calPads" <<
1657                   (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapIROCArray[ivalue];
1658             else
1659                cstream << "calPads" <<
1660                   (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapOROCArray[ivalue];
1661          }
1662       }
1663
1664       cstream << "calPads" <<
1665          "sector=" << isector;
1666
1667       cstream << "calPads" <<
1668          "row.=" << &posArray[0] <<
1669          "pad.=" << &posArray[1] <<
1670          "lx.=" << &posArray[2] <<
1671          "ly.=" << &posArray[3] <<
1672          "gx.=" << &posArray[4] <<
1673          "gy.=" << &posArray[5] <<
1674          "rpad.=" << &posArray[6] <<
1675          "channel.=" << &posArray[7];
1676
1677       cstream << "calPads" <<
1678          "\n";
1679
1680       delete[] posArray;
1681       delete[] vectorArray;
1682    } //for (UInt_t isector = 0; isector < tpcROCinstance->GetNSectors(); isector++)
1683
1684    delete[] names;
1685    if (mapFileName) {
1686       delete mapIROCs;
1687       delete mapOROCs;
1688       delete[] mapIROCArray;
1689       delete[] mapOROCArray;
1690       delete[] mapNames;
1691    }
1692 }
1693
1694
1695 void AliTPCCalibViewer::MakeTree(const char * fileName, TObjArray * array, const char * mapFileName, AliTPCCalPad* outlierPad, Float_t ltmFraction) {
1696   //
1697   // Write a tree with all available information
1698   // if mapFileName is speciefied, the Map information are also written to the tree
1699   // pads specified in outlierPad are not used for calculating statistics
1700   // The following statistical information on the basis of a ROC are calculated:
1701   // "_Median", "_Mean", "_LTM", "_RMS_LTM"
1702   // "_Median_OutlierCutted", "_Mean_OutlierCutted", "_RMS_OutlierCutted", "_LTM_OutlierCutted", "_RMS_LTM_OutlierCutted"
1703   // The following position variables are available:
1704   // "row", "pad", "lx", "ly", "gx", "gy", "rpad", "channel"
1705   //
1706   // The tree out of this function is the basis for the AliTPCCalibViewer and the AliTPCCalibViewerGUI.
1707
1708    AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
1709
1710    TObjArray* mapIROCs = 0;
1711    TObjArray* mapOROCs = 0;
1712    TVectorF *mapIROCArray = 0;
1713    TVectorF *mapOROCArray = 0;
1714    Int_t mapEntries = 0;
1715    TString* mapNames = 0;
1716
1717    if (mapFileName) {
1718       TFile mapFile(mapFileName, "read");
1719
1720       TList* listOfROCs = mapFile.GetListOfKeys();
1721       mapEntries = listOfROCs->GetEntries()/2;
1722       mapIROCs = new TObjArray(mapEntries*2);
1723       mapOROCs = new TObjArray(mapEntries*2);
1724       mapIROCArray = new TVectorF[mapEntries];
1725       mapOROCArray = new TVectorF[mapEntries];
1726
1727       mapNames = new TString[mapEntries];
1728       for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1729          TString rocName(((TKey*)(listOfROCs->At(ivalue*2)))->GetName());
1730          rocName.Remove(rocName.Length()-4, 4);
1731          mapIROCs->AddAt((AliTPCCalROC*)mapFile.Get((rocName + "IROC").Data()), ivalue);
1732          mapOROCs->AddAt((AliTPCCalROC*)mapFile.Get((rocName + "OROC").Data()), ivalue);
1733          mapNames[ivalue].Append(rocName);
1734       }
1735
1736       for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1737          mapIROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(0));
1738          mapOROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(36));
1739
1740          for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(0); ichannel++)
1741             (mapIROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapIROCs->At(ivalue)))->GetValue(ichannel);
1742          for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(36); ichannel++)
1743             (mapOROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapOROCs->At(ivalue)))->GetValue(ichannel);
1744       }
1745
1746    } //  if (mapFileName)
1747
1748    TTreeSRedirector cstream(fileName);
1749    Int_t arrayEntries = 0;
1750    if (array) arrayEntries = array->GetEntries();
1751
1752    TString* names = new TString[arrayEntries];
1753    for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
1754       names[ivalue].Append(((AliTPCCalPad*)array->At(ivalue))->GetName());
1755
1756    for (UInt_t isector = 0; isector < tpcROCinstance->GetNSectors(); isector++) {
1757       //
1758       // get statistic for given sector
1759       //
1760       TVectorF median(arrayEntries);
1761       TVectorF mean(arrayEntries);
1762       TVectorF rms(arrayEntries);
1763       TVectorF ltm(arrayEntries);
1764       TVectorF ltmrms(arrayEntries);
1765       TVectorF medianWithOut(arrayEntries);
1766       TVectorF meanWithOut(arrayEntries);
1767       TVectorF rmsWithOut(arrayEntries);
1768       TVectorF ltmWithOut(arrayEntries);
1769       TVectorF ltmrmsWithOut(arrayEntries);
1770
1771       TVectorF *vectorArray = new TVectorF[arrayEntries];
1772       for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
1773          vectorArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
1774
1775       for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1776          AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue);
1777          AliTPCCalROC* calROC = calPad->GetCalROC(isector);
1778          AliTPCCalROC* outlierROC = 0;
1779          if (outlierPad) outlierROC = outlierPad->GetCalROC(isector);
1780          if (calROC) {
1781             median[ivalue] = calROC->GetMedian();
1782             mean[ivalue] = calROC->GetMean();
1783             rms[ivalue] = calROC->GetRMS();
1784             Double_t ltmrmsValue = 0;
1785             ltm[ivalue] = calROC->GetLTM(&ltmrmsValue, ltmFraction);
1786             ltmrms[ivalue] = ltmrmsValue;
1787             if (outlierROC) {
1788                medianWithOut[ivalue] = calROC->GetMedian(outlierROC);
1789                meanWithOut[ivalue] = calROC->GetMean(outlierROC);
1790                rmsWithOut[ivalue] = calROC->GetRMS(outlierROC);
1791                ltmrmsValue = 0;
1792                ltmWithOut[ivalue] = calROC->GetLTM(&ltmrmsValue, ltmFraction, outlierROC);
1793                ltmrmsWithOut[ivalue] = ltmrmsValue;
1794             }
1795          }
1796          else {
1797             median[ivalue] = 0.;
1798             mean[ivalue] = 0.;
1799             rms[ivalue] = 0.;
1800             ltm[ivalue] = 0.;
1801             ltmrms[ivalue] = 0.;
1802             medianWithOut[ivalue] = 0.;
1803             meanWithOut[ivalue] = 0.;
1804             rmsWithOut[ivalue] = 0.;
1805             ltmWithOut[ivalue] = 0.;
1806             ltmrmsWithOut[ivalue] = 0.;
1807          }
1808       }
1809
1810       //
1811       // fill vectors of variable per pad
1812       //
1813       TVectorF *posArray = new TVectorF[8];
1814       for (Int_t ivalue = 0; ivalue < 8; ivalue++)
1815          posArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
1816
1817       Float_t posG[3] = {0};
1818       Float_t posL[3] = {0};
1819       Int_t ichannel = 0;
1820       for (UInt_t irow = 0; irow < tpcROCinstance->GetNRows(isector); irow++) {
1821          for (UInt_t ipad = 0; ipad < tpcROCinstance->GetNPads(isector, irow); ipad++) {
1822             tpcROCinstance->GetPositionLocal(isector, irow, ipad, posL);
1823             tpcROCinstance->GetPositionGlobal(isector, irow, ipad, posG);
1824             posArray[0][ichannel] = irow;
1825             posArray[1][ichannel] = ipad;
1826             posArray[2][ichannel] = posL[0];
1827             posArray[3][ichannel] = posL[1];
1828             posArray[4][ichannel] = posG[0];
1829             posArray[5][ichannel] = posG[1];
1830             posArray[6][ichannel] = (Int_t)(ipad - (Double_t)(tpcROCinstance->GetNPads(isector, irow))/2);
1831             posArray[7][ichannel] = ichannel;
1832
1833             // loop over array containing AliTPCCalPads
1834             for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1835                AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue);
1836                AliTPCCalROC* calROC = calPad->GetCalROC(isector);
1837                if (calROC)
1838                   (vectorArray[ivalue])[ichannel] = calROC->GetValue(irow, ipad);
1839                else
1840                   (vectorArray[ivalue])[ichannel] = 0;
1841             }
1842             ichannel++;
1843          }
1844       }
1845
1846       cstream << "calPads" <<
1847          "sector=" << isector;
1848
1849       for  (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1850          cstream << "calPads" <<
1851             (Char_t*)((names[ivalue] + "_Median=").Data()) << median[ivalue] <<
1852             (Char_t*)((names[ivalue] + "_Mean=").Data()) << mean[ivalue] <<
1853             (Char_t*)((names[ivalue] + "_RMS=").Data()) << rms[ivalue] <<
1854             (Char_t*)((names[ivalue] + "_LTM=").Data()) << ltm[ivalue] <<
1855             (Char_t*)((names[ivalue] + "_RMS_LTM=").Data()) << ltmrms[ivalue];
1856          if (outlierPad) {
1857             cstream << "calPads" <<
1858                (Char_t*)((names[ivalue] + "_Median_OutlierCutted=").Data()) << medianWithOut[ivalue] <<
1859                (Char_t*)((names[ivalue] + "_Mean_OutlierCutted=").Data()) << meanWithOut[ivalue] <<
1860                (Char_t*)((names[ivalue] + "_RMS_OutlierCutted=").Data()) << rmsWithOut[ivalue] <<
1861                (Char_t*)((names[ivalue] + "_LTM_OutlierCutted=").Data()) << ltmWithOut[ivalue] <<
1862                (Char_t*)((names[ivalue] + "_RMS_LTM_OutlierCutted=").Data()) << ltmrmsWithOut[ivalue];
1863          }
1864       }
1865
1866       for  (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1867          cstream << "calPads" <<
1868             (Char_t*)((names[ivalue] + ".=").Data()) << &vectorArray[ivalue];
1869       }
1870
1871       if (mapFileName) {
1872          for  (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1873             if (isector < 36)
1874                cstream << "calPads" <<
1875                   (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapIROCArray[ivalue];
1876             else
1877                cstream << "calPads" <<
1878                   (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapOROCArray[ivalue];
1879          }
1880       }
1881
1882       cstream << "calPads" <<
1883          "row.=" << &posArray[0] <<
1884          "pad.=" << &posArray[1] <<
1885          "lx.=" << &posArray[2] <<
1886          "ly.=" << &posArray[3] <<
1887          "gx.=" << &posArray[4] <<
1888          "gy.=" << &posArray[5] <<
1889          "rpad.=" << &posArray[6] <<
1890          "channel.=" << &posArray[7];
1891
1892       cstream << "calPads" <<
1893          "\n";
1894
1895       delete[] posArray;
1896       delete[] vectorArray;
1897    }
1898
1899
1900    delete[] names;
1901    if (mapFileName) {
1902       delete mapIROCs;
1903       delete mapOROCs;
1904       delete[] mapIROCArray;
1905       delete[] mapOROCArray;
1906       delete[] mapNames;
1907    }
1908 }
1909
1910
1911 void AliTPCCalibViewer::MakeTree(const char *outPutFileName, const Char_t *inputFileName, AliTPCCalPad *outlierPad, Float_t ltmFraction, const char *mapFileName ){
1912    //
1913    // Function to create a calibration Tree with all available information.
1914    // See also documentation to MakeTree
1915    // the file "inputFileName" must be in the following format (see also CreateObjectList):
1916    // (each colum separated by tabs, "dependingOnType" can have 2 or 3 colums)
1917    //
1918    // type      path    dependingOnType
1919    //
1920    // type == "CE":
1921    // dependingOnType = CETmean CEQmean CETrms
1922    //
1923    // type == "Pulser":
1924    // dependingOnType = PulserTmean     PulsterQmean    PulserTrms
1925    //
1926    // type == "Pedestals":
1927    // dependingOnType = Pedestals       Noise
1928    //
1929    // type == "CalPad":
1930    // dependingOnType = NameInFile      NameToWriteToFile
1931    //
1932    //
1933    TObjArray objArray;
1934    CreateObjectList(inputFileName, &objArray);
1935    MakeTree(outPutFileName, &objArray, mapFileName, outlierPad, ltmFraction);
1936 }
1937
1938 void AliTPCCalibViewer::CreateObjectList(const Char_t *filename, TObjArray *calibObjects){
1939    //
1940    // Function to create a TObjArray out of a given file
1941    // the file must be in the following format:
1942    // (each colum separated by tabs, "dependingOnType" can have 2 or 3 colums)
1943    //
1944    //
1945    // type      path    dependingOnType
1946    //
1947    // type == "CE":
1948    // dependingOnType = CETmean CEQmean CETrms
1949    //
1950    // type == "Pulser":
1951    // dependingOnType = PulserTmean     PulsterQmean    PulserTrms
1952    //
1953    // type == "Pedestals":
1954    // dependingOnType = Pedestals       Noise
1955    //
1956    // type == "CalPad":
1957    // dependingOnType = NameInFile      NameToWriteToFile
1958    //
1959    //
1960    //
1961    if ( calibObjects == 0x0 ) return;
1962    ifstream in;
1963    in.open(filename);
1964    if ( !in.is_open() ){
1965       fprintf(stderr,"Error: cannot open list file '%s'", filename);
1966       return;
1967    }
1968
1969    AliTPCCalPad *calPad=0x0;
1970
1971    TString sFile;
1972    sFile.ReadFile(in);
1973    in.close();
1974
1975    TObjArray *arrFileLine = sFile.Tokenize("\n");
1976    TIter nextLine(arrFileLine);
1977
1978    TObjString *sObjLine = 0x0;
1979    while ( (sObjLine = (TObjString*)nextLine()) ){
1980       TString sLine(sObjLine->GetString());
1981
1982       TObjArray *arrCol = sLine.Tokenize("\t");
1983       Int_t nCols = arrCol->GetEntriesFast();
1984
1985       TObjString *sObjType     = (TObjString*)(arrCol->At(0));
1986       TObjString *sObjFileName = (TObjString*)(arrCol->At(1));
1987       TObjString *sObjName = 0x0;
1988
1989       if ( !sObjType || !sObjFileName ) continue;
1990       TString sType(sObjType->GetString());
1991       TString sFileName(sObjFileName->GetString());
1992       printf("Type %s, opening %s \n", sType.Data(), sFileName.Data());
1993       TFile *fIn = TFile::Open(sFileName);
1994       if ( !fIn ){
1995          fprintf(stderr,"File not found: '%s'", sFileName.Data());
1996          continue;
1997       }
1998
1999       if ( sType == "CE" ){  // next three colums are the names for CETmean, CEQmean and CETrms
2000          AliTPCCalibCE *ce = (AliTPCCalibCE*)fIn->Get("AliTPCCalibCE");
2001          calPad = new AliTPCCalPad((TObjArray*)ce->GetCalPadT0());
2002          if (nCols > 2) {  // check, if the name is provided
2003             sObjName = (TObjString*)(arrCol->At(2));
2004             calPad->SetNameTitle(sObjName->GetString().Data(), sObjName->GetString().Data());
2005          }
2006          else calPad->SetNameTitle("CETmean","CETmean");
2007          calibObjects->Add(calPad);
2008
2009          calPad = new AliTPCCalPad((TObjArray*)ce->GetCalPadQ());
2010          if (nCols > 3) {  // check, if the name is provided
2011             sObjName = (TObjString*)(arrCol->At(3));
2012             calPad->SetNameTitle(sObjName->GetString().Data(), sObjName->GetString().Data());
2013          }
2014          else calPad->SetNameTitle("CEQmean","CEQmean");
2015          calibObjects->Add(calPad);
2016
2017          calPad = new AliTPCCalPad((TObjArray*)ce->GetCalPadRMS());
2018          if (nCols > 4) {  // check, if the name is provided
2019             sObjName = (TObjString*)(arrCol->At(4));
2020             calPad->SetNameTitle(sObjName->GetString().Data(), sObjName->GetString().Data());
2021          }
2022          else calPad->SetNameTitle("CETrms","CETrms");
2023          calibObjects->Add(calPad);
2024
2025       } else if ( sType == "Pulser") {
2026          AliTPCCalibPulser *sig = (AliTPCCalibPulser*)fIn->Get("AliTPCCalibPulser");
2027
2028          calPad = new AliTPCCalPad((TObjArray*)sig->GetCalPadT0());
2029          if (nCols > 2) {
2030             sObjName = (TObjString*)(arrCol->At(2));
2031             calPad->SetNameTitle(sObjName->GetString().Data(), sObjName->GetString().Data());
2032          }
2033          else calPad->SetNameTitle("PulserTmean","PulserTmean");
2034          calibObjects->Add(calPad);
2035
2036          calPad = new AliTPCCalPad((TObjArray*)sig->GetCalPadQ());
2037          if (nCols > 3) {
2038             sObjName = (TObjString*)(arrCol->At(3));
2039             calPad->SetNameTitle(sObjName->GetString().Data(), sObjName->GetString().Data());
2040          }
2041          else calPad->SetNameTitle("PulserQmean","PulserQmean");
2042          calibObjects->Add(calPad);
2043
2044          calPad = new AliTPCCalPad((TObjArray*)sig->GetCalPadRMS());
2045          if (nCols > 4) {
2046             sObjName = (TObjString*)(arrCol->At(4));
2047             calPad->SetNameTitle(sObjName->GetString().Data(), sObjName->GetString().Data());
2048          }
2049          else calPad->SetNameTitle("PulserTrms","PulserTrms");
2050          calibObjects->Add(calPad);
2051
2052       } else if ( sType == "Pedestals") {
2053          AliTPCCalibPedestal *ped = (AliTPCCalibPedestal*)fIn->Get("AliTPCCalibPedestal");
2054
2055          calPad = new AliTPCCalPad((TObjArray*)ped->GetCalPadPedestal());
2056          if (nCols > 2) {
2057             sObjName = (TObjString*)(arrCol->At(2));
2058             calPad->SetNameTitle(sObjName->GetString().Data(), sObjName->GetString().Data());
2059          }
2060          else calPad->SetNameTitle("Pedestals","Pedestals");
2061          calibObjects->Add(calPad);
2062
2063          calPad = new AliTPCCalPad((TObjArray*)ped->GetCalPadRMS());
2064          if (nCols > 3) {
2065             sObjName = (TObjString*)(arrCol->At(3));
2066             calPad->SetNameTitle(sObjName->GetString().Data(), sObjName->GetString().Data());
2067          }
2068          else calPad->SetNameTitle("Noise","Noise");
2069          calibObjects->Add(calPad);
2070
2071       } else if ( sType == "CalPad") {
2072          if (nCols > 2) sObjName = (TObjString*)(arrCol->At(2));
2073          else continue;
2074          calPad = new AliTPCCalPad(*(AliTPCCalPad*)fIn->Get(sObjName->GetString().Data()));
2075          if (nCols > 3) {
2076             sObjName = (TObjString*)(arrCol->At(3));
2077             calPad->SetNameTitle(sObjName->GetString().Data(), sObjName->GetString().Data());
2078          }
2079          calibObjects->Add(calPad);
2080       } else {
2081          fprintf(stderr,"Undefined Type: '%s'",sType.Data());
2082       }
2083       delete fIn;
2084    }
2085 }
2086
2087
2088