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