]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCalibViewer.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removal of obsolete files from build system
[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>
40#include <TString.h>
41#include <TRandom.h>
42#include <TLegend.h>
43#include <TLine.h>
44#include <TCanvas.h>
45#include <TROOT.h>
46#include <TStyle.h>
72d0ab7e 47#include <TH1.h>
39bcd65d 48#include <TH1F.h>
49#include <THashTable.h>
50#include <TObjString.h>
51#include "TTreeStream.h"
52#include "TFile.h"
53#include "TKey.h"
54
55
56//
57// AliRoot includes
58//
59#include "AliTPCCalibViewer.h"
60
61ClassImp(AliTPCCalibViewer)
62
63AliTPCCalibViewer::AliTPCCalibViewer()
64 :TObject(),
65 fTree(0),
66 fFile(0),
a6d2bd0c 67 fListOfObjectsToBeDeleted(0),
68 fTreeMustBeDeleted(0)
39bcd65d 69{
70 //
71 // Default constructor
72 //
73
74}
75
76//_____________________________________________________________________________
77AliTPCCalibViewer::AliTPCCalibViewer(const AliTPCCalibViewer &c)
78 :TObject(c),
79 fTree(0),
80 fFile(0),
a6d2bd0c 81 fListOfObjectsToBeDeleted(0),
82 fTreeMustBeDeleted(0)
39bcd65d 83{
84 //
85 // dummy AliTPCCalibViewer copy constructor
86 // not yet working!!!
87 //
88 fTree = c.fTree;
a6d2bd0c 89 fTreeMustBeDeleted = c.fTreeMustBeDeleted;
39bcd65d 90 //fFile = new TFile(*(c.fFile));
91 fListOfObjectsToBeDeleted = c.fListOfObjectsToBeDeleted;
92}
93
94//_____________________________________________________________________________
95AliTPCCalibViewer::AliTPCCalibViewer(TTree* tree)
96 :TObject(),
97 fTree(0),
98 fFile(0),
a6d2bd0c 99 fListOfObjectsToBeDeleted(0),
100 fTreeMustBeDeleted(0)
39bcd65d 101{
102 //
103 // Constructor that initializes the calibration viewer
104 //
105 fTree = tree;
a6d2bd0c 106 fTreeMustBeDeleted = kFALSE;
39bcd65d 107 fListOfObjectsToBeDeleted = new TObjArray();
108}
109
110//_____________________________________________________________________________
111AliTPCCalibViewer::AliTPCCalibViewer(char* fileName, char* treeName)
112 :TObject(),
113 fTree(0),
114 fFile(0),
a6d2bd0c 115 fListOfObjectsToBeDeleted(0),
116 fTreeMustBeDeleted(0)
39bcd65d 117{
118 //
119 // Constructor to initialize the calibration viewer
120 // the file 'fileName' contains the tree 'treeName'
121 //
122 fFile = new TFile(fileName, "read");
123 fTree = (TTree*) fFile->Get(treeName);
a6d2bd0c 124 fTreeMustBeDeleted = kTRUE;
39bcd65d 125 fListOfObjectsToBeDeleted = new TObjArray();
126}
127
128//____________________________________________________________________________
129AliTPCCalibViewer & AliTPCCalibViewer::operator =(const AliTPCCalibViewer & param)
130{
131 //
132 // assignment operator - dummy
133 // not yet working!!!
134 //
135 fTree = param.fTree;
a6d2bd0c 136 fTreeMustBeDeleted = param.fTreeMustBeDeleted;
39bcd65d 137 //fFile = new TFile(*(param.fFile));
138 fListOfObjectsToBeDeleted = param.fListOfObjectsToBeDeleted;
139 return (*this);
140}
141
142//_____________________________________________________________________________
143AliTPCCalibViewer::~AliTPCCalibViewer()
144{
145 //
146 // AliTPCCalibViewer destructor
a6d2bd0c 147 // all objects will be deleted, the file will be closed, the pictures will disappear
39bcd65d 148 //
a6d2bd0c 149 if (fTree && fTreeMustBeDeleted) {
150 fTree->SetCacheSize(0);
151 fTree->Delete();
152 //delete fTree;
153 }
39bcd65d 154 if (fFile) {
155 fFile->Close();
156 fFile = 0;
157 }
158
159 for (Int_t i = fListOfObjectsToBeDeleted->GetEntriesFast()-1; i >= 0; i--) {
160 //cout << "Index " << i << " trying to delete the following object: " << fListOfObjectsToBeDeleted->At(i)->GetName() << "..."<< endl;
161 delete fListOfObjectsToBeDeleted->At(i);
162 }
163 delete fListOfObjectsToBeDeleted;
164}
165
a6d2bd0c 166//_____________________________________________________________________________
167void AliTPCCalibViewer::Delete(Option_t* option) {
168 //
169 // Should be called from AliTPCCalibViewerGUI class only.
170 // If you use Delete() do not call the destructor.
171 // All objects (except those contained in fListOfObjectsToBeDeleted) will be deleted, the file will be closed.
172 //
173
174 option = option; // to avoid warnings on compiling
175 if (fTree && fTreeMustBeDeleted) {
176 fTree->SetCacheSize(0);
177 fTree->Delete();
178 }
179 if (fFile)
180 delete fFile;
181 delete fListOfObjectsToBeDeleted;
182}
183
39bcd65d 184//_____________________________________________________________________________
185Int_t AliTPCCalibViewer::EasyDraw(const char* drawCommand, const char* sector, const char* cuts, const char* drawOptions, Bool_t writeDrawCommand) const {
186 //
187 // easy drawing of data, use '~' for abbreviation of '.fElements'
188 // example: EasyDraw("CETmean~-CETmean_mean", "A", "(CETmean~-CETmean_mean)>0")
189 // sector: sector-number - only the specified sector will be drwawn
190 // 'A'/'C' or 'a'/'c' - side A/C will be drawn
191 // 'ALL' - whole TPC will be drawn, projected on one side
192 // cuts: specifies cuts
193 // drawOptions: draw options like 'same'
194 // writeDrawCommand: write the command, that is passed to TTree::Draw
195 //
72d0ab7e 196
39bcd65d 197 TString drawStr(drawCommand);
198 TString sectorStr(sector);
199 sectorStr.ToUpper();
200 TString cutStr("");
72d0ab7e 201 //TString drawOptionsStr("profcolz ");
202 TString drawOptionsStr("");
39bcd65d 203 TRandom rnd(0);
204 Int_t rndNumber = rnd.Integer(10000);
72d0ab7e 205
206 if (drawOptions && strcmp(drawOptions, "") != 0)
39bcd65d 207 drawOptionsStr += drawOptions;
72d0ab7e 208 else
209 drawOptionsStr += "profcolz";
39bcd65d 210
211 if (sectorStr == "A") {
212 drawStr += ":gy.fElements:gx.fElements>>prof";
213 drawStr += rndNumber;
214 drawStr += "(330,-250,250,330,-250,250)";
215 cutStr += "(sector/18)%2==0 ";
216 }
217 else if (sectorStr == "C") {
218 drawStr += ":gy.fElements:gx.fElements>>prof";
219 drawStr += rndNumber;
220 drawStr += "(330,-250,250,330,-250,250)";
221 cutStr += "(sector/18)%2==1 ";
222 }
223 else if (sectorStr == "ALL") {
224 drawStr += ":gy.fElements:gx.fElements>>prof";
225 drawStr += rndNumber;
226 drawStr += "(330,-250,250,330,-250,250)";
227 }
228 else if (sectorStr.IsDigit()) {
229 Int_t isec = sectorStr.Atoi();
230 drawStr += ":rpad.fElements:row.fElements>>prof";
231 drawStr += rndNumber;
232 if (isec < 36 && isec >= 0)
233 drawStr += "(63,0,63,108,-54,54)";
234 else if (isec < 72 && isec >= 36)
235 drawStr += "(96,0,96,140,-70,70)";
236 else {
237 Error("EasyDraw","The TPC contains only sectors between 0 and 71.");
238 return -1;
239 }
240 cutStr += "(sector==";
241 cutStr += isec;
242 cutStr += ") ";
243 }
244
245 if (cuts && cuts[0] != 0) {
246 if (cutStr.Length() != 0) cutStr += "&& ";
247 cutStr += "(";
248 cutStr += cuts;
249 cutStr += ")";
250 }
251 drawStr.ReplaceAll("~", ".fElements");
252 cutStr.ReplaceAll("~", ".fElements");
253 if (writeDrawCommand) std::cout << "fTree->Draw(\"" << drawStr << "\", \"" << cutStr << "\", \"" << drawOptionsStr << "\");" << std::endl;
254 return fTree->Draw(drawStr.Data(), cutStr.Data(), drawOptionsStr.Data());
255}
256
72d0ab7e 257
39bcd65d 258Int_t AliTPCCalibViewer::EasyDraw(const char* drawCommand, Int_t sector, const char* cuts, const char* drawOptions, Bool_t writeDrawCommand) const {
259 //
260 // easy drawing of data, use '~' for abbreviation of '.fElements'
261 // example: EasyDraw("CETmean~-CETmean_mean", 34, "(CETmean~-CETmean_mean)>0")
262 // sector: sector-number - only the specified sector will be drwawn
263 // cuts: specifies cuts
264 // drawOptions: draw options like 'same'
265 // writeDrawCommand: write the command, that is passed to TTree::Draw
266 //
267 if (sector >= 0 && sector < 72) {
268 char sectorChr[3];
269 sprintf(sectorChr, "%i", sector);
270 return EasyDraw(drawCommand, sectorChr, cuts, drawOptions, writeDrawCommand);
271 }
272 Error("EasyDraw","The TPC contains only sectors between 0 and 71.");
273 return -1;
274}
275
72d0ab7e 276
39bcd65d 277//_____________________________________________________________________________
278Int_t AliTPCCalibViewer::EasyDraw1D(const char* drawCommand, const char* sector, const char* cuts, const char* drawOptions, Bool_t writeDrawCommand) const {
279 //
280 // easy drawing of data, use '~' for abbreviation of '.fElements'
281 // example: EasyDraw("CETmean~-CETmean_mean", "A", "(CETmean~-CETmean_mean)>0")
282 // sector: sector-number - the specified sector will be drwawn
283 // 'A'/'C' or 'a'/'c' - side A/C will be drawn
284 // 'ALL' - whole TPC will be drawn, projected on one side
285 // cuts: specifies cuts
286 // drawOptions: draw options like 'same'
287 // writeDrawCommand: write the command, that is passed to TTree::Draw
288 //
289
290 TString drawStr(drawCommand);
291 TString sectorStr(sector);
292 TString drawOptionsStr(drawOptions);
293 sectorStr.ToUpper();
294 TString cutStr("");
295
296 if (sectorStr == "A")
297 cutStr += "(sector/18)%2==0 ";
298 else if (sectorStr == "C")
299 cutStr += "(sector/18)%2==1 ";
300 else if (sectorStr.IsDigit()) {
301 Int_t isec = sectorStr.Atoi();
302 if (isec < 0 || isec > 71) {
303 Error("EasyDraw","The TPC contains only sectors between 0 and 71.");
304 return -1;
305 }
306 cutStr += "(sector==";
307 cutStr += isec;
308 cutStr += ") ";
309 }
310
311 if (cuts && cuts[0] != 0) {
312 if (cutStr.Length() != 0) cutStr += "&& ";
313 cutStr += "(";
314 cutStr += cuts;
315 cutStr += ")";
316 }
317
318 drawStr.ReplaceAll("~", ".fElements");
319 cutStr.ReplaceAll("~", ".fElements");
320 if (writeDrawCommand) std::cout << "fTree->Draw(\"" << drawStr << "\", \"" << cutStr << "\", \"" << drawOptionsStr << "\");" << std::endl;
321 return fTree->Draw(drawStr.Data(), cutStr.Data(), drawOptionsStr.Data());
322}
323
72d0ab7e 324
39bcd65d 325Int_t AliTPCCalibViewer::EasyDraw1D(const char* drawCommand, Int_t sector, const char* cuts, const char* drawOptions, Bool_t writeDrawCommand) const {
326 //
327 // easy drawing of data, use '~' for abbreviation of '.fElements'
328 // example: EasyDraw("CETmean~-CETmean_mean", 34, "(CETmean~-CETmean_mean)>0")
329 // sector: sector-number - the specified sector will be drwawn
330 // cuts: specifies cuts
331 // drawOptions: draw options like 'same'
332 // writeDrawCommand: write the command, that is passed to TTree::Draw
333 //
334
335 if (sector >= 0 && sector < 72) {
336 char sectorChr[3];
337 sprintf(sectorChr, "%i", sector);
338 return EasyDraw1D(drawCommand, sectorChr, cuts, drawOptions, writeDrawCommand);
339 }
340 Error("EasyDraw","The TPC contains only sectors between 0 and 71.");
341 return -1;
342}
343
39bcd65d 344
72d0ab7e 345Int_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 {
346 //
347 // Easy drawing of data, in principle the same as EasyDraw1D
348 // Difference: A line for the mean / median / LTM is drawn
349 // in 'sigmas' you can specify in which distance to the mean/median/LTM you want to see a line in sigma-units, separated by ';'
350 // 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.
351 // "plotMean", "plotMedian" and "plotLTM": what kind of lines do you want to see?
352 //
353 if (sector >= 0 && sector < 72) {
354 char sectorChr[3];
355 sprintf(sectorChr, "%i", sector);
356 return DrawHisto1D(drawCommand, sectorChr, cuts, sigmas, plotMean, plotMedian, plotLTM);
357 }
358 Error("DrawHisto1D","The TPC contains only sectors between 0 and 71.");
359 return -1;
360}
39bcd65d 361
72d0ab7e 362
363Int_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 {
364 //
365 // Easy drawing of data, in principle the same as EasyDraw1D
366 // Difference: A line for the mean / median / LTM is drawn
367 // in 'sigmas' you can specify in which distance to the mean/median/LTM you want to see a line in sigma-units, separated by ';'
368 // 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.
369 // "plotMean", "plotMedian" and "plotLTM": what kind of lines do you want to see?
370 //
39bcd65d 371 Int_t oldOptStat = gStyle->GetOptStat();
372 gStyle->SetOptStat(0000000);
72d0ab7e 373 Double_t ltmFraction = 0.8;
374
375 TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");
376 TVectorF nsigma(sigmasTokens->GetEntriesFast());
377 for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
378 TString str(((TObjString*)sigmasTokens->At(i))->GetString());
379 Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
380 nsigma[i] = sig;
39bcd65d 381 }
382
72d0ab7e 383 TString drawStr(drawCommand);
384 drawStr += " >> tempHist";
385 Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts);
386 TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
387 // FIXME is this histogram deleted automatically?
388 Double_t *values = fTree->GetV1(); // value is the array containing 'entries' numbers
389
390 Double_t mean = TMath::Mean(entries, values);
391 Double_t median = TMath::Median(entries, values);
392 Double_t sigma = TMath::RMS(entries, values);
393 Double_t maxY = htemp->GetMaximum();
394
39bcd65d 395 char c[500];
72d0ab7e 396 TLegend * legend = new TLegend(.7,.7, .99, .99, "Statistical information");
397// sprintf(c, "%s, sector: %i", type, sector);
39bcd65d 398 fListOfObjectsToBeDeleted->Add(legend);
399
39bcd65d 400 if (plotMean) {
72d0ab7e 401 // draw Mean
402 TLine* line = new TLine(mean, 0, mean, maxY);
39bcd65d 403 fListOfObjectsToBeDeleted->Add(line);
404 line->SetLineColor(kRed);
405 line->SetLineWidth(2);
406 line->SetLineStyle(1);
407 line->Draw();
72d0ab7e 408 sprintf(c, "Mean: %f", mean);
39bcd65d 409 legend->AddEntry(line, c, "l");
72d0ab7e 410 // draw sigma lines
39bcd65d 411 for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
72d0ab7e 412 TLine* linePlusSigma = new TLine(mean + nsigma[i] * sigma, 0, mean + nsigma[i] * sigma, maxY);
39bcd65d 413 fListOfObjectsToBeDeleted->Add(linePlusSigma);
414 linePlusSigma->SetLineColor(kRed);
72d0ab7e 415 linePlusSigma->SetLineStyle(2 + i);
39bcd65d 416 linePlusSigma->Draw();
72d0ab7e 417 TLine* lineMinusSigma = new TLine(mean - nsigma[i] * sigma, 0, mean - nsigma[i] * sigma, maxY);
39bcd65d 418 fListOfObjectsToBeDeleted->Add(lineMinusSigma);
419 lineMinusSigma->SetLineColor(kRed);
72d0ab7e 420 lineMinusSigma->SetLineStyle(2 + i);
39bcd65d 421 lineMinusSigma->Draw();
72d0ab7e 422 sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
39bcd65d 423 legend->AddEntry(lineMinusSigma, c, "l");
424 }
425 }
39bcd65d 426 if (plotMedian) {
72d0ab7e 427 // draw median
428 TLine* line = new TLine(median, 0, median, maxY);
39bcd65d 429 fListOfObjectsToBeDeleted->Add(line);
430 line->SetLineColor(kBlue);
431 line->SetLineWidth(2);
432 line->SetLineStyle(1);
433 line->Draw();
72d0ab7e 434 sprintf(c, "Median: %f", median);
39bcd65d 435 legend->AddEntry(line, c, "l");
72d0ab7e 436 // draw sigma lines
39bcd65d 437 for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
72d0ab7e 438 TLine* linePlusSigma = new TLine(median + nsigma[i] * sigma, 0, median + nsigma[i]*sigma, maxY);
39bcd65d 439 fListOfObjectsToBeDeleted->Add(linePlusSigma);
440 linePlusSigma->SetLineColor(kBlue);
72d0ab7e 441 linePlusSigma->SetLineStyle(2 + i);
39bcd65d 442 linePlusSigma->Draw();
72d0ab7e 443 TLine* lineMinusSigma = new TLine(median - nsigma[i] * sigma, 0, median - nsigma[i]*sigma, maxY);
39bcd65d 444 fListOfObjectsToBeDeleted->Add(lineMinusSigma);
445 lineMinusSigma->SetLineColor(kBlue);
72d0ab7e 446 lineMinusSigma->SetLineStyle(2 + i);
39bcd65d 447 lineMinusSigma->Draw();
72d0ab7e 448 sprintf(c, "%i #sigma = %f",(Int_t)(nsigma[i]), (Float_t)(nsigma[i] * sigma));
39bcd65d 449 legend->AddEntry(lineMinusSigma, c, "l");
450 }
451 }
39bcd65d 452 if (plotLTM) {
72d0ab7e 453 // draw LTM
454 Double_t ltmRms = 0;
455 Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
456 TLine* line = new TLine(ltm, 0, ltm, maxY);
39bcd65d 457 fListOfObjectsToBeDeleted->Add(line);
458 line->SetLineColor(kGreen+2);
459 line->SetLineWidth(2);
460 line->SetLineStyle(1);
461 line->Draw();
72d0ab7e 462 sprintf(c, "LTM: %f", ltm);
39bcd65d 463 legend->AddEntry(line, c, "l");
72d0ab7e 464 // draw sigma lines
39bcd65d 465 for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
72d0ab7e 466 TLine* linePlusSigma = new TLine(ltm + nsigma[i] * ltmRms, 0, ltm + nsigma[i] * ltmRms, maxY);
39bcd65d 467 fListOfObjectsToBeDeleted->Add(linePlusSigma);
468 linePlusSigma->SetLineColor(kGreen+2);
469 linePlusSigma->SetLineStyle(2+i);
470 linePlusSigma->Draw();
471
72d0ab7e 472 TLine* lineMinusSigma = new TLine(ltm - nsigma[i] * ltmRms, 0, ltm - nsigma[i] * ltmRms, maxY);
39bcd65d 473 fListOfObjectsToBeDeleted->Add(lineMinusSigma);
474 lineMinusSigma->SetLineColor(kGreen+2);
475 lineMinusSigma->SetLineStyle(2+i);
476 lineMinusSigma->Draw();
72d0ab7e 477 sprintf(c, "%i #sigma = %f", (Int_t)(nsigma[i]), (Float_t)(nsigma[i] * ltmRms));
39bcd65d 478 legend->AddEntry(lineMinusSigma, c, "l");
479 }
480 }
72d0ab7e 481 if (!plotMean && !plotMedian && !plotLTM) return -1;
39bcd65d 482 legend->Draw();
483 gStyle->SetOptStat(oldOptStat);
72d0ab7e 484 return 1;
39bcd65d 485}
486
39bcd65d 487
72d0ab7e 488Int_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 {
489 //
490 // 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
491 // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'array', 'n' specifies the length of the array
492 // '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
493 // 'nbins': number of bins, 'binLow': first bin, 'binUp': last bin
494 // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
495 // sigmaStep: the binsize of the generated histogram
496 // Begin_Latex
497 // 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 }
498 // End_Latex
499 //
500 //
501 // Creates a histogram, where you can see, how much of the data are inside sigma-intervals
502 // around the mean/median/LTM
503 // with drawCommand, sector and cuts you specify your input data, see EasyDraw
504 // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma)
505 // sigmaStep: the binsize of the generated histogram
506 // plotMean/plotMedian/plotLTM: specifies where to put the center
507 //
508 if (sector >= 0 && sector < 72) {
509 char sectorChr[3];
510 sprintf(sectorChr, "%i", sector);
511 return SigmaCut(drawCommand, sectorChr, cuts, sigmaMax, plotMean, plotMedian, plotLTM, pm, sigmas, sigmaStep);
512 }
513 Error("SigmaCut","The TPC contains only sectors between 0 and 71.");
514 return -1;
515}
516
517
518Int_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 {
519 //
520 // Creates a histogram, where you can see, how much of the data are inside sigma-intervals
521 // around the mean/median/LTM
522 // with drawCommand, sector and cuts you specify your input data, see EasyDraw
523 // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma)
524 // sigmaStep: the binsize of the generated histogram
525 // plotMean/plotMedian/plotLTM: specifies where to put the center
526 //
527
528 Double_t ltmFraction = 0.8;
529
530 TString drawStr(drawCommand);
531 drawStr += " >> tempHist";
532
533 Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
534 TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
535 // FIXME is this histogram deleted automatically?
536 Double_t *values = fTree->GetV1(); // value is the array containing 'entries' numbers
537
538 Double_t mean = TMath::Mean(entries, values);
539 Double_t median = TMath::Median(entries, values);
540 Double_t sigma = TMath::RMS(entries, values);
541
542 TLegend * legend = new TLegend(.7,.7, .99, .99, "Cumulative");
543 fListOfObjectsToBeDeleted->Add(legend);
544 TH1F *cutHistoMean = 0;
545 TH1F *cutHistoMedian = 0;
546 TH1F *cutHistoLTM = 0;
547
548 TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");
549 TVectorF nsigma(sigmasTokens->GetEntriesFast());
550 for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
551 TString str(((TObjString*)sigmasTokens->At(i))->GetString());
552 Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
553 nsigma[i] = sig;
554 }
555
556 if (plotMean) {
557 cutHistoMean = AliTPCCalibViewer::SigmaCut(htemp, mean, sigma, sigmaMax, sigmaStep, pm);
558 if (cutHistoMean) {
559 fListOfObjectsToBeDeleted->Add(cutHistoMean);
560 cutHistoMean->SetLineColor(kRed);
561 legend->AddEntry(cutHistoMean, "Mean", "l");
562 cutHistoMean->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
563 cutHistoMean->Draw();
564 DrawLines(cutHistoMean, nsigma, legend, kRed, pm);
565 } // if (cutHistoMean)
566
567 }
568 if (plotMedian) {
569 cutHistoMedian = AliTPCCalibViewer::SigmaCut(htemp, median, sigma, sigmaMax, sigmaStep, pm);
570 if (cutHistoMedian) {
571 fListOfObjectsToBeDeleted->Add(cutHistoMedian);
572 cutHistoMedian->SetLineColor(kBlue);
573 legend->AddEntry(cutHistoMedian, "Median", "l");
574 cutHistoMedian->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
575 if (plotMean && cutHistoMean) cutHistoMedian->Draw("same");
576 else cutHistoMedian->Draw();
577 DrawLines(cutHistoMedian, nsigma, legend, kBlue, pm);
578 } // if (cutHistoMedian)
579 }
580 if (plotLTM) {
581 Double_t ltmRms = 0;
582 Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
583 cutHistoLTM = AliTPCCalibViewer::SigmaCut(htemp, ltm, ltmRms, sigmaMax, sigmaStep, pm);
584 if (cutHistoLTM) {
585 fListOfObjectsToBeDeleted->Add(cutHistoLTM);
586 cutHistoLTM->SetLineColor(kGreen+2);
587 legend->AddEntry(cutHistoLTM, "LTM", "l");
588 cutHistoLTM->SetTitle(Form("%s, cumulative; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
589 if (plotMean && cutHistoMean || plotMedian && cutHistoMedian) cutHistoLTM->Draw("same");
590 else cutHistoLTM->Draw();
591 DrawLines(cutHistoLTM, nsigma, legend, kGreen+2, pm);
592 }
593 }
594 if (!plotMean && !plotMedian && !plotLTM) return -1;
595 legend->Draw();
596 return 1;
597}
598
599
600
601
602Int_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 {
603 //
604 // 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"
605 // "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
606 // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
607 // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
608 // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
609 // The actual work is done on the array.
610 /* Begin_Latex
611 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 }
612 End_Latex
613 */
614 if (sector >= 0 && sector < 72) {
615 char sectorChr[3];
616 sprintf(sectorChr, "%i", sector);
617 return Integrate(drawCommand, sectorChr, cuts, sigmaMax, plotMean, plotMedian, plotLTM, sigmas, sigmaStep);
618 }
619 Error("Integrate","The TPC contains only sectors between 0 and 71.");
620 return -1;
621
622}
623
624
625Int_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 {
626 //
627 // 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"
628 // "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
629 // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
630 // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
631 // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
632 // The actual work is done on the array.
633 /* Begin_Latex
634 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 }
635 End_Latex
636 */
637
638 Double_t ltmFraction = 0.8;
639
640 TString drawStr(drawCommand);
641 drawStr += " >> tempHist";
642
643 Int_t entries = EasyDraw1D(drawStr.Data(), sector, cuts, "goff");
644 TH1F *htemp = (TH1F*)gDirectory->Get("tempHist");
645 // FIXME is this histogram deleted automatically?
646 Double_t *values = fTree->GetV1(); // value is the array containing 'entries' numbers
647
648 Double_t mean = TMath::Mean(entries, values);
649 Double_t median = TMath::Median(entries, values);
650 Double_t sigma = TMath::RMS(entries, values);
651
652 TObjArray *sigmasTokens = TString(sigmas).Tokenize(";");
653 TVectorF nsigma(sigmasTokens->GetEntriesFast());
654 for (Int_t i = 0; i < sigmasTokens->GetEntriesFast(); i++) {
655 TString str(((TObjString*)sigmasTokens->At(i))->GetString());
656 Double_t sig = (str.IsFloat()) ? str.Atof() : 0;
657 nsigma[i] = sig;
658 }
659
660 TLegend * legend = new TLegend(.7,.7, .99, .99, "Integrated histogram");
661 fListOfObjectsToBeDeleted->Add(legend);
662 TH1F *integralHistoMean = 0;
663 TH1F *integralHistoMedian = 0;
664 TH1F *integralHistoLTM = 0;
665
666 if (plotMean) {
667 integralHistoMean = AliTPCCalibViewer::Integrate(htemp, mean, sigma, sigmaMax, sigmaStep);
668 if (integralHistoMean) {
669 fListOfObjectsToBeDeleted->Add(integralHistoMean);
670 integralHistoMean->SetLineColor(kRed);
671 legend->AddEntry(integralHistoMean, "Mean", "l");
672 integralHistoMean->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
673 integralHistoMean->Draw();
674 DrawLines(integralHistoMean, nsigma, legend, kRed, kTRUE);
39bcd65d 675 }
72d0ab7e 676 }
677 if (plotMedian) {
678 integralHistoMedian = AliTPCCalibViewer::Integrate(htemp, median, sigma, sigmaMax, sigmaStep);
679 if (integralHistoMedian) {
680 fListOfObjectsToBeDeleted->Add(integralHistoMedian);
681 integralHistoMedian->SetLineColor(kBlue);
682 legend->AddEntry(integralHistoMedian, "Median", "l");
683 integralHistoMedian->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
684 if (plotMean && integralHistoMean) integralHistoMedian->Draw("same");
685 else integralHistoMedian->Draw();
686 DrawLines(integralHistoMedian, nsigma, legend, kBlue, kTRUE);
39bcd65d 687 }
72d0ab7e 688 }
689 if (plotLTM) {
690 Double_t ltmRms = 0;
691 Double_t ltm = GetLTM(entries, values, &ltmRms, ltmFraction);
692 integralHistoLTM = AliTPCCalibViewer::Integrate(htemp, ltm, ltmRms, sigmaMax, sigmaStep);
693 if (integralHistoLTM) {
694 fListOfObjectsToBeDeleted->Add(integralHistoLTM);
695 integralHistoLTM->SetLineColor(kGreen+2);
696 legend->AddEntry(integralHistoLTM, "LTM", "l");
697 integralHistoLTM->SetTitle(Form("%s, integrated; Multiples of #sigma; Fraction of included data", htemp->GetTitle()));
698 if (plotMean && integralHistoMean || plotMedian && integralHistoMedian) integralHistoLTM->Draw("same");
699 else integralHistoLTM->Draw();
700 DrawLines(integralHistoLTM, nsigma, legend, kGreen+2, kTRUE);
39bcd65d 701 }
702 }
72d0ab7e 703 if (!plotMean && !plotMedian && !plotLTM) return -1;
704 legend->Draw();
705 return 1;
706}
707
708
709void AliTPCCalibViewer::DrawLines(TH1F *histogram, TVectorF nsigma, TLegend *legend, Int_t color, Bool_t pm) const {
710 //
711 // Private function for SigmaCut(...) and Integrate(...)
712 // Draws lines into the given histogram, specified by "nsigma", the lines are addeed to the legend
713 //
39bcd65d 714
72d0ab7e 715 // start to draw the lines, loop over requested sigmas
39bcd65d 716 char c[500];
72d0ab7e 717 for (Int_t i = 0; i < nsigma.GetNoElements(); i++) {
718 if (!pm) {
719 Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
720 TLine* lineUp = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
721 fListOfObjectsToBeDeleted->Add(lineUp);
722 lineUp->SetLineColor(color);
723 lineUp->SetLineStyle(2 + i);
724 lineUp->Draw();
725 TLine* lineLeft = new TLine(nsigma[i], histogram->GetBinContent(bin), 0, histogram->GetBinContent(bin));
726 fListOfObjectsToBeDeleted->Add(lineLeft);
727 lineLeft->SetLineColor(color);
728 lineLeft->SetLineStyle(2 + i);
729 lineLeft->Draw();
730 sprintf(c, "Fraction(%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
731 legend->AddEntry(lineLeft, c, "l");
732 }
733 else { // if (pm)
734 Int_t bin = histogram->GetXaxis()->FindBin(nsigma[i]);
735 TLine* lineUp1 = new TLine(nsigma[i], 0, nsigma[i], histogram->GetBinContent(bin));
736 fListOfObjectsToBeDeleted->Add(lineUp1);
737 lineUp1->SetLineColor(color);
738 lineUp1->SetLineStyle(2 + i);
739 lineUp1->Draw();
740 TLine* lineLeft1 = new TLine(nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
741 fListOfObjectsToBeDeleted->Add(lineLeft1);
742 lineLeft1->SetLineColor(color);
743 lineLeft1->SetLineStyle(2 + i);
744 lineLeft1->Draw();
745 sprintf(c, "Fraction(+%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
746 legend->AddEntry(lineLeft1, c, "l");
747 bin = histogram->GetXaxis()->FindBin(-nsigma[i]);
748 TLine* lineUp2 = new TLine(-nsigma[i], 0, -nsigma[i], histogram->GetBinContent(bin));
749 fListOfObjectsToBeDeleted->Add(lineUp2);
750 lineUp2->SetLineColor(color);
751 lineUp2->SetLineStyle(2 + i);
752 lineUp2->Draw();
753 TLine* lineLeft2 = new TLine(-nsigma[i], histogram->GetBinContent(bin), histogram->GetBinLowEdge(0)+histogram->GetBinWidth(0), histogram->GetBinContent(bin));
754 fListOfObjectsToBeDeleted->Add(lineLeft2);
755 lineLeft2->SetLineColor(color);
756 lineLeft2->SetLineStyle(2 + i);
757 lineLeft2->Draw();
758 sprintf(c, "Fraction(-%f #sigma) = %f",nsigma[i], histogram->GetBinContent(bin));
759 legend->AddEntry(lineLeft2, c, "l");
760 }
761 } // for (Int_t i = 0; i < nsigma.GetNoElements(); i++)
762}
763
764
765
766
767
768
769/////////////////
770// Array tools //
771/////////////////
772
773
774Int_t AliTPCCalibViewer::GetBin(Float_t value, Int_t nbins, Double_t binLow, Double_t binUp){
775 // Returns the 'bin' for 'value'
776 // The interval between 'binLow' and 'binUp' is divided into 'nbins' equidistant bins
777 // avoid index out of bounds error: 'if (bin < binLow) bin = binLow' and vice versa
778 /* Begin_Latex
779 GetBin(value) = #frac{nbins - 1}{binUp - binLow} #upoint (value - binLow) +1
780 End_Latex
781 */
782
783 Int_t bin = TMath::Nint( (Float_t)(value - binLow) / (Float_t)(binUp - binLow) * (nbins-1) ) + 1;
784 // avoid index out of bounds:
785 if (value < binLow) bin = 0;
786 if (value > binUp) bin = nbins + 1;
787 return bin;
39bcd65d 788
72d0ab7e 789}
39bcd65d 790
72d0ab7e 791
792Double_t AliTPCCalibViewer::GetLTM(Int_t n, Double_t *array, Double_t *sigma, Double_t fraction){
793 //
794 // returns the LTM and sigma
795 //
796 Double_t *ddata = new Double_t[n];
797 Double_t mean = 0, lsigma = 0;
798 UInt_t nPoints = 0;
799 for (UInt_t i = 0; i < (UInt_t)n; i++) {
800 ddata[nPoints]= array[nPoints];
801 nPoints++;
802 }
803 Int_t hh = TMath::Min(TMath::Nint(fraction * nPoints), Int_t(n));
804 AliMathBase::EvaluateUni(nPoints, ddata, mean, lsigma, hh);
805 if (sigma) *sigma = lsigma;
806 delete [] ddata;
807 return mean;
39bcd65d 808}
809
810
72d0ab7e 811TH1F* AliTPCCalibViewer::SigmaCut(TH1F *histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep, Bool_t pm) {
812 //
813 // 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
814 // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'histogram'
815 // '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
816 // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
817 // sigmaStep: the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
818 // pm: Decide weather Begin_Latex t > 0 End_Latex (first case) or Begin_Latex t End_Latex arbitrary (secound case)
819 // The actual work is done on the array.
820 /* Begin_Latex
821 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
822 or
823 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 }
824 End_Latex
825 begin_macro(source)
826 {
827 Float_t mean = 0;
828 Float_t sigma = 1.5;
829 Float_t sigmaMax = 4;
830 gROOT->SetStyle("Plain");
831 TH1F *distribution = new TH1F("Distribution", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
832 TRandom rand(23);
833 for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
834 Float_t *ar = distribution->GetArray();
835
836 TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas", "", 350, 350);
837 macro_example_canvas->Divide(0,3);
838 TVirtualPad *pad1 = macro_example_canvas->cd(1);
839 pad1->SetGridy();
840 pad1->SetGridx();
841 distribution->Draw();
842 TVirtualPad *pad2 = macro_example_canvas->cd(2);
843 pad2->SetGridy();
844 pad2->SetGridx();
845
846 TH1F *shist = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax);
847 shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
848 shist->Draw();
849 TVirtualPad *pad3 = macro_example_canvas->cd(3);
850 pad3->SetGridy();
851 pad3->SetGridx();
852 TH1F *shistPM = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax, -1, kTRUE);
853 shistPM->Draw();
854 return macro_example_canvas;
855 }
856 end_macro
857 */
858
859 Float_t *array = histogram->GetArray();
860 Int_t nbins = histogram->GetXaxis()->GetNbins();
861 Float_t binLow = histogram->GetXaxis()->GetXmin();
862 Float_t binUp = histogram->GetXaxis()->GetXmax();
863 return AliTPCCalibViewer::SigmaCut(nbins, array, mean, sigma, nbins, binLow, binUp, sigmaMax, sigmaStep, pm);
864}
865
866
867
868TH1F* 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){
869 //
870 // 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
871 // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'array', 'n' specifies the length of the array
872 // '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
873 // 'nbins': number of bins, 'binLow': first bin, 'binUp': last bin
874 // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
875 // sigmaStep: the binsize of the generated histogram
876 // Here the actual work is done.
877
878 if (sigma == 0) return 0;
879 Float_t binWidth = (binUp-binLow)/(nbins - 1);
880 if (sigmaStep <= 0) sigmaStep = binWidth;
881 Int_t kbins = (Int_t)(sigmaMax * sigma / sigmaStep) + 1; // + 1 due to overflow bin in histograms
882 if (pm) kbins = 2 * (Int_t)(sigmaMax * sigma / sigmaStep) + 1;
883 Float_t kbinLow = !pm ? 0 : -sigmaMax;
884 Float_t kbinUp = sigmaMax;
885 TH1F *hist = new TH1F("sigmaCutHisto","Cumulative; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp);
886 hist->SetDirectory(0);
887 hist->Reset();
888
889 // calculate normalization
890 Double_t normalization = 0;
891 for (Int_t i = 0; i <= n; i++) {
892 normalization += array[i];
893 }
894
895 // given units: units from given histogram
896 // sigma units: in units of sigma
897 // iDelta: integrate in interval (mean +- iDelta), given units
898 // x: ofset from mean for integration, given units
899 // hist: needs
900
901// printf("nbins: %i, binLow: %f, binUp: %f \n", nbins, binLow, binUp);
902 // fill histogram
903 for (Float_t iDelta = 0; iDelta <= sigmaMax * sigma; iDelta += sigmaStep) {
904 // integrate array
905 Double_t valueP = array[GetBin(mean, nbins, binLow, binUp)];
906 Double_t valueM = array[GetBin(mean-binWidth, nbins, binLow, binUp)];
907 // add bin of mean value only once to the histogram
908// printf("++ adding bins: ");
909 for (Float_t x = binWidth; x <= iDelta; x += binWidth) {
910 valueP += (mean + x <= binUp) ? array[GetBin(mean + x, nbins, binLow, binUp)] : 0;
911 valueM += (mean-binWidth - x >= binLow) ? array[GetBin(mean-binWidth - x, nbins, binLow, binUp)] : 0;
912// printf("%i, ", GetBin(mean + x, nbins, binLow, binUp));
913 }
914// printf("\n");
915 if (valueP / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail +++ \n", valueP, normalization);
916 if (valueP / normalization > 100) return hist;
917 if (valueM / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail +++ \n", valueM, normalization);
918 if (valueM / normalization > 100) return hist;
919 valueP = (valueP / normalization);
920 valueM = (valueM / normalization);
921 if (pm) {
922 Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
923 hist->SetBinContent(bin, valueP);
924 bin = GetBin(-iDelta/sigma, kbins, kbinLow, kbinUp);
925 hist->SetBinContent(bin, valueM);
926 }
927 else { // if (!pm)
928 Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
929 hist->SetBinContent(bin, valueP + valueM);
930// printf(" first integration bin: %i, last integration bin in + direction: %i \n", GetBin(mean+binWidth, nbins, binLow, binUp), GetBin(iDelta, nbins, binLow, binUp));
931// printf(" first integration bin: %i, last integration bin in - direction: %i \n", GetBin(mean+binWidth, nbins, binLow, binUp), GetBin(-iDelta, nbins, binLow, binUp));
932// printf(" value: %f, normalization: %f, iDelta: %f, Bin: %i \n", valueP+valueM, normalization, iDelta, bin);
933 }
934 }
935 //hist->SetMaximum(0.7);
936 if (!pm) hist->SetMaximum(1.2);
937 return hist;
938}
939
940
941TH1F* AliTPCCalibViewer::SigmaCut(Int_t n, Double_t *array, Double_t mean, Double_t sigma, Int_t nbins, Double_t *xbins, Double_t sigmaMax){
942 //
943 // SigmaCut for variable binsize
944 // NOT YET IMPLEMENTED !!!
945 //
946 printf("SigmaCut with variable binsize, Not yet implemented\n");
947 // avoid compiler warnings:
948 n=n;
949 mean=mean;
950 sigma=sigma;
951 nbins=nbins;
952 sigmaMax=sigmaMax;
953 array=array;
954 xbins=xbins;
955
956 return 0;
957}
958
959
960TH1F* AliTPCCalibViewer::Integrate(TH1F *histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
961 //
962 // 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"
963 // "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
964 // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
965 // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
966 // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
967 // The actual work is done on the array.
968 /* Begin_Latex
969 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 }
970 End_Latex
971 begin_macro(source)
972 {
973 Float_t mean = 0;
974 Float_t sigma = 1.5;
975 Float_t sigmaMax = 4;
976 gROOT->SetStyle("Plain");
977 TH1F *distribution = new TH1F("Distribution", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
978 TRandom rand(23);
979 for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
980 Float_t *ar = distribution->GetArray();
981
982 TCanvas* macro_example_canvas = new TCanvas("macro_example_canvas", "", 350, 350);
983 macro_example_canvas->Divide(0,2);
984 TVirtualPad *pad1 = macro_example_canvas->cd(1);
985 pad1->SetGridy();
986 pad1->SetGridx();
987 distribution->Draw();
988 TVirtualPad *pad2 = macro_example_canvas->cd(2);
989 pad2->SetGridy();
990 pad2->SetGridx();
991 TH1F *shist = AliTPCCalibViewer::Integrate(distribution, mean, sigma, sigmaMax);
992 shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
993 shist->Draw();
994
995 return macro_example_canvas;
996 }
997 end_macro
998 */
999
1000
1001 Float_t *array = histogram->GetArray();
1002 Int_t nbins = histogram->GetXaxis()->GetNbins();
1003 Float_t binLow = histogram->GetXaxis()->GetXmin();
1004 Float_t binUp = histogram->GetXaxis()->GetXmax();
1005 return AliTPCCalibViewer::Integrate(nbins, array, nbins, binLow, binUp, mean, sigma, sigmaMax, sigmaStep);
1006}
1007
1008
1009TH1F* 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){
1010 // 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"
1011 // "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
1012 // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
1013 // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
1014 // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
1015 // Here the actual work is done.
1016
1017 Bool_t givenUnits = kTRUE;
1018 if (sigma != 0 && sigmaMax != 0) givenUnits = kFALSE;
1019 if (givenUnits) {
1020 sigma = 1;
1021 sigmaMax = (binUp - binLow) / 2.;
1022 }
1023
1024 Float_t binWidth = (binUp-binLow)/(nbins - 1);
1025 if (sigmaStep <= 0) sigmaStep = binWidth;
1026 Int_t kbins = (Int_t)(sigmaMax * sigma / sigmaStep) + 1; // + 1 due to overflow bin in histograms
1027 Float_t kbinLow = givenUnits ? binLow : -sigmaMax;
1028 Float_t kbinUp = givenUnits ? binUp : sigmaMax;
1029 TH1F *hist = 0;
1030 if (givenUnits) hist = new TH1F("integratedHisto","Integrated Histogram; Given x; Fraction of included data", kbins, kbinLow, kbinUp);
1031 if (!givenUnits) hist = new TH1F("integratedHisto","Integrated Histogram; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp);
1032 hist->SetDirectory(0);
1033 hist->Reset();
1034
1035 // calculate normalization
1036 // printf("calculating normalization, integrating from bin 1 to %i \n", n);
1037 Double_t normalization = 0;
1038 for (Int_t i = 1; i <= n; i++) {
1039 normalization += array[i];
1040 }
1041 // printf("normalization: %f \n", normalization);
1042
1043 // given units: units from given histogram
1044 // sigma units: in units of sigma
1045 // iDelta: integrate in interval (mean +- iDelta), given units
1046 // x: ofset from mean for integration, given units
1047 // hist: needs
1048
1049 // fill histogram
1050 for (Float_t iDelta = mean - sigmaMax * sigma; iDelta <= mean + sigmaMax * sigma; iDelta += sigmaStep) {
1051 // integrate array
1052 Double_t value = 0;
1053 for (Float_t x = mean - sigmaMax * sigma; x <= iDelta; x += binWidth) {
1054 value += (x <= binUp && x >= binLow) ? array[GetBin(x, nbins, binLow, binUp)] : 0;
1055 }
1056 if (value / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail +++ \n", value, normalization);
1057 if (value / normalization > 100) return hist;
1058 Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
1059 // printf("first integration bin: %i, last integration bin: %i \n", GetBin(mean - sigmaMax * sigma, nbins, binLow, binUp), GetBin(iDelta, nbins, binLow, binUp));
1060 // printf("value: %f, normalization: %f, normalized value: %f, iDelta: %f, Bin: %i \n", value, normalization, value/normalization, iDelta, bin);
1061 value = (value / normalization);
1062 hist->SetBinContent(bin, value);
1063 }
1064 return hist;
1065}
1066
1067
1068
1069
1070
1071////////////////////////
1072// end of Array tools //
1073////////////////////////
1074
1075
1076
39bcd65d 1077//_____________________________________________________________________________
1078AliTPCCalPad* AliTPCCalibViewer::GetCalPad(const char* desiredData, char* cuts, char* calPadName) const {
1079 //
1080 // creates a AliTPCCalPad out of the 'desiredData'
1081 // the functionality of EasyDraw1D is used
1082 // calPadName specifies the name of the created AliTPCCalPad
1083 // - this takes a while -
1084 //
1085 TString drawStr(desiredData);
1086 drawStr.Append(":channel~");
1087 AliTPCCalPad * createdCalPad = new AliTPCCalPad(calPadName, calPadName);
1088 Int_t entries = 0;
1089 for (Int_t sec = 0; sec < 72; sec++) {
1090 entries = EasyDraw1D(drawStr.Data(), (Int_t)sec, cuts, "goff");
a6d2bd0c 1091 if (entries == -1) return 0;
39bcd65d 1092 for (Int_t i = 0; i < entries; i++)
1093 createdCalPad->GetCalROC(sec)->SetValue((UInt_t)(fTree->GetV2()[i]), (Float_t)(fTree->GetV1()[i]));
1094 }
1095 return createdCalPad;
1096}
1097
1098//_____________________________________________________________________________
1099AliTPCCalROC* AliTPCCalibViewer::GetCalROC(const char* desiredData, UInt_t sector, char* cuts) const {
1100 //
1101 // creates a AliTPCCalROC out of the desiredData
1102 // the functionality of EasyDraw1D is used
1103 // sector specifies the sector of the created AliTPCCalROC
1104 //
1105 TString drawStr(desiredData);
1106 drawStr.Append(":channel~");
1107 Int_t entries = EasyDraw1D(drawStr.Data(), (Int_t)sector, cuts, "goff");
a6d2bd0c 1108 if (entries == -1) return 0;
39bcd65d 1109 AliTPCCalROC * createdROC = new AliTPCCalROC(sector);
1110 for (Int_t i = 0; i < entries; i++)
1111 createdROC->SetValue((UInt_t)(fTree->GetV2()[i]), fTree->GetV1()[i]);
1112 return createdROC;
1113}
1114
1115
1116TObjArray* AliTPCCalibViewer::GetListOfVariables(Bool_t printList) {
1117 //
1118 // scan the tree - produces a list of available variables in the tree
1119 // printList: print the list to the screen, after the scan is done
1120 //
1121 TObjArray* arr = new TObjArray();
1122 TObjString* str = 0;
1123 Int_t nentries = fTree->GetListOfBranches()->GetEntries();
1124 for (Int_t i = 0; i < nentries; i++) {
1125 str = new TObjString(fTree->GetListOfBranches()->At(i)->GetName());
1126 str->String().ReplaceAll("_Median", "");
1127 str->String().ReplaceAll("_Mean", "");
1128 str->String().ReplaceAll("_RMS", "");
1129 str->String().ReplaceAll("_LTM", "");
1130 str->String().ReplaceAll("_OutlierCutted", "");
1131 str->String().ReplaceAll(".", "");
1132 if (!arr->FindObject(str) &&
1133 !(str->String() == "channel" || str->String() == "gx" || str->String() == "gy" ||
1134 str->String() == "lx" || str->String() == "ly" || str->String() == "pad" ||
1135 str->String() == "row" || str->String() == "rpad" || str->String() == "sector" ))
1136 arr->Add(str);
1137 }
1138 arr->Sort();
1139
1140 if (printList) {
1141 TIterator* iter = arr->MakeIterator();
1142 iter->Reset();
1143 TObjString* currentStr = 0;
1144 while ( (currentStr = (TObjString*)(iter->Next())) ) {
1145 std::cout << currentStr->GetString().Data() << std::endl;
1146 }
1147 delete iter;
1148 }
1149 return arr;
1150}
1151
1152
72d0ab7e 1153TObjArray* AliTPCCalibViewer::GetListOfNormalizationVariables(Bool_t printList) const{
39bcd65d 1154 //
1155 // produces a list of available variables for normalization in the tree
1156 // printList: print the list to the screen, after the scan is done
1157 //
1158 TObjArray* arr = new TObjArray();
1159 arr->Add(new TObjString("_Mean"));
1160 arr->Add(new TObjString("_Mean_OutlierCutted"));
1161 arr->Add(new TObjString("_Median"));
1162 arr->Add(new TObjString("_Median_OutlierCutted"));
1163 arr->Add(new TObjString("_LTM"));
1164 arr->Add(new TObjString("_LTM_OutlierCutted"));
1165 arr->Add(new TObjString("LFitIntern_4_8.fElements"));
1166 arr->Add(new TObjString("GFitIntern_Lin.fElements"));
1167 arr->Add(new TObjString("GFitIntern_Par.fElements"));
a6d2bd0c 1168 arr->Add(new TObjString("FitLinLocal"));
1169 arr->Add(new TObjString("FitLinGlobal"));
72d0ab7e 1170 arr->Add(new TObjString("FitParLocal"));
1171 arr->Add(new TObjString("FitParGlobal"));
39bcd65d 1172
1173 if (printList) {
1174 TIterator* iter = arr->MakeIterator();
1175 iter->Reset();
1176 TObjString* currentStr = 0;
1177 while ((currentStr = (TObjString*)(iter->Next()))) {
1178 std::cout << currentStr->GetString().Data() << std::endl;
1179 }
1180 delete iter;
1181 }
1182 return arr;
1183}
1184
1185
1186TFriendElement* AliTPCCalibViewer::AddReferenceTree(const char* filename, const char* treename, const char* refname){
1187 //
1188 // add a reference tree to the current tree
1189 // by default the treename is 'calPads' and the reference treename is 'R'
1190 //
1191 TFile *file = new TFile(filename);
1192 fListOfObjectsToBeDeleted->Add(file);
1193 TTree * tree = (TTree*)file->Get(treename);
1194 return AddFriend(tree, refname);
1195}
1196
1197
1198TObjArray* AliTPCCalibViewer::GetArrayOfCalPads(){
1199 //
1200 // Returns a TObjArray with all AliTPCCalPads that are stored in the tree
1201 // - this takes a while -
1202 //
1203 TObjArray *listOfCalPads = GetListOfVariables();
1204 TObjArray *calPadsArray = new TObjArray();
1205 Int_t numberOfCalPads = listOfCalPads->GetEntries();
1206 for (Int_t i = 0; i < numberOfCalPads; i++) {
1207 std::cout << "Creating calPad " << (i+1) << " of " << numberOfCalPads << "\r" << std::flush;
1208 char* calPadName = (char*)((TObjString*)(listOfCalPads->At(i)))->GetString().Data();
1209 TString drawCommand = ((TObjString*)(listOfCalPads->At(i)))->GetString();
1210 drawCommand.Append("~");
1211 AliTPCCalPad* calPad = GetCalPad(drawCommand.Data(), "", calPadName);
1212 calPadsArray->Add(calPad);
1213 }
1214 std::cout << std::endl;
1215 listOfCalPads->Delete();
1216 delete listOfCalPads;
1217 return calPadsArray;
1218}
1219
1220
a6d2bd0c 1221TString* AliTPCCalibViewer::Fit(const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, TVectorD &fitParam, TMatrixD &covMatrix){
1222 //
1223 // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
1224 // returns chi2, fitParam and covMatrix
1225 // returns TString with fitted formula
1226 //
1227
1228 TString formulaStr(formula);
1229 TString drawStr(drawCommand);
1230 TString cutStr(cuts);
1231
1232 // abbreviations:
1233 drawStr.ReplaceAll("~",".fElements");
1234 cutStr.ReplaceAll("~",".fElements");
1235 formulaStr.ReplaceAll("~", ".fElements");
1236
1237 formulaStr.ReplaceAll("++", "~");
1238 TObjArray* formulaTokens = formulaStr.Tokenize("~");
1239 Int_t dim = formulaTokens->GetEntriesFast();
1240
1241 fitParam.ResizeTo(dim);
1242 covMatrix.ResizeTo(dim,dim);
1243
1244 TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
1245 fitter->StoreData(kTRUE);
1246 fitter->ClearPoints();
1247
1248 Int_t entries = Draw(drawStr.Data(), cutStr.Data(), "goff");
1249 if (entries == -1) return new TString("An ERROR has occured during fitting!");
1250 Double_t **values = new Double_t*[dim+1] ;
1251
1252 for (Int_t i = 0; i < dim + 1; i++){
1253 Int_t centries = 0;
1254 if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff");
1255 else centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff");
1256
1257 if (entries != centries) return new TString("An ERROR has occured during fitting!");
1258 values[i] = new Double_t[entries];
1259 memcpy(values[i], fTree->GetV1(), entries*sizeof(Double_t));
1260 }
1261
1262 // add points to the fitter
1263 for (Int_t i = 0; i < entries; i++){
1264 Double_t x[1000];
1265 for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
1266 fitter->AddPoint(x, values[dim][i], 1);
1267 }
1268
1269 fitter->Eval();
1270 fitter->GetParameters(fitParam);
1271 fitter->GetCovarianceMatrix(covMatrix);
1272 chi2 = fitter->GetChisquare();
1273 chi2 = chi2;
1274
1275 TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
1276
1277 for (Int_t iparam = 0; iparam < dim; iparam++) {
1278 returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
1279 if (iparam < dim-1) returnFormula.Append("+");
1280 }
1281 returnFormula.Append(" )");
1282 delete formulaTokens;
1283 delete fitter;
1284 delete[] values;
1285 return preturnFormula;
1286}
1287
1288
39bcd65d 1289void AliTPCCalibViewer::MakeTreeWithObjects(const char * fileName, TObjArray * array, const char * mapFileName) {
1290 //
1291 // Write tree with all available information
1292 // im mapFileName is speciefied, the Map information are also written to the tree
1293 // AliTPCCalPad-Objects are written directly to the tree, so that they can be accessd later on
1294 // (does not work!!!)
1295 //
1296 AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
1297
1298 TObjArray* mapIROCs = 0;
1299 TObjArray* mapOROCs = 0;
1300 TVectorF *mapIROCArray = 0;
1301 TVectorF *mapOROCArray = 0;
1302 Int_t mapEntries = 0;
1303 TString* mapNames = 0;
1304
1305 if (mapFileName) {
1306 TFile mapFile(mapFileName, "read");
1307
1308 TList* listOfROCs = mapFile.GetListOfKeys();
1309 mapEntries = listOfROCs->GetEntries()/2;
1310 mapIROCs = new TObjArray(mapEntries*2);
1311 mapOROCs = new TObjArray(mapEntries*2);
1312 mapIROCArray = new TVectorF[mapEntries];
1313 mapOROCArray = new TVectorF[mapEntries];
1314
1315 mapNames = new TString[mapEntries];
1316 for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
72d0ab7e 1317 TString rocName(((TKey*)(listOfROCs->At(ivalue*2)))->GetName());
1318 rocName.Remove(rocName.Length()-4, 4);
1319 mapIROCs->AddAt((AliTPCCalROC*)mapFile.Get((rocName + "IROC").Data()), ivalue);
1320 mapOROCs->AddAt((AliTPCCalROC*)mapFile.Get((rocName + "OROC").Data()), ivalue);
1321 mapNames[ivalue].Append(rocName);
39bcd65d 1322 }
1323
1324 for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1325 mapIROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(0));
1326 mapOROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(36));
1327
1328 for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(0); ichannel++)
1329 (mapIROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapIROCs->At(ivalue)))->GetValue(ichannel);
1330 for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(36); ichannel++)
1331 (mapOROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapOROCs->At(ivalue)))->GetValue(ichannel);
1332 }
1333
1334 } // if (mapFileName)
1335
1336 TTreeSRedirector cstream(fileName);
1337 Int_t arrayEntries = array->GetEntries();
1338
1339 // Read names of AliTPCCalPads and save them in names[]
1340 TString* names = new TString[arrayEntries];
1341 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
1342 names[ivalue].Append(((AliTPCCalPad*)array->At(ivalue))->GetName());
1343
1344 for (UInt_t isector = 0; isector < tpcROCinstance->GetNSectors(); isector++) {
1345
1346 TVectorF *vectorArray = new TVectorF[arrayEntries];
1347 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
1348 vectorArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
1349
1350
1351 //
1352 // fill vectors of variable per pad
1353 //
1354 TVectorF *posArray = new TVectorF[8];
1355 for (Int_t ivalue = 0; ivalue < 8; ivalue++)
1356 posArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
1357
1358 Float_t posG[3] = {0};
1359 Float_t posL[3] = {0};
1360 Int_t ichannel = 0;
1361 for (UInt_t irow = 0; irow < tpcROCinstance->GetNRows(isector); irow++) {
1362 for (UInt_t ipad = 0; ipad < tpcROCinstance->GetNPads(isector, irow); ipad++) {
1363 tpcROCinstance->GetPositionLocal(isector, irow, ipad, posL);
1364 tpcROCinstance->GetPositionGlobal(isector, irow, ipad, posG);
1365 posArray[0][ichannel] = irow;
1366 posArray[1][ichannel] = ipad;
1367 posArray[2][ichannel] = posL[0];
1368 posArray[3][ichannel] = posL[1];
1369 posArray[4][ichannel] = posG[0];
1370 posArray[5][ichannel] = posG[1];
1371 posArray[6][ichannel] = (Int_t)(ipad - (Double_t)(tpcROCinstance->GetNPads(isector, irow))/2);
1372 posArray[7][ichannel] = ichannel;
1373
1374 // loop over array containing AliTPCCalPads
1375 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1376 AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue);
1377 AliTPCCalROC* calROC = calPad->GetCalROC(isector);
1378 if (calROC)
1379 (vectorArray[ivalue])[ichannel] = calROC->GetValue(irow, ipad);
1380 else
1381 (vectorArray[ivalue])[ichannel] = 0;
1382 }
1383 ichannel++;
1384 }
1385 }
1386 AliTPCCalROC dummyROC(0);
1387 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1388 AliTPCCalROC *roc = ((AliTPCCalPad*)array->At(ivalue))->GetCalROC(isector);
1389 if (!roc) roc = &dummyROC;
1390 cstream << "calPads" <<
1391 (Char_t*)((names[ivalue] + ".=").Data()) << &vectorArray[ivalue];
1392 cstream << "calPads" <<
1393 (Char_t*)((names[ivalue] + "Pad.=").Data()) << roc;
1394 }
1395
1396 if (mapFileName) {
1397 for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1398 if (isector < 36)
1399 cstream << "calPads" <<
1400 (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapIROCArray[ivalue];
1401 else
1402 cstream << "calPads" <<
1403 (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapOROCArray[ivalue];
1404 }
1405 }
1406
1407 cstream << "calPads" <<
1408 "sector=" << isector;
1409
1410 cstream << "calPads" <<
1411 "row.=" << &posArray[0] <<
1412 "pad.=" << &posArray[1] <<
1413 "lx.=" << &posArray[2] <<
1414 "ly.=" << &posArray[3] <<
1415 "gx.=" << &posArray[4] <<
1416 "gy.=" << &posArray[5] <<
1417 "rpad.=" << &posArray[6] <<
1418 "channel.=" << &posArray[7];
1419
1420 cstream << "calPads" <<
1421 "\n";
1422
1423 delete[] posArray;
1424 delete[] vectorArray;
1425 } //for (UInt_t isector = 0; isector < tpcROCinstance->GetNSectors(); isector++)
1426
1427 delete[] names;
1428 if (mapFileName) {
1429 delete mapIROCs;
1430 delete mapOROCs;
1431 delete[] mapIROCArray;
1432 delete[] mapOROCArray;
1433 delete[] mapNames;
1434 }
1435}
1436
72d0ab7e 1437
39bcd65d 1438void AliTPCCalibViewer::MakeTree(const char * fileName, TObjArray * array, const char * mapFileName, AliTPCCalPad* outlierPad, Float_t ltmFraction) {
1439 //
1440 // Write a tree with all available information
72d0ab7e 1441 // if mapFileName is speciefied, the Map information are also written to the tree
39bcd65d 1442 // pads specified in outlierPad are not used for calculating statistics
72d0ab7e 1443 // The following statistical information on the basis of a ROC are calculated:
1444 // "_Median", "_Mean", "_LTM", "_RMS_LTM"
1445 // "_Median_OutlierCutted", "_Mean_OutlierCutted", "_RMS_OutlierCutted", "_LTM_OutlierCutted", "_RMS_LTM_OutlierCutted"
1446 // The following position variables are available:
1447 // "row", "pad", "lx", "ly", "gx", "gy", "rpad", "channel"
1448 //
1449 // The tree out of this function is the basis for the AliTPCCalibViewer and the AliTPCCalibViewerGUI.
1450
39bcd65d 1451 AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
1452
1453 TObjArray* mapIROCs = 0;
1454 TObjArray* mapOROCs = 0;
1455 TVectorF *mapIROCArray = 0;
1456 TVectorF *mapOROCArray = 0;
1457 Int_t mapEntries = 0;
1458 TString* mapNames = 0;
1459
1460 if (mapFileName) {
1461 TFile mapFile(mapFileName, "read");
1462
1463 TList* listOfROCs = mapFile.GetListOfKeys();
1464 mapEntries = listOfROCs->GetEntries()/2;
1465 mapIROCs = new TObjArray(mapEntries*2);
1466 mapOROCs = new TObjArray(mapEntries*2);
1467 mapIROCArray = new TVectorF[mapEntries];
1468 mapOROCArray = new TVectorF[mapEntries];
1469
1470 mapNames = new TString[mapEntries];
1471 for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
72d0ab7e 1472 TString rocName(((TKey*)(listOfROCs->At(ivalue*2)))->GetName());
1473 rocName.Remove(rocName.Length()-4, 4);
1474 mapIROCs->AddAt((AliTPCCalROC*)mapFile.Get((rocName + "IROC").Data()), ivalue);
1475 mapOROCs->AddAt((AliTPCCalROC*)mapFile.Get((rocName + "OROC").Data()), ivalue);
1476 mapNames[ivalue].Append(rocName);
39bcd65d 1477 }
1478
1479 for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1480 mapIROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(0));
1481 mapOROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(36));
1482
1483 for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(0); ichannel++)
1484 (mapIROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapIROCs->At(ivalue)))->GetValue(ichannel);
1485 for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(36); ichannel++)
1486 (mapOROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapOROCs->At(ivalue)))->GetValue(ichannel);
1487 }
1488
1489 } // if (mapFileName)
1490
1491 TTreeSRedirector cstream(fileName);
72d0ab7e 1492 Int_t arrayEntries = 0;
1493 if (array) arrayEntries = array->GetEntries();
39bcd65d 1494
1495 TString* names = new TString[arrayEntries];
1496 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
1497 names[ivalue].Append(((AliTPCCalPad*)array->At(ivalue))->GetName());
1498
1499 for (UInt_t isector = 0; isector < tpcROCinstance->GetNSectors(); isector++) {
1500 //
1501 // get statistic for given sector
1502 //
1503 TVectorF median(arrayEntries);
1504 TVectorF mean(arrayEntries);
1505 TVectorF rms(arrayEntries);
1506 TVectorF ltm(arrayEntries);
1507 TVectorF ltmrms(arrayEntries);
1508 TVectorF medianWithOut(arrayEntries);
1509 TVectorF meanWithOut(arrayEntries);
1510 TVectorF rmsWithOut(arrayEntries);
1511 TVectorF ltmWithOut(arrayEntries);
1512 TVectorF ltmrmsWithOut(arrayEntries);
1513
1514 TVectorF *vectorArray = new TVectorF[arrayEntries];
1515 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
1516 vectorArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
1517
1518 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1519 AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue);
1520 AliTPCCalROC* calROC = calPad->GetCalROC(isector);
1521 AliTPCCalROC* outlierROC = 0;
1522 if (outlierPad) outlierROC = outlierPad->GetCalROC(isector);
1523 if (calROC) {
1524 median[ivalue] = calROC->GetMedian();
1525 mean[ivalue] = calROC->GetMean();
1526 rms[ivalue] = calROC->GetRMS();
1527 Double_t ltmrmsValue = 0;
1528 ltm[ivalue] = calROC->GetLTM(&ltmrmsValue, ltmFraction);
1529 ltmrms[ivalue] = ltmrmsValue;
1530 if (outlierROC) {
1531 medianWithOut[ivalue] = calROC->GetMedian(outlierROC);
1532 meanWithOut[ivalue] = calROC->GetMean(outlierROC);
1533 rmsWithOut[ivalue] = calROC->GetRMS(outlierROC);
1534 ltmrmsValue = 0;
1535 ltmWithOut[ivalue] = calROC->GetLTM(&ltmrmsValue, ltmFraction, outlierROC);
1536 ltmrmsWithOut[ivalue] = ltmrmsValue;
1537 }
1538 }
1539 else {
1540 median[ivalue] = 0.;
1541 mean[ivalue] = 0.;
1542 rms[ivalue] = 0.;
1543 ltm[ivalue] = 0.;
1544 ltmrms[ivalue] = 0.;
1545 medianWithOut[ivalue] = 0.;
1546 meanWithOut[ivalue] = 0.;
1547 rmsWithOut[ivalue] = 0.;
1548 ltmWithOut[ivalue] = 0.;
1549 ltmrmsWithOut[ivalue] = 0.;
1550 }
1551 }
1552
1553 //
1554 // fill vectors of variable per pad
1555 //
1556 TVectorF *posArray = new TVectorF[8];
1557 for (Int_t ivalue = 0; ivalue < 8; ivalue++)
1558 posArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
1559
1560 Float_t posG[3] = {0};
1561 Float_t posL[3] = {0};
1562 Int_t ichannel = 0;
1563 for (UInt_t irow = 0; irow < tpcROCinstance->GetNRows(isector); irow++) {
1564 for (UInt_t ipad = 0; ipad < tpcROCinstance->GetNPads(isector, irow); ipad++) {
1565 tpcROCinstance->GetPositionLocal(isector, irow, ipad, posL);
1566 tpcROCinstance->GetPositionGlobal(isector, irow, ipad, posG);
1567 posArray[0][ichannel] = irow;
1568 posArray[1][ichannel] = ipad;
1569 posArray[2][ichannel] = posL[0];
1570 posArray[3][ichannel] = posL[1];
1571 posArray[4][ichannel] = posG[0];
1572 posArray[5][ichannel] = posG[1];
1573 posArray[6][ichannel] = (Int_t)(ipad - (Double_t)(tpcROCinstance->GetNPads(isector, irow))/2);
1574 posArray[7][ichannel] = ichannel;
1575
1576 // loop over array containing AliTPCCalPads
1577 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1578 AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue);
1579 AliTPCCalROC* calROC = calPad->GetCalROC(isector);
1580 if (calROC)
1581 (vectorArray[ivalue])[ichannel] = calROC->GetValue(irow, ipad);
1582 else
1583 (vectorArray[ivalue])[ichannel] = 0;
1584 }
1585 ichannel++;
1586 }
1587 }
1588
1589 cstream << "calPads" <<
1590 "sector=" << isector;
1591
1592 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1593 cstream << "calPads" <<
1594 (Char_t*)((names[ivalue] + "_Median=").Data()) << median[ivalue] <<
1595 (Char_t*)((names[ivalue] + "_Mean=").Data()) << mean[ivalue] <<
1596 (Char_t*)((names[ivalue] + "_RMS=").Data()) << rms[ivalue] <<
1597 (Char_t*)((names[ivalue] + "_LTM=").Data()) << ltm[ivalue] <<
1598 (Char_t*)((names[ivalue] + "_RMS_LTM=").Data()) << ltmrms[ivalue];
1599 if (outlierPad) {
1600 cstream << "calPads" <<
1601 (Char_t*)((names[ivalue] + "_Median_OutlierCutted=").Data()) << medianWithOut[ivalue] <<
1602 (Char_t*)((names[ivalue] + "_Mean_OutlierCutted=").Data()) << meanWithOut[ivalue] <<
1603 (Char_t*)((names[ivalue] + "_RMS_OutlierCutted=").Data()) << rmsWithOut[ivalue] <<
1604 (Char_t*)((names[ivalue] + "_LTM_OutlierCutted=").Data()) << ltmWithOut[ivalue] <<
1605 (Char_t*)((names[ivalue] + "_RMS_LTM_OutlierCutted=").Data()) << ltmrmsWithOut[ivalue];
1606 }
1607 }
1608
1609 for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
1610 cstream << "calPads" <<
1611 (Char_t*)((names[ivalue] + ".=").Data()) << &vectorArray[ivalue];
1612 }
1613
1614 if (mapFileName) {
1615 for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
1616 if (isector < 36)
1617 cstream << "calPads" <<
1618 (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapIROCArray[ivalue];
1619 else
1620 cstream << "calPads" <<
1621 (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapOROCArray[ivalue];
1622 }
1623 }
1624
1625 cstream << "calPads" <<
1626 "row.=" << &posArray[0] <<
1627 "pad.=" << &posArray[1] <<
1628 "lx.=" << &posArray[2] <<
1629 "ly.=" << &posArray[3] <<
1630 "gx.=" << &posArray[4] <<
1631 "gy.=" << &posArray[5] <<
1632 "rpad.=" << &posArray[6] <<
1633 "channel.=" << &posArray[7];
1634
1635 cstream << "calPads" <<
1636 "\n";
1637
1638 delete[] posArray;
1639 delete[] vectorArray;
1640 }
1641
1642
1643 delete[] names;
1644 if (mapFileName) {
1645 delete mapIROCs;
1646 delete mapOROCs;
1647 delete[] mapIROCArray;
1648 delete[] mapOROCArray;
1649 delete[] mapNames;
1650 }
1651}
1652