Adding README file + using proper Makefile (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibViewer.cxx
CommitLineData
39bcd65d 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 //
72d0ab7e 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// //
39bcd65d 34///////////////////////////////////////////////////////////////////////////////
35
36//
37// ROOT includes
38//
39#include <iostream>
4f3934a1 40#include <fstream>
39bcd65d 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>
72d0ab7e 48#include <TH1.h>
39bcd65d 49#include <TH1F.h>
50#include <THashTable.h>
51#include <TObjString.h>
52#include "TTreeStream.h"
53#include "TFile.h"
54#include "TKey.h"
4f3934a1 55#include "TGraph.h"
832bb12c 56#include "TDirectory.h"
4f3934a1 57#include "AliTPCCalibPulser.h"
58#include "AliTPCCalibPedestal.h"
59#include "AliTPCCalibCE.h"
832bb12c 60#include "TFriendElement.h"
4f3934a1 61// #include "TObjArray.h"
62// #include "TObjString.h"
63// #include "TString.h"
64// #include "AliTPCCalPad.h"
39bcd65d 65
66
67//
68// AliRoot includes
69//
70#include "AliTPCCalibViewer.h"
71
72ClassImp(AliTPCCalibViewer)
73
832bb12c 74
39bcd65d 75AliTPCCalibViewer::AliTPCCalibViewer()
76 :TObject(),
77 fTree(0),
78 fFile(0),
a6d2bd0c 79 fListOfObjectsToBeDeleted(0),
832bb12c 80 fTreeMustBeDeleted(0),
81 fAbbreviation(0),
82 fAppendString(0)
39bcd65d 83{
84 //
85 // Default constructor
86 //
87
88}
89
90//_____________________________________________________________________________
91AliTPCCalibViewer::AliTPCCalibViewer(const AliTPCCalibViewer &c)
92 :TObject(c),
93 fTree(0),
94 fFile(0),
a6d2bd0c 95 fListOfObjectsToBeDeleted(0),
832bb12c 96 fTreeMustBeDeleted(0),
97 fAbbreviation(0),
98 fAppendString(0)
39bcd65d 99{
100 //
101 // dummy AliTPCCalibViewer copy constructor
102 // not yet working!!!
103 //
104 fTree = c.fTree;
a6d2bd0c 105 fTreeMustBeDeleted = c.fTreeMustBeDeleted;
39bcd65d 106 //fFile = new TFile(*(c.fFile));
107 fListOfObjectsToBeDeleted = c.fListOfObjectsToBeDeleted;
832bb12c 108 fAbbreviation = c.fAbbreviation;
109 fAppendString = c.fAppendString;
39bcd65d 110}
111
112//_____________________________________________________________________________
113AliTPCCalibViewer::AliTPCCalibViewer(TTree* tree)
114 :TObject(),
115 fTree(0),
116 fFile(0),
a6d2bd0c 117 fListOfObjectsToBeDeleted(0),
832bb12c 118 fTreeMustBeDeleted(0),
119 fAbbreviation(0),
120 fAppendString(0)
39bcd65d 121{
122 //
123 // Constructor that initializes the calibration viewer
124 //
125 fTree = tree;
a6d2bd0c 126 fTreeMustBeDeleted = kFALSE;
39bcd65d 127 fListOfObjectsToBeDeleted = new TObjArray();
832bb12c 128 fAbbreviation = "~";
129 fAppendString = ".fElements";
39bcd65d 130}
131
132//_____________________________________________________________________________
133AliTPCCalibViewer::AliTPCCalibViewer(char* fileName, char* treeName)
134 :TObject(),
135 fTree(0),
136 fFile(0),
a6d2bd0c 137 fListOfObjectsToBeDeleted(0),
832bb12c 138 fTreeMustBeDeleted(0),
139 fAbbreviation(0),
140 fAppendString(0)
141
39bcd65d 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);
a6d2bd0c 149 fTreeMustBeDeleted = kTRUE;
39bcd65d 150 fListOfObjectsToBeDeleted = new TObjArray();
832bb12c 151 fAbbreviation = "~";
152 fAppendString = ".fElements";
39bcd65d 153}
154
155//____________________________________________________________________________
156AliTPCCalibViewer & AliTPCCalibViewer::operator =(const AliTPCCalibViewer & param)
157{
158 //
159 // assignment operator - dummy
160 // not yet working!!!
161 //
162 fTree = param.fTree;
a6d2bd0c 163 fTreeMustBeDeleted = param.fTreeMustBeDeleted;
39bcd65d 164 //fFile = new TFile(*(param.fFile));
165 fListOfObjectsToBeDeleted = param.fListOfObjectsToBeDeleted;
832bb12c 166 fAbbreviation = param.fAbbreviation;
167 fAppendString = param.fAppendString;
39bcd65d 168 return (*this);
169}
170
171//_____________________________________________________________________________
172AliTPCCalibViewer::~AliTPCCalibViewer()
173{
174 //
175 // AliTPCCalibViewer destructor
a6d2bd0c 176 // all objects will be deleted, the file will be closed, the pictures will disappear
39bcd65d 177 //
a6d2bd0c 178 if (fTree && fTreeMustBeDeleted) {
179 fTree->SetCacheSize(0);
180 fTree->Delete();
181 //delete fTree;
182 }
39bcd65d 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//_____________________________________________________________________________
a6d2bd0c 196void 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//_____________________________________________________________________________
39bcd65d 214Int_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 //
72d0ab7e 225
39bcd65d 226 TString drawStr(drawCommand);
227 TString sectorStr(sector);
228 sectorStr.ToUpper();
229 TString cutStr("");
72d0ab7e 230 //TString drawOptionsStr("profcolz ");
231 TString drawOptionsStr("");
39bcd65d 232 TRandom rnd(0);
233 Int_t rndNumber = rnd.Integer(10000);
72d0ab7e 234
235 if (drawOptions && strcmp(drawOptions, "") != 0)
39bcd65d 236 drawOptionsStr += drawOptions;
72d0ab7e 237 else
238 drawOptionsStr += "profcolz";
39bcd65d 239
240 if (sectorStr == "A") {
832bb12c 241 drawStr += Form(":gy%s:gx%s>>prof", fAppendString.Data(), fAppendString.Data());
39bcd65d 242 drawStr += rndNumber;
243 drawStr += "(330,-250,250,330,-250,250)";
244 cutStr += "(sector/18)%2==0 ";
245 }
246 else if (sectorStr == "C") {
832bb12c 247 drawStr += Form(":gy%s:gx%s>>prof", fAppendString.Data(), fAppendString.Data());
39bcd65d 248 drawStr += rndNumber;
249 drawStr += "(330,-250,250,330,-250,250)";
250 cutStr += "(sector/18)%2==1 ";
251 }
252 else if (sectorStr == "ALL") {
832bb12c 253 drawStr += Form(":gy%s:gx%s>>prof", fAppendString.Data(), fAppendString.Data());
39bcd65d 254 drawStr += rndNumber;
255 drawStr += "(330,-250,250,330,-250,250)";
256 }
257 else if (sectorStr.IsDigit()) {
258 Int_t isec = sectorStr.Atoi();
832bb12c 259 drawStr += Form(":rpad%s:row%s>>prof", fAppendString.Data(), fAppendString.Data());
39bcd65d 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 }
832bb12c 280 drawStr.ReplaceAll(fAbbreviation, fAppendString);
281 cutStr.ReplaceAll(fAbbreviation, fAppendString);
39bcd65d 282 if (writeDrawCommand) std::cout << "fTree->Draw(\"" << drawStr << "\", \"" << cutStr << "\", \"" << drawOptionsStr << "\");" << std::endl;
832bb12c 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;
39bcd65d 289}
290
72d0ab7e 291
39bcd65d 292Int_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
72d0ab7e 310
39bcd65d 311//_____________________________________________________________________________
312Int_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
832bb12c 352 drawStr.ReplaceAll(fAbbreviation, fAppendString);
353 cutStr.ReplaceAll(fAbbreviation, fAppendString);
39bcd65d 354 if (writeDrawCommand) std::cout << "fTree->Draw(\"" << drawStr << "\", \"" << cutStr << "\", \"" << drawOptionsStr << "\");" << std::endl;
832bb12c 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;
39bcd65d 365}
366
72d0ab7e 367
39bcd65d 368Int_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
39bcd65d 387
832bb12c 388void 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
72d0ab7e 420Int_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}
39bcd65d 436
72d0ab7e 437
438Int_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 //
39bcd65d 446 Int_t oldOptStat = gStyle->GetOptStat();
447 gStyle->SetOptStat(0000000);
72d0ab7e 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;
39bcd65d 456 }
457
72d0ab7e 458 TString drawStr(drawCommand);
832bb12c 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 }
72d0ab7e 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
39bcd65d 475 char c[500];
72d0ab7e 476 TLegend * legend = new TLegend(.7,.7, .99, .99, "Statistical information");
477// sprintf(c, "%s, sector: %i", type, sector);
832bb12c 478 //fListOfObjectsToBeDeleted->Add(legend);
39bcd65d 479
39bcd65d 480 if (plotMean) {
72d0ab7e 481 // draw Mean
482 TLine* line = new TLine(mean, 0, mean, maxY);
832bb12c 483 //fListOfObjectsToBeDeleted->Add(line);
39bcd65d 484 line->SetLineColor(kRed);
485 line->SetLineWidth(2);
486 line->SetLineStyle(1);
487 line->Draw();
72d0ab7e 488 sprintf(c, "Mean: %f", mean);
39bcd65d 489 legend->AddEntry(line, c, "l");
72d0ab7e 490 // draw sigma lines
39bcd65d 491 for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
72d0ab7e 492 TLine* linePlusSigma = new TLine(mean + nsigma[i] * sigma, 0, mean + nsigma[i] * sigma, maxY);
832bb12c 493 //fListOfObjectsToBeDeleted->Add(linePlusSigma);
39bcd65d 494 linePlusSigma->SetLineColor(kRed);
72d0ab7e 495 linePlusSigma->SetLineStyle(2 + i);
39bcd65d 496 linePlusSigma->Draw();
72d0ab7e 497 TLine* lineMinusSigma = new TLine(mean - nsigma[i] * sigma, 0, mean - nsigma[i] * sigma, maxY);
832bb12c 498 //fListOfObjectsToBeDeleted->Add(lineMinusSigma);
39bcd65d 499 lineMinusSigma->SetLineColor(kRed);
72d0ab7e 500 lineMinusSigma->SetLineStyle(2 + i);
39bcd65d 501 lineMinusSigma->Draw();
72d0ab7e 502 sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
39bcd65d 503 legend->AddEntry(lineMinusSigma, c, "l");
504 }
505 }
39bcd65d 506 if (plotMedian) {
72d0ab7e 507 // draw median
508 TLine* line = new TLine(median, 0, median, maxY);
832bb12c 509 //fListOfObjectsToBeDeleted->Add(line);
39bcd65d 510 line->SetLineColor(kBlue);
511 line->SetLineWidth(2);
512 line->SetLineStyle(1);
513 line->Draw();
72d0ab7e 514 sprintf(c, "Median: %f", median);
39bcd65d 515 legend->AddEntry(line, c, "l");
72d0ab7e 516 // draw sigma lines
39bcd65d 517 for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
72d0ab7e 518 TLine* linePlusSigma = new TLine(median + nsigma[i] * sigma, 0, median + nsigma[i]*sigma, maxY);
832bb12c 519 //fListOfObjectsToBeDeleted->Add(linePlusSigma);
39bcd65d 520 linePlusSigma->SetLineColor(kBlue);
72d0ab7e 521 linePlusSigma->SetLineStyle(2 + i);
39bcd65d 522 linePlusSigma->Draw();
72d0ab7e 523 TLine* lineMinusSigma = new TLine(median - nsigma[i] * sigma, 0, median - nsigma[i]*sigma, maxY);
832bb12c 524 //fListOfObjectsToBeDeleted->Add(lineMinusSigma);
39bcd65d 525 lineMinusSigma->SetLineColor(kBlue);
72d0ab7e 526 lineMinusSigma->SetLineStyle(2 + i);
39bcd65d 527 lineMinusSigma->Draw();
72d0ab7e 528 sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
39bcd65d 529 legend->AddEntry(lineMinusSigma, c, "l");
530 }
531 }
39bcd65d 532 if (plotLTM) {
72d0ab7e 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);
832bb12c 537 //fListOfObjectsToBeDeleted->Add(line);
39bcd65d 538 line->SetLineColor(kGreen+2);
539 line->SetLineWidth(2);
540 line->SetLineStyle(1);
541 line->Draw();
72d0ab7e 542 sprintf(c, "LTM: %f", ltm);
39bcd65d 543 legend->AddEntry(line, c, "l");
72d0ab7e 544 // draw sigma lines
39bcd65d 545 for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
72d0ab7e 546 TLine* linePlusSigma = new TLine(ltm + nsigma[i] * ltmRms, 0, ltm + nsigma[i] * ltmRms, maxY);
832bb12c 547 //fListOfObjectsToBeDeleted->Add(linePlusSigma);
39bcd65d 548 linePlusSigma->SetLineColor(kGreen+2);
549 linePlusSigma->SetLineStyle(2+i);
550 linePlusSigma->Draw();
551
72d0ab7e 552 TLine* lineMinusSigma = new TLine(ltm - nsigma[i] * ltmRms, 0, ltm - nsigma[i] * ltmRms, maxY);
832bb12c 553 //fListOfObjectsToBeDeleted->Add(lineMinusSigma);
39bcd65d 554 lineMinusSigma->SetLineColor(kGreen+2);
555 lineMinusSigma->SetLineStyle(2+i);
556 lineMinusSigma->Draw();
72d0ab7e 557 sprintf(c, "%i #sigma = %f", (Int_t)(nsigma[i]), (Float_t)(nsigma[i] * ltmRms));
39bcd65d 558 legend->AddEntry(lineMinusSigma, c, "l");
559 }
560 }
72d0ab7e 561 if (!plotMean && !plotMedian && !plotLTM) return -1;
39bcd65d 562 legend->Draw();
563 gStyle->SetOptStat(oldOptStat);
72d0ab7e 564 return 1;
39bcd65d 565}
566
39bcd65d 567
72d0ab7e 568Int_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
598Int_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);
832bb12c 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 }
72d0ab7e 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");
832bb12c 628 //fListOfObjectsToBeDeleted->Add(legend);
72d0ab7e 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) {
832bb12c 644 //fListOfObjectsToBeDeleted->Add(cutHistoMean);
72d0ab7e 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) {
832bb12c 656 //fListOfObjectsToBeDeleted->Add(cutHistoMedian);
72d0ab7e 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) {
832bb12c 670 //fListOfObjectsToBeDeleted->Add(cutHistoLTM);
72d0ab7e 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
832bb12c 684
4f3934a1 685Int_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");
832bb12c 720 //fListOfObjectsToBeDeleted->Add(legend);
4f3934a1 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) {
832bb12c 738 //fListOfObjectsToBeDeleted->Add(cutGraphMean);
4f3934a1 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
72d0ab7e 778
72d0ab7e 779Int_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
4f3934a1 802Int_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 {
72d0ab7e 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");
832bb12c 838 //fListOfObjectsToBeDeleted->Add(legend);
72d0ab7e 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) {
832bb12c 846 //fListOfObjectsToBeDeleted->Add(integralHistoMean);
72d0ab7e 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);
39bcd65d 852 }
72d0ab7e 853 }
854 if (plotMedian) {
855 integralHistoMedian = AliTPCCalibViewer::Integrate(htemp, median, sigma, sigmaMax, sigmaStep);
856 if (integralHistoMedian) {
832bb12c 857 //fListOfObjectsToBeDeleted->Add(integralHistoMedian);
72d0ab7e 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);
39bcd65d 864 }
72d0ab7e 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) {
832bb12c 871 //fListOfObjectsToBeDeleted->Add(integralHistoLTM);
72d0ab7e 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);
39bcd65d 878 }
879 }
72d0ab7e 880 if (!plotMean && !plotMedian && !plotLTM) return -1;
881 legend->Draw();
882 return 1;
883}
884
885
4f3934a1 886Int_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);
832bb12c 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 }
4f3934a1 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");
832bb12c 937 //fListOfObjectsToBeDeleted->Add(legend);
4f3934a1 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) {
832bb12c 946 //fListOfObjectsToBeDeleted->Add(integralGraphMean);
4f3934a1 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) {
832bb12c 961 //fListOfObjectsToBeDeleted->Add(integralGraphMedian);
4f3934a1 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) {
832bb12c 979 //fListOfObjectsToBeDeleted->Add(integralGraphLTM);
4f3934a1 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
72d0ab7e 994void 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 //
39bcd65d 999
72d0ab7e 1000 // start to draw the lines, loop over requested sigmas
39bcd65d 1001 char c[500];
72d0ab7e 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));
832bb12c 1006 //fListOfObjectsToBeDeleted->Add(lineUp);
72d0ab7e 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));
832bb12c 1011 //fListOfObjectsToBeDeleted->Add(lineLeft);
72d0ab7e 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));
832bb12c 1021 //fListOfObjectsToBeDeleted->Add(lineUp1);
72d0ab7e 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));
832bb12c 1026 //fListOfObjectsToBeDeleted->Add(lineLeft1);
72d0ab7e 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));
832bb12c 1034 //fListOfObjectsToBeDeleted->Add(lineUp2);
72d0ab7e 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));
832bb12c 1039 //fListOfObjectsToBeDeleted->Add(lineLeft2);
72d0ab7e 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
4f3934a1 1050void 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]));
832bb12c 1061 //fListOfObjectsToBeDeleted->Add(lineUp);
4f3934a1 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]));
832bb12c 1066 //fListOfObjectsToBeDeleted->Add(lineLeft);
4f3934a1 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]));
832bb12c 1075 //fListOfObjectsToBeDeleted->Add(lineUp1);
4f3934a1 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]));
832bb12c 1080 //fListOfObjectsToBeDeleted->Add(lineLeft1);
4f3934a1 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]));
832bb12c 1087 //fListOfObjectsToBeDeleted->Add(lineUp2);
4f3934a1 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]));
832bb12c 1092 //fListOfObjectsToBeDeleted->Add(lineLeft2);
4f3934a1 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
72d0ab7e 1102
1103
1104
1105
1106/////////////////
1107// Array tools //
1108/////////////////
1109
1110
1111Int_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;
39bcd65d 1125
72d0ab7e 1126}
39bcd65d 1127
72d0ab7e 1128
1129Double_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;
39bcd65d 1145}
1146
1147
72d0ab7e 1148TH1F* 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");
5547c0ab 1168 TH1F *distribution = new TH1F("Distribution1", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
72d0ab7e 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
5547c0ab 1173 TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas_SigmaCut", "", 350, 350);
72d0ab7e 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
72d0ab7e 1204TH1F* 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
1277TH1F* 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
1296TH1F* 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");
5547c0ab 1313 TH1F *distribution = new TH1F("Distribution2", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
72d0ab7e 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
5547c0ab 1318 TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas_Integrate", "", 350, 350);
72d0ab7e 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
7072d276 1331 return macro_example_canvas_Integrate;
72d0ab7e 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
1345TH1F* 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
39bcd65d 1413//_____________________________________________________________________________
1414AliTPCCalPad* 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);
832bb12c 1422 drawStr.Append(":channel");
1423 drawStr.Append(fAbbreviation);
39bcd65d 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");
a6d2bd0c 1428 if (entries == -1) return 0;
39bcd65d 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//_____________________________________________________________________________
1436AliTPCCalROC* 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);
08eef367 1443 drawStr.Append(Form(":channel%s", fAbbreviation.Data()));
39bcd65d 1444 Int_t entries = EasyDraw1D(drawStr.Data(), (Int_t)sector, cuts, "goff");
a6d2bd0c 1445 if (entries == -1) return 0;
39bcd65d 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
1453TObjArray* 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 }
832bb12c 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
39bcd65d 1507 arr->Sort();
832bb12c 1508// ((TFriendElement*)gui->GetViewer()->GetTree()->GetListOfFriends()->At(0))->GetTree()->GetListOfBranches()->At(0)->GetName()
1509// ((TFriendElement*)gui->GetViewer()->GetTree()->GetListOfFriends()->At(0))->GetTree()->GetListOfBranches()
1510
39bcd65d 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
72d0ab7e 1525TObjArray* AliTPCCalibViewer::GetListOfNormalizationVariables(Bool_t printList) const{
39bcd65d 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"));
832bb12c 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())));
a6d2bd0c 1540 arr->Add(new TObjString("FitLinLocal"));
1541 arr->Add(new TObjString("FitLinGlobal"));
72d0ab7e 1542 arr->Add(new TObjString("FitParLocal"));
1543 arr->Add(new TObjString("FitParGlobal"));
39bcd65d 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
1558TFriendElement* 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
1570TObjArray* 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++) {
832bb12c 1579 std::cout << "Creating calPad " << (i+1) << " of " << numberOfCalPads << "\r" << std::flush;
39bcd65d 1580 char* calPadName = (char*)((TObjString*)(listOfCalPads->At(i)))->GetString().Data();
1581 TString drawCommand = ((TObjString*)(listOfCalPads->At(i)))->GetString();
832bb12c 1582 drawCommand.Append(fAbbreviation.Data());
39bcd65d 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
a6d2bd0c 1593TString* 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:
832bb12c 1605 drawStr.ReplaceAll(fAbbreviation, fAppendString);
1606 cutStr.ReplaceAll(fAbbreviation, fAppendString);
1607 formulaStr.ReplaceAll(fAbbreviation, fAppendString);
a6d2bd0c 1608
832bb12c 1609 formulaStr.ReplaceAll("++", fAbbreviation);
1610 TObjArray* formulaTokens = formulaStr.Tokenize(fAbbreviation.Data());
a6d2bd0c 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
39bcd65d 1661void 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++) {
72d0ab7e 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);
39bcd65d 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
72d0ab7e 1809
39bcd65d 1810void AliTPCCalibViewer::MakeTree(const char * fileName, TObjArray * array, const char * mapFileName, AliTPCCalPad* outlierPad, Float_t ltmFraction) {
1811 //
1812 // Write a tree with all available information
72d0ab7e 1813 // if mapFileName is speciefied, the Map information are also written to the tree
39bcd65d 1814 // pads specified in outlierPad are not used for calculating statistics
72d0ab7e 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
39bcd65d 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++) {
72d0ab7e 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);
39bcd65d 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);
72d0ab7e 1864 Int_t arrayEntries = 0;
1865 if (array) arrayEntries = array->GetEntries();
39bcd65d 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
4f3934a1 2025
2026void 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
832bb12c 2053
4f3934a1 2054void 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