]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix for infinite loop.
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTracks.cxx
CommitLineData
10757ee9 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 to analyse tracks for calibration //
20// to be used as a component in AliTPCSelectorTracks //
21// In the constructor you have to specify name and title //
22// to get the Object out of a file. //
23// The parameter 'clusterParam', a AliTPCClusterParam object //
24// (needed for TPC cluster error and shape parameterization) //
25// Normally you get this object out of the file 'TPCClusterParam.root' //
26// In the parameter 'cuts' the cuts are specified, that decide //
27// weather a track will be accepted for calibration or not. //
28// //
29//
30// The data flow:
31//
32/*
33 Raw Data -> Local Reconstruction -> Tracking -> Calibration -> RefData (component itself)
34 Offline/HLT Offline/HLT OCDB entries (AliTPCClusterParam)
35*/
36
9f0beaf7 37/*
38
39How to retrive it from file (created using calibration task):
40
41gSystem->Load("libANALYSIS");
42gSystem->Load("libTPCcalib");
43TFile fcalib("CalibObjects.root");
44TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
45AliTPCcalibTracks * calibTracks = ( AliTPCcalibTracks *)array->FindObject("calibTracks");
46
9c171b96 47
48//USAGE of debug stream example
49 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
50 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
51 AliXRDPROOFtoolkit tool;
52 TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
53 chainres->Lookup();
9f0beaf7 54*/
55
56
10757ee9 57 //
58///////////////////////////////////////////////////////////////////////////////
59
60//
61// ROOT includes
62//
63#include <iostream>
64#include <fstream>
65using namespace std;
66
67#include <TROOT.h>
68#include <TChain.h>
69#include <TFile.h>
70#include <TH3F.h>
5b00528f 71#include <TProfile.h>
72
10757ee9 73//
74//#include <TPDGCode.h>
75#include <TStyle.h>
76#include "TLinearFitter.h"
77//#include "TMatrixD.h"
78#include "TTreeStream.h"
79#include "TF1.h"
80#include <TCanvas.h>
81#include <TGraph2DErrors.h>
82#include "TPostScript.h"
83#include "TCint.h"
84
85#include <TH2D.h>
86#include <TF2.h>
87#include <TSystem.h>
88#include <TCollection.h>
89#include <iostream>
90#include <TLinearFitter.h>
2e5bcb67 91#include <TString.h>
10757ee9 92
93//
94// AliROOT includes
95//
96#include "AliMagF.h"
97#include "AliTracker.h"
98#include "AliESD.h"
99#include "AliESDtrack.h"
100#include "AliESDfriend.h"
101#include "AliESDfriendTrack.h"
102#include "AliTPCseed.h"
103#include "AliTPCclusterMI.h"
104#include "AliTPCROC.h"
105
106#include "AliTPCParamSR.h"
107#include "AliTrackPointArray.h"
108#include "AliTPCcalibTracks.h"
109#include "AliTPCClusterParam.h"
110#include "AliTPCcalibTracksCuts.h"
111#include "AliTPCCalPadRegion.h"
112#include "AliTPCCalPad.h"
113#include "AliTPCCalROC.h"
114#include "TText.h"
115#include "TPaveText.h"
116#include "TSystem.h"
2e5bcb67 117#include "TStatToolkit.h"
118#include "TCut.h"
af6a50bb 119#include "THnSparse.h"
10757ee9 120
10757ee9 121
122
123ClassImp(AliTPCcalibTracks)
124
125
126AliTPCcalibTracks::AliTPCcalibTracks():
c32da879 127 AliTPCcalibBase(),
2c632057 128 fClusterParam(0),
2c632057 129 fROC(0),
af6a50bb 130 fHisDeltaY(0), // THnSparse - delta Y
131 fHisDeltaZ(0), // THnSparse - delta Z
132 fHisRMSY(0), // THnSparse - rms Y
133 fHisRMSZ(0), // THnSparse - rms Z
134 fHisQmax(0), // THnSparse - qmax
135 fHisQtot(0), // THnSparse - qtot
2c632057 136 fArrayAmpRow(0),
137 fArrayAmp(0),
138 fArrayQDY(0),
139 fArrayQDZ(0),
140 fArrayQRMSY(0),
141 fArrayQRMSZ(0),
142 fArrayChargeVsDriftlength(0),
143 fcalPadRegionChargeVsDriftlength(0),
144 fDeltaY(0),
145 fDeltaZ(0),
146 fResolY(0),
147 fResolZ(0),
148 fRMSY(0),
149 fRMSZ(0),
150 fCuts(0),
151 fHclus(0),
152 fRejectedTracksHisto(0),
153 fHclusterPerPadrow(0),
154 fHclusterPerPadrowRaw(0),
155 fClusterCutHisto(0),
156 fCalPadClusterPerPad(0),
9f0beaf7 157 fCalPadClusterPerPadRaw(0)
10757ee9 158{
159 //
160 // AliTPCcalibTracks default constructor
161 //
ae28e92e 162 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;
10757ee9 163}
164
165
b42cf9ed 166
167AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
ae28e92e 168 AliTPCcalibBase(calibTracks),
2c632057 169 fClusterParam(0),
2c632057 170 fROC(0),
af6a50bb 171 fHisDeltaY(0), // THnSparse - delta Y
172 fHisDeltaZ(0), // THnSparse - delta Z
173 fHisRMSY(0), // THnSparse - rms Y
174 fHisRMSZ(0), // THnSparse - rms Z
175 fHisQmax(0), // THnSparse - qmax
176 fHisQtot(0), // THnSparse - qtot
2c632057 177 fArrayAmpRow(0),
178 fArrayAmp(0),
179 fArrayQDY(0),
180 fArrayQDZ(0),
181 fArrayQRMSY(0),
182 fArrayQRMSZ(0),
183 fArrayChargeVsDriftlength(0),
184 fcalPadRegionChargeVsDriftlength(0),
185 fDeltaY(0),
186 fDeltaZ(0),
187 fResolY(0),
188 fResolZ(0),
189 fRMSY(0),
190 fRMSZ(0),
191 fCuts(0),
192 fHclus(0),
193 fRejectedTracksHisto(0),
194 fHclusterPerPadrow(0),
195 fHclusterPerPadrowRaw(0),
196 fClusterCutHisto(0),
197 fCalPadClusterPerPad(0),
9f0beaf7 198 fCalPadClusterPerPadRaw(0)
2c632057 199{
10757ee9 200 //
201 // AliTPCcalibTracks copy constructor
202 //
ae28e92e 203 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
10757ee9 204
205 Bool_t dirStatus = TH1::AddDirectoryStatus();
206 TH1::AddDirectory(kFALSE);
207
208 Int_t length = -1;
209 // backward compatibility: if the data member doesn't yet exist, it will not be merged
b42cf9ed 210 (calibTracks.fArrayAmpRow) ? length = calibTracks.fArrayAmpRow->GetEntriesFast() : length = -1;
10757ee9 211 fArrayAmpRow = new TObjArray(length);
212 fArrayAmp = new TObjArray(length);
213 for (Int_t i = 0; i < length; i++) {
b42cf9ed 214 fArrayAmpRow->AddAt( (TProfile*)calibTracks.fArrayAmpRow->At(i)->Clone(), i);
215 fArrayAmp->AddAt( ((TProfile*)calibTracks.fArrayAmp->At(i)->Clone()), i);
10757ee9 216 }
217
b42cf9ed 218 (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
10757ee9 219 fArrayQDY= new TObjArray(length);
220 fArrayQDZ= new TObjArray(length);
221 fArrayQRMSY= new TObjArray(length);
222 fArrayQRMSZ= new TObjArray(length);
223 for (Int_t i = 0; i < length; i++) {
b42cf9ed 224 fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
225 fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
226 fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
227 fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
10757ee9 228 }
229
b42cf9ed 230 (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
10757ee9 231 fResolY = new TObjArray(length);
232 fResolZ = new TObjArray(length);
233 fRMSY = new TObjArray(length);
234 fRMSZ = new TObjArray(length);
235 for (Int_t i = 0; i < length; i++) {
b42cf9ed 236 fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
237 fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
238 fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
239 fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
10757ee9 240 }
241
b42cf9ed 242 (calibTracks.fArrayChargeVsDriftlength) ? length = calibTracks.fArrayChargeVsDriftlength->GetEntriesFast() : length = -1;
243 (calibTracks.fArrayChargeVsDriftlength) ? fArrayChargeVsDriftlength = new TObjArray(length) : fArrayChargeVsDriftlength = 0;
10757ee9 244 for (Int_t i = 0; i < length; i++) {
b42cf9ed 245 fArrayChargeVsDriftlength->AddAt( ((TProfile*)calibTracks.fArrayChargeVsDriftlength->At(i)->Clone()), i);
10757ee9 246 }
247
b42cf9ed 248 fDeltaY = (TH1F*)calibTracks.fDeltaY->Clone();
249 fDeltaZ = (TH1F*)calibTracks.fDeltaZ->Clone();
250 fHclus = (TH1I*)calibTracks.fHclus->Clone();
251 fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
252 fRejectedTracksHisto = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
253 fHclusterPerPadrow = (TH1I*)calibTracks.fHclusterPerPadrow->Clone();
254 fHclusterPerPadrowRaw = (TH1I*)calibTracks.fHclusterPerPadrowRaw->Clone();
255 fcalPadRegionChargeVsDriftlength = (AliTPCCalPadRegion*)calibTracks.fcalPadRegionChargeVsDriftlength->Clone();
256 fCalPadClusterPerPad = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
257 fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
258
259 fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(),
260 calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
b42cf9ed 261 SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
10757ee9 262 TH1::AddDirectory(dirStatus); // set status back to original status
263// cout << "+++++ end of copy constructor +++++" << endl; // TO BE REMOVED
264}
265
266
b42cf9ed 267AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
268 //
269 // assgnment operator
270 //
271 if (this != &calibTracks) {
272 new (this) AliTPCcalibTracks(calibTracks);
273 }
274 return *this;
275
276}
277
278
10757ee9 279AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts, Int_t logLevel) :
c32da879 280 AliTPCcalibBase(),
2c632057 281 fClusterParam(0),
2c632057 282 fROC(0),
af6a50bb 283 fHisDeltaY(0), // THnSparse - delta Y
284 fHisDeltaZ(0), // THnSparse - delta Z
285 fHisRMSY(0), // THnSparse - rms Y
286 fHisRMSZ(0), // THnSparse - rms Z
287 fHisQmax(0), // THnSparse - qmax
288 fHisQtot(0), // THnSparse - qtot
2c632057 289 fArrayAmpRow(0),
290 fArrayAmp(0),
291 fArrayQDY(0),
292 fArrayQDZ(0),
293 fArrayQRMSY(0),
294 fArrayQRMSZ(0),
295 fArrayChargeVsDriftlength(0),
296 fcalPadRegionChargeVsDriftlength(0),
297 fDeltaY(0),
298 fDeltaZ(0),
299 fResolY(0),
300 fResolZ(0),
301 fRMSY(0),
302 fRMSZ(0),
303 fCuts(0),
304 fHclus(0),
305 fRejectedTracksHisto(0),
306 fHclusterPerPadrow(0),
307 fHclusterPerPadrowRaw(0),
308 fClusterCutHisto(0),
309 fCalPadClusterPerPad(0),
9f0beaf7 310 fCalPadClusterPerPadRaw(0)
2c632057 311 {
10757ee9 312 //
313 // AliTPCcalibTracks constructor
314 // specify 'name' and 'title' of your object
315 // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
316 // In the parameter 'cuts' the cuts are specified, that decide
317 // weather a track will be accepted for calibration or not.
ae28e92e 318 //
319 // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
10757ee9 320 //
321 // All histograms are instatiated in this constructor.
322 //
c32da879 323 this->SetName(name);
324 this->SetTitle(title);
325
ae28e92e 326 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
10757ee9 327 G__SetCatchException(0);
328
329 fClusterParam = clusterParam;
330 if (fClusterParam){
331 fClusterParam->SetInstance(fClusterParam);
332 }
333 else {
334 Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
335 }
336 fCuts = cuts;
ae28e92e 337 SetDebugLevel(logLevel);
af6a50bb 338 MakeHistos();
10757ee9 339
340 TH1::AddDirectory(kFALSE);
341
342 char chname[1000];
343 TProfile * prof1=0;
344 TH1F * his1 =0;
345 fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160); // valgrind 3
8b9b3187 346 fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
10757ee9 347 fHclusterPerPadrow = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160);
348 fHclusterPerPadrowRaw = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160);
349 fCalPadClusterPerPad = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
350 fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
351 fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
352
353 // Amplitude - sector - row histograms
354 fArrayAmpRow = new TObjArray(72);
355 fArrayAmp = new TObjArray(72);
356 fArrayChargeVsDriftlength = new TObjArray(72);
357
358 for (Int_t i = 0; i < 36; i++){
359 sprintf(chname,"Amp_row_Sector%d",i);
360 prof1 = new TProfile(chname,chname,63,0,64); // valgrind 3 193,536 bytes in 354 blocks are still reachable
361 prof1->SetXTitle("Pad row");
362 prof1->SetYTitle("Mean Max amplitude");
363 fArrayAmpRow->AddAt(prof1,i);
364 sprintf(chname,"Amp_row_Sector%d",i+36);
365 prof1 = new TProfile(chname,chname,96,0,97); // valgrind 3 3,912 bytes in 6 blocks are possibly lost
366 prof1->SetXTitle("Pad row");
367 prof1->SetYTitle("Mean Max amplitude");
368 fArrayAmpRow->AddAt(prof1,i+36);
369
370 // amplitude
371 sprintf(chname,"Amp_Sector%d",i);
15e48021 372 his1 = new TH1F(chname,chname,100,0,32); // valgrind
10757ee9 373 his1->SetXTitle("Max Amplitude (ADC)");
374 fArrayAmp->AddAt(his1,i);
375 sprintf(chname,"Amp_Sector%d",i+36);
15e48021 376 his1 = new TH1F(chname,chname,100,0,32); // valgrind 3 13,408,208 bytes in 229 blocks are still reachable
10757ee9 377 his1->SetXTitle("Max Amplitude (ADC)");
378 fArrayAmp->AddAt(his1,i+36);
379
380 // driftlength
381 sprintf(chname, "driftlengt vs. charge, ROC %i", i);
15e48021 382 prof1 = new TProfile(chname, chname, 25, 0, 250);
10757ee9 383 prof1->SetYTitle("Charge");
384 prof1->SetXTitle("Driftlength");
385 fArrayChargeVsDriftlength->AddAt(prof1,i);
386 sprintf(chname, "driftlengt vs. charge, ROC %i", i+36);
15e48021 387 prof1 = new TProfile(chname, chname, 25, 0, 250);
10757ee9 388 prof1->SetYTitle("Charge");
389 prof1->SetXTitle("Driftlength");
390 fArrayChargeVsDriftlength->AddAt(prof1,i+36);
391 }
392
393 TH1::AddDirectory(kFALSE);
394
395 fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1);
396 fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1);
397
398 fResolY = new TObjArray(3);
399 fResolZ = new TObjArray(3);
400 fRMSY = new TObjArray(3);
401 fRMSZ = new TObjArray(3);
402 TH3F * his3D;
403 //
404 his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
405 fResolY->AddAt(his3D,0);
406 his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
407 fResolY->AddAt(his3D,1);
408 his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
409 fResolY->AddAt(his3D,2);
410 //
411 his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
412 fResolZ->AddAt(his3D,0);
413 his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
414 fResolZ->AddAt(his3D,1);
415 his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
416 fResolZ->AddAt(his3D,2);
417 //
418 his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
419 fRMSY->AddAt(his3D,0);
420 his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
421 fRMSY->AddAt(his3D,1);
422 his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
423 fRMSY->AddAt(his3D,2);
424 //
425 his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
426 fRMSZ->AddAt(his3D,0);
427 his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
428 fRMSZ->AddAt(his3D,1);
429 his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
430 fRMSZ->AddAt(his3D,2);
431 //
432
433 TH1::AddDirectory(kFALSE);
434
435 fArrayQDY = new TObjArray(300);
436 fArrayQDZ = new TObjArray(300);
437 fArrayQRMSY = new TObjArray(300);
438 fArrayQRMSZ = new TObjArray(300);
439 for (Int_t iq = 0; iq <= 10; iq++){
440 for (Int_t ipad = 0; ipad < 3; ipad++){
441 Int_t bin = GetBin(iq, ipad);
442 Float_t qmean = GetQ(bin);
d23650ed 443 char hname[200];
444 sprintf(hname,"ResolY Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 445 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
10757ee9 446 fArrayQDY->AddAt(his3D, bin);
d23650ed 447 sprintf(hname,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 448 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
10757ee9 449 fArrayQDZ->AddAt(his3D, bin);
d23650ed 450 sprintf(hname,"RMSY Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 451 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
10757ee9 452 fArrayQRMSY->AddAt(his3D, bin);
d23650ed 453 sprintf(hname,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 454 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
10757ee9 455 fArrayQRMSZ->AddAt(his3D, bin);
456 }
457 }
458
459 fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region");
460 TProfile *tempProf;
461 for (UInt_t padSize = 0; padSize < 3; padSize++) {
462 for (UInt_t isector = 0; isector < 36; isector++) {
463 if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector);
464 if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium pads", isector);
465 if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long pads", isector);
466 tempProf = new TProfile(chname, chname, 500, 0, 250);
467 tempProf->SetYTitle("Charge");
468 tempProf->SetXTitle("Driftlength");
469 fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize);
470 }
471 }
472
10757ee9 473
ae28e92e 474 if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl;
10757ee9 475 cout << "end of main constructor" << endl; // TO BE REMOVED
476}
477
478
479AliTPCcalibTracks::~AliTPCcalibTracks() {
480 //
481 // AliTPCcalibTracks destructor
482 //
483
ae28e92e 484 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
10757ee9 485 Int_t length = 0;
486 if (fArrayAmpRow) length = fArrayAmpRow->GetEntriesFast();
487 for (Int_t i = 0; i < length; i++){
488 delete fArrayAmpRow->At(i);
489 delete fArrayAmp->At(i);
490 }
491 delete fArrayAmpRow;
492 delete fArrayAmp;
493
494 delete fDeltaY;
495 delete fDeltaZ;
496
497 if (fResolY) length = fResolY->GetEntriesFast();
498 for (Int_t i = 0; i < length; i++){
499 delete fResolY->At(i);
500 delete fResolZ->At(i);
501 delete fRMSY->At(i);
502 delete fRMSZ->At(i);
503 }
504 delete fResolY;
505 delete fResolZ;
506 delete fRMSY;
507 delete fRMSZ;
508
509 if (fArrayQDY) length = fArrayQDY->GetEntriesFast();
510 for (Int_t i = 0; i < length; i++){
511 delete fArrayQDY->At(i);
512 delete fArrayQDZ->At(i);
513 delete fArrayQRMSY->At(i);
514 delete fArrayQRMSZ->At(i);
515 }
516
517 if (fArrayChargeVsDriftlength) length = fArrayChargeVsDriftlength->GetEntriesFast();
518 for (Int_t i = 0; i < length; i++){
519 delete fArrayChargeVsDriftlength->At(i);
520 }
521
9f0beaf7 522
10757ee9 523 delete fArrayQDY;
524 delete fArrayQDZ;
525 delete fArrayQRMSY;
526 delete fArrayQRMSZ;
527 delete fArrayChargeVsDriftlength;
528
529 delete fHclus;
530 delete fRejectedTracksHisto;
531 delete fClusterCutHisto;
532 delete fHclusterPerPadrow;
533 delete fHclusterPerPadrowRaw;
534 if (fCalPadClusterPerPad) delete fCalPadClusterPerPad;
535 if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
9a214fd5 536 if(fcalPadRegionChargeVsDriftlength) {
537 fcalPadRegionChargeVsDriftlength->Delete();
538 delete fcalPadRegionChargeVsDriftlength;
539 }
af6a50bb 540 delete fHisDeltaY; // THnSparse - delta Y
541 delete fHisDeltaZ; // THnSparse - delta Z
542 delete fHisRMSY; // THnSparse - rms Y
543 delete fHisRMSZ; // THnSparse - rms Z
544 delete fHisQmax; // THnSparse - qmax
545 delete fHisQtot; // THnSparse - qtot
546
10757ee9 547}
548
549
10757ee9 550
c32da879 551void AliTPCcalibTracks::Process(AliTPCseed *track){
10757ee9 552 //
553 // To be called in the selector
554 // first AcceptTrack is evaluated, then calls all the following analyse functions:
555 // FillResolutionHistoLocal(track)
af6a50bb 556
10757ee9 557 //
ae28e92e 558 if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
10757ee9 559 Int_t accpetStatus = AcceptTrack(track);
560 if (accpetStatus == 0) {
561 FillResolutionHistoLocal(track);
10757ee9 562 }
563 else fRejectedTracksHisto->Fill(accpetStatus);
564}
565
566
567
568Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
569 //
570 // calculate bins for given q and pad type
571 // used in TObjArray
572 //
573 Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );
574 res *= 3;
575 res += pad;
576 return res;
577}
578
579
580Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
581 //
582 // calculate bins for given iq and pad type
583 // used in TObjArray
584 //
585 return iq * 3 + pad;;
586}
587
588
589Float_t AliTPCcalibTracks::GetQ(Int_t bin){
590 //
591 // returns to bin belonging charge
592 // (bin / 3 + 3)^2
593 //
594 Int_t bin0 = bin / 3;
595 bin0 += 3;
596 return bin0 * bin0;
597}
598
599
600Float_t AliTPCcalibTracks::GetPad(Int_t bin){
601 //
602 // returns to bin belonging pad
603 // bin % 3
604 //
605 return bin % 3;
606}
607
608
609
610Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
611 //
612 // Function, that decides wheather a given track is accepted for
613 // the analysis or not.
614 // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
615 // Returns 0 if a track is accepted or an integer different from 0
616 // to indicate the failed cut
617 //
618 const Int_t kMinClusters = fCuts->GetMinClusters();
619 const Float_t kMinRatio = fCuts->GetMinRatio();
620 const Float_t kMax1pt = fCuts->GetMax1pt();
621 const Float_t kEdgeYXCutNoise = fCuts->GetEdgeYXCutNoise();
622 const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
623
624 //
625 // edge induced noise tracks - NEXT RELEASE will be removed during tracking
626 if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
627 if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
628 if (track->GetNumberOfClusters() < kMinClusters) return 2;
629 Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
630 if (ratio < kMinRatio) return 3;
631// Float_t mpt = track->Get1Pt(); // Get1Pt() doesn't exist any more
632 Float_t mpt = track->GetSigned1Pt();
633 if (TMath::Abs(mpt) > kMax1pt) return 4;
634 //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
635 //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
636 //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
637
430a3403 638 if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");
10757ee9 639 return 0;
640}
641
642
643void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
644 //
645 // fill resolution histograms - localy - tracklet in the neighborhood
646 // write debug information to 'TPCSelectorDebug.root'
647 //
648 // _ the main function, called during track analysis _
649 //
650 // loop over all padrows along the track
651 // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
652 //
653 // loop again over all padrows along the track
654 // fit tracklet (clusters in given padrow +- kDelta padrows)
655 // with polynom of 2nd order and two polynoms of 1st order
656 // take both polynoms of 1st order, calculate difference of their parameters
657 // add covariance matrixes and calculate chi2 of this difference
658 // if this chi2 is bigger than a given threshold, assume that the current cluster is
659 // a kink an goto next padrow
660 // if not:
661 // fill fArrayAmpRow, array with amplitudes vs. row for given sector
662 // fill fArrayAmp, array with amplitude histograms for give sector
663 // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fDeltaY, fDeltaZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
664 //
665 // write debug information to 'TPCSelectorDebug.root'
666 // only for every kDeltaWriteDebugStream'th padrow to reduce data volume
667 // and to avoid redundant data
668 //
669
9f0beaf7 670 static TLinearFitter fFitterLinY1(2,"pol1"); //
671 static TLinearFitter fFitterLinZ1(2,"pol1"); //
672 static TLinearFitter fFitterLinY2(2,"pol1"); //
673 static TLinearFitter fFitterLinZ2(2,"pol1"); //
674 static TLinearFitter fFitterParY(3,"pol2"); //
675 static TLinearFitter fFitterParZ(3,"pol2"); //
676
677 fFitterLinY1.StoreData(kFALSE);
678 fFitterLinZ1.StoreData(kFALSE);
679 fFitterLinY2.StoreData(kFALSE);
680 fFitterLinZ2.StoreData(kFALSE);
681 fFitterParY.StoreData(kFALSE);
682 fFitterParZ.StoreData(kFALSE);
683
684
ae28e92e 685 if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
10757ee9 686 const Int_t kDelta = 10; // delta rows to fit
687 const Float_t kMinRatio = 0.75; // minimal ratio
2e5bcb67 688 // const Float_t kCutChi2 = 6.; // cut chi2 - left right - kink removal
10757ee9 689 const Float_t kErrorFraction = 0.5; // use only clusters with small interpolation error - for error param
690 const Int_t kFirstLargePad = 127; // medium pads -> long pads
691 const Float_t kLargePadSize = 1.5; // factor between medium and long pads' area
692 const Int_t kDeltaWriteDebugStream = 5; // only for every kDeltaWriteDebugStream'th padrow debug information is calulated and written to debugstream
10757ee9 693 TVectorD paramY0(2);
694 TVectorD paramZ0(2);
695 TVectorD paramY1(2);
696 TVectorD paramZ1(2);
697 TVectorD paramY2(3);
698 TVectorD paramZ2(3);
699 TMatrixD matrixY0(2,2);
700 TMatrixD matrixZ0(2,2);
701 TMatrixD matrixY1(2,2);
702 TMatrixD matrixZ1(2,2);
703
704 // estimate mean error
705 Int_t nTrackletsAll = 0; // number of tracklets for given track
706 Float_t csigmaY = 0; // mean sigma for tracklet refit in Y direction
707 Float_t csigmaZ = 0; // mean sigma for tracklet refit in Z direction
708 Int_t nClusters = 0; // working variable, number of clusters per tracklet
709 Int_t sectorG = -1; // working variable, sector of tracklet, has to stay constant for one tracklet
710
711 fHclus->Fill(track->GetNumberOfClusters()); // for statistics overview
712 // ---------------------------------------------------------------------
713 for (Int_t irow = 0; irow < 159; irow++){
714 // loop over all rows along the track
715 // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
716 // calculate mean chi^2 for this track-fit in Y and Z direction
717 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
718 if (!cluster0) continue; // no cluster found
719 Int_t sector = cluster0->GetDetector();
720 fHclusterPerPadrowRaw->Fill(irow);
721
722 Int_t ipad = TMath::Nint(cluster0->GetPad());
723 Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
724 fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
725
726 if (sector != sectorG){
727 // track leaves sector before it crossed enough rows to fit / initialization
728 nClusters = 0;
9f0beaf7 729 fFitterParY.ClearPoints();
730 fFitterParZ.ClearPoints();
10757ee9 731 sectorG = sector;
732 }
733 else {
734 nClusters++;
735 Double_t x = cluster0->GetX();
9f0beaf7 736 fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
737 fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
10757ee9 738 //
739 if ( nClusters >= kDelta + 3 ){
740 // if more than 13 (kDelta+3) clusters were added to the fitters
741 // fit the tracklet, increase trackletCounter
9f0beaf7 742 fFitterParY.Eval();
743 fFitterParZ.Eval();
10757ee9 744 nTrackletsAll++;
9f0beaf7 745 csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
746 csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
10757ee9 747 nClusters = -1;
9f0beaf7 748 fFitterParY.ClearPoints();
749 fFitterParZ.ClearPoints();
10757ee9 750 }
751 }
752 } // for (Int_t irow = 0; irow < 159; irow++)
753 // mean chi^2 for all tracklet fits in Y and in Z direction:
754 csigmaY = TMath::Sqrt(csigmaY / nTrackletsAll);
755 csigmaZ = TMath::Sqrt(csigmaZ / nTrackletsAll);
756 // ---------------------------------------------------------------------
757
758 for (Int_t irow = 0; irow < 159; irow++){
759 // loop again over all rows along the track
760 // do analysis
761 //
762 Int_t nclFound = 0; // number of clusters in the neighborhood
763 Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster
764 Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster
765 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
766 if (!cluster0) continue;
767 Int_t sector = cluster0->GetDetector();
768 Float_t xref = cluster0->GetX();
769
770 // Make Fit
9f0beaf7 771 fFitterParY.ClearPoints();
772 fFitterParZ.ClearPoints();
773 fFitterLinY1.ClearPoints();
774 fFitterLinZ1.ClearPoints();
775 fFitterLinY2.ClearPoints();
776 fFitterLinZ2.ClearPoints();
10757ee9 777
778 // fit tracklet (clusters in given padrow +- kDelta padrows)
779 // with polynom of 2nd order and two polynoms of 1st order
780 // take both polynoms of 1st order, calculate difference of their parameters
781 // add covariance matrixes and calculate chi2 of this difference
782 // if this chi2 is bigger than a given threshold, assume that the current cluster is
783 // a kink an goto next padrow
784
785 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
786 // loop over irow +- kDelta rows (neighboured rows)
787 //
788 //
789 if (idelta == 0) continue; // don't use center cluster
790 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
791 AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
792 if (!currentCluster) continue;
793 if (currentCluster->GetType() < 0) continue;
794 if (currentCluster->GetDetector() != sector) continue;
795 Double_t x = currentCluster->GetX() - xref; // x = differece: current cluster - cluster @ irow
796 nclFound++;
797 if (idelta < 0){
798 ncl0++;
9f0beaf7 799 fFitterLinY1.AddPoint(&x, currentCluster->GetY(), csigmaY);
800 fFitterLinZ1.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
10757ee9 801 }
802 if (idelta > 0){
803 ncl1++;
9f0beaf7 804 fFitterLinY2.AddPoint(&x, currentCluster->GetY(), csigmaY);
805 fFitterLinZ2.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
10757ee9 806 }
9f0beaf7 807 fFitterParY.AddPoint(&x, currentCluster->GetY(), csigmaY);
808 fFitterParZ.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
10757ee9 809 } // loop over neighbourhood for fitter filling
9f0beaf7 810
811
10757ee9 812
813 if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
814 if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
815 if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow
9f0beaf7 816 fFitterParY.Eval();
817 fFitterParZ.Eval();
818 Double_t chi2 = (fFitterParY.GetChisquare() + fFitterParZ.GetChisquare()) / (2. * nclFound - 6.);
9f0beaf7 819 TTreeSRedirector *cstream = GetDebugStreamer();
820 if (cstream){
821 (*cstream)<<"Cut9"<<
822 "chi2="<<chi2<<
823 "\n";
824 }
10757ee9 825 // REMOVE KINK
826 // only when there are enough clusters (4) in each direction
827 if (ncl0 > 4){
9f0beaf7 828 fFitterLinY1.Eval();
829 fFitterLinZ1.Eval();
10757ee9 830 }
831 if (ncl1 > 4){
9f0beaf7 832 fFitterLinY2.Eval();
833 fFitterLinZ2.Eval();
10757ee9 834 }
835
836 if (ncl0 > 4 && ncl1 > 4){
9f0beaf7 837 fFitterLinY1.GetCovarianceMatrix(matrixY0);
838 fFitterLinY2.GetCovarianceMatrix(matrixY1);
839 fFitterLinZ1.GetCovarianceMatrix(matrixZ0);
840 fFitterLinZ2.GetCovarianceMatrix(matrixZ1);
841 fFitterLinY2.GetParameters(paramY1);
842 fFitterLinZ2.GetParameters(paramZ1);
843 fFitterLinY1.GetParameters(paramY0);
844 fFitterLinZ1.GetParameters(paramZ0);
10757ee9 845 paramY0 -= paramY1;
846 paramZ0 -= paramZ1;
847 matrixY0 += matrixY1;
848 matrixZ0 += matrixZ1;
d23650ed 849 Double_t cchi2 = 0;
10757ee9 850
851 TMatrixD difY(2, 1, paramY0.GetMatrixArray());
852 TMatrixD difYT(1, 2, paramY0.GetMatrixArray());
853 matrixY0.Invert();
854 TMatrixD mulY(matrixY0, TMatrixD::kMult, difY);
855 TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY);
d23650ed 856 cchi2 += chi2Y(0, 0);
10757ee9 857
858 TMatrixD difZ(2, 1, paramZ0.GetMatrixArray());
859 TMatrixD difZT(1, 2, paramZ0.GetMatrixArray());
860 matrixZ0.Invert();
861 TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ);
862 TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ);
d23650ed 863 cchi2 += chi2Z(0, 0);
10757ee9 864
9f0beaf7 865 if (cstream){
866 (*cstream)<<"Cut8"<<
d23650ed 867 "chi2="<<cchi2<<
9f0beaf7 868 "\n";
869 }
10757ee9 870 }
871
872 // current padrow has no kink
873
874 // get fit parameters from pol2 fit:
875 Double_t paramY[4], paramZ[4];
9f0beaf7 876 paramY[0] = fFitterParY.GetParameter(0);
877 paramY[1] = fFitterParY.GetParameter(1);
878 paramY[2] = fFitterParY.GetParameter(2);
879 paramZ[0] = fFitterParZ.GetParameter(0);
880 paramZ[1] = fFitterParZ.GetParameter(1);
881 paramZ[2] = fFitterParZ.GetParameter(2);
10757ee9 882
883 Double_t tracky = paramY[0];
884 Double_t trackz = paramZ[0];
885 Float_t deltay = tracky - cluster0->GetY();
886 Float_t deltaz = trackz - cluster0->GetZ();
887 Float_t angley = paramY[1] - paramY[0] / xref;
888 Float_t anglez = paramZ[1];
889
890 Float_t max = cluster0->GetMax();
891 UInt_t isegment = cluster0->GetDetector() % 36;
892 Int_t padSize = 0; // short pads
893 if (cluster0->GetDetector() >= 36) {
894 padSize = 1; // medium pads
895 if (cluster0->GetRow() > 63) padSize = 2; // long pads
896 }
897
898 // =========================================
899 // wirte collected information to histograms
900 // =========================================
901
902 TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector);
903 if ( irow >= kFirstLargePad) max /= kLargePadSize;
15e48021 904 Double_t smax = TMath::Sqrt(max);
905 profAmpRow->Fill( (Double_t)cluster0->GetRow(), smax );
10757ee9 906 TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector);
15e48021 907 hisAmp->Fill(smax);
10757ee9 908
909 // remove the following two lines one day:
910 TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
15e48021 911 profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
10757ee9 912
913 TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
15e48021 914 profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
10757ee9 915
916 fHclusterPerPadrow->Fill(irow); // fill histogram showing clusters per padrow
917 Int_t ipad = TMath::Nint(cluster0->GetPad());
918 Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
919 fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
920
921
922 TH3F * his3 = 0;
923 his3 = (TH3F*)fRMSY->At(padSize);
924 if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
925 his3 = (TH3F*)fRMSZ->At(padSize);
926 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
927
928 his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
929 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
930 his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
931 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
932
933
934 // Fill resolution histograms
935 Bool_t useForResol = kTRUE;
9f0beaf7 936 if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
10757ee9 937
9f0beaf7 938 if (cstream){
939 Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
940 Float_t sy = cluster0->GetSigmaY2();
941 Float_t sz = cluster0->GetSigmaZ2();
942 (*cstream)<<"Resol0"<<
108953e9 943 "run="<<fRun<< // run number
944 "event="<<fEvent<< // event number
945 "time="<<fTime<< // time stamp of event
946 "trigger="<<fTrigger<< // trigger
947 "mag="<<fMagF<< // magnetic field
9f0beaf7 948 "padSize="<<padSize<<
949 "angley="<<angley<<
950 "anglez="<<anglez<<
951 "zdr="<<zdrift<<
952 "dy="<<deltay<<
953 "dz="<<deltaz<<
954 "sy="<<sy<<
955 "sz="<<sz<<
956 "\n";
957 }
958
10757ee9 959 if (useForResol){
960 fDeltaY->Fill(deltay);
961 fDeltaZ->Fill(deltaz);
962 his3 = (TH3F*)fResolY->At(padSize);
963 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
964 his3 = (TH3F*)fResolZ->At(padSize);
965 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
966 his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
967 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
968 his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
969 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
970 }
971
972 //=============================================================================================
973
974 if (useForResol && nclFound > 2 * kMinRatio * kDelta
ae28e92e 975 && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
2acad464 976 if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
10757ee9 977 FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
978 } // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
af6a50bb 979 //
980 // Fill THN histograms
981 //
982 Double_t xvar[9];
983 xvar[1]=padSize;
984 xvar[2]=1.-TMath::Abs(cluster0->GetZ())/250.;
985 xvar[3]=cluster0->GetMax();
986 xvar[5]=angley;
987 xvar[6]=anglez;
988
989 xvar[4]=cluster0->GetPad()-TMath::Nint(cluster0->GetPad());
990 xvar[0]=deltay;
991 fHisDeltaY->Fill(xvar);
992 xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
993 fHisRMSY->Fill(xvar);
994
995 xvar[4]=cluster0->GetTimeBin()-TMath::Nint(cluster0->GetTimeBin());
996 xvar[0]=deltaz;
997 fHisDeltaZ->Fill(xvar);
998 xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
999 fHisRMSZ->Fill(xvar);
10757ee9 1000
1001 } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
1002} // FillResolutionHistoLocal(...)
1003
1004
1005
1006void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t angley, Float_t anglez, Int_t nclFound, Int_t kDelta) {
1007 //
1008 // - debug part of FillResolutionHistoLocal -
1009 // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
ae28e92e 1010 // called only for GetStreamLevel() > 4
10757ee9 1011 // fill resolution trees
1012 //
1013
1014 Int_t sector = cluster0->GetDetector();
1015 Float_t xref = cluster0->GetX();
1016 Int_t padSize = 0; // short pads
1017 if (cluster0->GetDetector() >= 36) {
1018 padSize = 1; // medium pads
1019 if (cluster0->GetRow() > 63) padSize = 2; // long pads
1020 }
1021
9f0beaf7 1022 static TLinearFitter fitY0(3, "pol2");
1023 static TLinearFitter fitZ0(3, "pol2");
1024 static TLinearFitter fitY2(5, "hyp4");
1025 static TLinearFitter fitZ2(5, "hyp4");
1026 static TLinearFitter fitY2Q(5, "hyp4");
1027 static TLinearFitter fitZ2Q(5, "hyp4");
1028 static TLinearFitter fitY2S(5, "hyp4");
1029 static TLinearFitter fitZ2S(5, "hyp4");
1030 fitY0.ClearPoints();
1031 fitZ0.ClearPoints();
1032 fitY2.ClearPoints();
1033 fitZ2.ClearPoints();
1034 fitY2Q.ClearPoints();
1035 fitZ2Q.ClearPoints();
1036 fitY2S.ClearPoints();
1037 fitZ2S.ClearPoints();
1038
1039 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
1040 // loop over irow +- kDelta rows (neighboured rows)
1041 //
1042 //
1043 if (idelta == 0) continue;
1044 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
1045 AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
1046 if (!cluster) continue;
1047 if (cluster->GetType() < 0) continue;
1048 if (cluster->GetDetector() != sector) continue;
1049 Double_t x = cluster->GetX() - xref;
1050 Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
1051 Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
1052 //
1053 Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1054 Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1055 Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
10757ee9 1056 TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
9f0beaf7 1057 Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
10757ee9 1058 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
9f0beaf7 1059 Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1060 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1061 TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
1062 Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1063 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1064 TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
1065 sigmaYS = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
1066 sigmaZS = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
1067 //
1068 if (kDelta != 0){
1069 fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
1070 fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
1071 }
1072 Double_t xxx[4];
1073 xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
1074 xxx[1] = x;
1075 xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
1076 xxx[3] = x * x;
1077 fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
1078 fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
1079 fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
1080 fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
1081 fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
1082 fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
1083 //
1084 } // neigbouhood-loop
10757ee9 1085 //
9f0beaf7 1086 fitY0.Eval();
1087 fitZ0.Eval();
1088 fitY2.Eval();
1089 fitZ2.Eval();
1090 fitY2Q.Eval();
1091 fitZ2Q.Eval();
1092 fitY2S.Eval();
1093 fitZ2S.Eval();
1094 Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
1095 Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
1096 Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
1097 Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
1098 Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
1099 Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
1100 Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
1101 Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
1102 //
1103 static TVectorD parY0(3);
1104 static TMatrixD matY0(3, 3);
1105 static TVectorD parZ0(3);
1106 static TMatrixD matZ0(3, 3);
1107 fitY0.GetParameters(parY0);
1108 fitY0.GetCovarianceMatrix(matY0);
1109 fitZ0.GetParameters(parZ0);
1110 fitZ0.GetCovarianceMatrix(matZ0);
1111 //
1112 static TVectorD parY2(5);
1113 static TMatrixD matY2(5,5);
1114 static TVectorD parZ2(5);
1115 static TMatrixD matZ2(5,5);
1116 fitY2.GetParameters(parY2);
1117 fitY2.GetCovarianceMatrix(matY2);
1118 fitZ2.GetParameters(parZ2);
1119 fitZ2.GetCovarianceMatrix(matZ2);
1120 //
1121 static TVectorD parY2Q(5);
1122 static TMatrixD matY2Q(5,5);
1123 static TVectorD parZ2Q(5);
1124 static TMatrixD matZ2Q(5,5);
1125 fitY2Q.GetParameters(parY2Q);
1126 fitY2Q.GetCovarianceMatrix(matY2Q);
1127 fitZ2Q.GetParameters(parZ2Q);
1128 fitZ2Q.GetCovarianceMatrix(matZ2Q);
1129 static TVectorD parY2S(5);
1130 static TMatrixD matY2S(5,5);
1131 static TVectorD parZ2S(5);
1132 static TMatrixD matZ2S(5,5);
1133 fitY2S.GetParameters(parY2S);
1134 fitY2S.GetCovarianceMatrix(matY2S);
1135 fitZ2S.GetParameters(parZ2S);
1136 fitZ2S.GetCovarianceMatrix(matZ2S);
1137 Float_t sigmaY0 = TMath::Sqrt(matY0(0,0));
1138 Float_t sigmaZ0 = TMath::Sqrt(matZ0(0,0));
1139 Float_t sigmaDY0 = TMath::Sqrt(matY0(1,1));
1140 Float_t sigmaDZ0 = TMath::Sqrt(matZ0(1,1));
1141 Float_t sigmaY2 = TMath::Sqrt(matY2(1,1));
1142 Float_t sigmaZ2 = TMath::Sqrt(matZ2(1,1));
1143 Float_t sigmaDY2 = TMath::Sqrt(matY2(3,3));
1144 Float_t sigmaDZ2 = TMath::Sqrt(matZ2(3,3));
1145 Float_t sigmaY2Q = TMath::Sqrt(matY2Q(1,1));
1146 Float_t sigmaZ2Q = TMath::Sqrt(matZ2Q(1,1));
1147 Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
1148 Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
1149 Float_t sigmaY2S = TMath::Sqrt(matY2S(1,1));
1150 Float_t sigmaZ2S = TMath::Sqrt(matZ2S(1,1));
1151 Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
1152 Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
1153
1154 // Error parameters
1155 Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
1156 Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
1157 //
1158 Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1159 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1160 Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1161 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1162 Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1163 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1164 Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1165 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1166
1167
1168 // RMS parameters
1169 Float_t meanRMSY = 0;
1170 Float_t meanRMSZ = 0;
1171 Int_t nclRMS = 0;
1172 for (Int_t idelta = -2; idelta <= 2; idelta++){
1173 // loop over neighbourhood
1174 if (idelta+irow < 0 || idelta+irow > 159) continue;
1175 // if (idelta+irow>159) continue;
1176 AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
1177 if (!cluster) continue;
1178 meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
1179 meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
1180 nclRMS++;
1181 }
1182 meanRMSY /= nclRMS;
1183 meanRMSZ /= nclRMS;
1184
1185 Float_t rmsY = TMath::Sqrt(cluster0->GetSigmaY2());
1186 Float_t rmsZ = TMath::Sqrt(cluster0->GetSigmaZ2());
1187 Float_t rmsYT = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1188 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1189 Float_t rmsZT = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1190 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1191 Float_t rmsYT0 = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1192 TMath::Abs(angley));
1193 Float_t rmsZT0 = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
10757ee9 1194 TMath::Abs(anglez));
9f0beaf7 1195 Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1196 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1197 Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1198 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1199 Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1200 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1201 rmsY,meanRMSY);
1202 Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1203 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1204 rmsZ,meanRMSZ);
1205 //
1206 // cluster debug
1207 TTreeSRedirector *cstream = GetDebugStreamer();
1208 if (cstream){
1209 (*cstream)<<"ResolCl"<< // valgrind 3 40,000 bytes in 1 blocks are possibly
108953e9 1210 "run="<<fRun<< // run number
1211 "event="<<fEvent<< // event number
1212 "time="<<fTime<< // time stamp of event
1213 "trigger="<<fTrigger<< // trigger
1214 "mag="<<fMagF<< // magnetic field
9f0beaf7 1215 "Sector="<<sector<<
1216 "Cl.="<<cluster0<<
1217 "CSigmaY0="<<csigmaY0<< // cluster errorY
1218 "CSigmaYQ="<<csigmaYQ<< // cluster errorY - q scaled
1219 "CSigmaYS="<<csigmaYS<< // cluster errorY - q scaled
1220 "CSigmaZ0="<<csigmaZ0<< //
1221 "CSigmaZQ="<<csigmaZQ<<
1222 "CSigmaZS="<<csigmaZS<<
1223 "shapeYF="<<rmsYFactor<<
1224 "shapeZF="<<rmsZFactor<<
1225 "rmsY="<<rmsY<<
1226 "rmsZ="<<rmsZ<<
1227 "rmsYM="<<meanRMSY<<
1228 "rmsZM="<<meanRMSZ<<
1229 "rmsYT="<<rmsYT<<
1230 "rmsZT="<<rmsZT<<
1231 "rmsYT0="<<rmsYT0<<
1232 "rmsZT0="<<rmsZT0<<
1233 "rmsYS="<<rmsYSigma<<
1234 "rmsZS="<<rmsZSigma<<
1235 "padSize="<<padSize<<
1236 "Ncl="<<nclFound<<
1237 "PY0.="<<&parY0<<
1238 "PZ0.="<<&parZ0<<
1239 "SigmaY0="<<sigmaY0<<
1240 "SigmaZ0="<<sigmaZ0<<
1241 "angley="<<angley<<
1242 "anglez="<<anglez<<
1243 "\n";
1244
1245 // tracklet dubug
1246 (*cstream)<<"ResolTr"<<
108953e9 1247 "run="<<fRun<< // run number
1248 "event="<<fEvent<< // event number
1249 "time="<<fTime<< // time stamp of event
1250 "trigger="<<fTrigger<< // trigger
1251 "mag="<<fMagF<< // magnetic field
9f0beaf7 1252 "padSize="<<padSize<<
1253 "IPad="<<padSize<<
1254 "Sector="<<sector<<
1255 "Ncl="<<nclFound<<
1256 "chi2Y0="<<chi2Y0<<
1257 "chi2Z0="<<chi2Z0<<
1258 "chi2Y2="<<chi2Y2<<
1259 "chi2Z2="<<chi2Z2<<
1260 "chi2Y2Q="<<chi2Y2Q<<
1261 "chi2Z2Q="<<chi2Z2Q<<
1262 "chi2Y2S="<<chi2Y2S<<
1263 "chi2Z2S="<<chi2Z2S<<
1264 "PY0.="<<&parY0<<
1265 "PZ0.="<<&parZ0<<
1266 "PY2.="<<&parY2<<
1267 "PZ2.="<<&parZ2<<
1268 "PY2Q.="<<&parY2Q<<
1269 "PZ2Q.="<<&parZ2Q<<
1270 "PY2S.="<<&parY2S<<
1271 "PZ2S.="<<&parZ2S<<
1272 "SigmaY0="<<sigmaY0<<
1273 "SigmaZ0="<<sigmaZ0<<
1274 "SigmaDY0="<<sigmaDY0<<
1275 "SigmaDZ0="<<sigmaDZ0<<
1276 "SigmaY2="<<sigmaY2<<
1277 "SigmaZ2="<<sigmaZ2<<
1278 "SigmaDY2="<<sigmaDY2<<
1279 "SigmaDZ2="<<sigmaDZ2<<
1280 "SigmaY2Q="<<sigmaY2Q<<
1281 "SigmaZ2Q="<<sigmaZ2Q<<
1282 "SigmaDY2Q="<<sigmaDY2Q<<
1283 "SigmaDZ2Q="<<sigmaDZ2Q<<
1284 "SigmaY2S="<<sigmaY2S<<
1285 "SigmaZ2S="<<sigmaZ2S<<
1286 "SigmaDY2S="<<sigmaDY2S<<
1287 "SigmaDZ2S="<<sigmaDZ2S<<
1288 "angley="<<angley<<
1289 "anglez="<<anglez<<
1290 "\n";
1291 }
10757ee9 1292}
1293
1294
1295
1296
1297
1298TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
1299 //
1300 // creates a new histogram which contains the difference between
1301 // the histogram hfit and the function func
1302 //
1303 TH2D * result = (TH2D*)hfit->Clone(); // valgrind 3 40,139 bytes in 11 blocks are still reachable
1304 result->SetTitle(Form("%s fit residuals",result->GetTitle()));
1305 result->SetName(Form("%s fit residuals",result->GetName()));
1306 TAxis *xaxis = hfit->GetXaxis();
1307 TAxis *yaxis = hfit->GetYaxis();
1308 Double_t x[2];
1309 for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
1310 x[1] = yaxis->GetBinCenter(biny);
1311 for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
1312 x[0] = xaxis->GetBinCenter(binx);
1313 Int_t bin = hfit->GetBin(binx, biny);
1314 Double_t val = hfit->GetBinContent(bin);
1315// result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
1316 result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
1317 }
1318 }
1319 return result;
1320}
1321
1322
1323void AliTPCcalibTracks::SetStyle() const {
1324 //
1325 // set style, can be called by all draw functions
1326 //
1327 gROOT->SetStyle("Plain");
1328 gStyle->SetFillColor(10);
1329 gStyle->SetPadColor(10);
1330 gStyle->SetCanvasColor(10);
1331 gStyle->SetStatColor(10);
1332 gStyle->SetPalette(1,0);
1333 gStyle->SetNumberContours(60);
1334}
1335
1336
1337void AliTPCcalibTracks::Draw(Option_t* opt){
1338 //
1339 // draw-function of AliTPCcalibTracks
1340 // will draws some exemplaric pictures
1341 //
1342
ae28e92e 1343 if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
10757ee9 1344 SetStyle();
1345 Double_t min = 0;
1346 Double_t max = 0;
1347 TCanvas *c1 = new TCanvas();
1348 c1->Divide(0, 3);
1349 TVirtualPad *upperThird = c1->GetPad(1);
1350 TVirtualPad *middleThird = c1->GetPad(2);
1351 TVirtualPad *lowerThird = c1->GetPad(3);
1352 upperThird->Divide(2,0);
1353 TVirtualPad *upleft = upperThird->GetPad(1);
1354 TVirtualPad *upright = upperThird->GetPad(2);
1355 middleThird->Divide(2,0);
1356 TVirtualPad *middleLeft = middleThird->GetPad(1);
1357 TVirtualPad *middleRight = middleThird->GetPad(2);
1358 lowerThird->Divide(2,0);
1359 TVirtualPad *downLeft = lowerThird->GetPad(1);
1360 TVirtualPad *downRight = lowerThird->GetPad(2);
1361
1362
1363 upleft->cd(0);
1364 min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1365 max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1366 fDeltaY->SetAxisRange(min, max);
1367 fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
1368 c1->Update();
1369
1370 upright->cd(0);
1371 max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1372 min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1373 fDeltaZ->SetAxisRange(min, max);
1374 fDeltaZ->Fit("gaus","q","",min, max);
1375 c1->Update();
1376
1377 middleLeft->cd();
1378 fHclus->Draw(opt);
1379
1380 middleRight->cd();
1381 fRejectedTracksHisto->Draw(opt);
1382 TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
1383 TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
1384 TText *t2 = pt->AddText("2: kMinClusters");
1385 TText *t3 = pt->AddText("3: kMinRatio");
1386 TText *t4 = pt->AddText("4: kMax1pt");
1387 t1 = t1; t2 = t2; t3 = t3; t4 = t4; // avoid compiler warnings
1388 pt->SetToolTipText("Legend for failed cuts");
1389 pt->Draw();
1390
1391 downLeft->cd();
1392 fHclusterPerPadrowRaw->Draw(opt);
1393
1394 downRight->cd();
1395 fHclusterPerPadrow->Draw(opt);
1396}
1397
1398
f78da5ae 1399void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
10757ee9 1400 //
1401 // all functions are called, that produce pictures
1402 // the histograms are written to the directory 'pathName'
1403 // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
1404 // 'stat' is also the number of minEntries for MakeResPlotsQTree
1405 //
1406
ae28e92e 1407 if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
10757ee9 1408 MakeAmpPlots(stat, pathName);
1409 MakeDeltaPlots(pathName);
1410 FitResolutionNew(pathName);
1411 FitRMSNew(pathName);
1412 MakeChargeVsDriftLengthPlots(pathName);
1413// MakeResPlotsQ(1, 1);
1414 MakeResPlotsQTree(stat, pathName);
1415}
1416
1417
f78da5ae 1418void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, const char* pathName){
10757ee9 1419 //
1420 // creates several plots:
1421 // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
1422 // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
1423 // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
1424 // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1425 // Empty histograms (sectors without data) are not written to file
1426 // the ps-files are written to the directory 'pathName', that is created if it does not exist
1427 // 'stat': only histograms with more than 'stat' entries are written to file.
1428 //
1429
1430 SetStyle();
1431 gSystem->MakeDirectory(pathName);
1432 gSystem->ChangeDirectory(pathName);
1433
1434 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1435 TPostScript *ps;
1436 // histograms with accumulated amplitude for all IROCs and OROCs
1437 TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
1438 allAmpHisIROC->SetName("Amp all IROCs");
1439 allAmpHisIROC->SetTitle("Amp all IROCs");
1440 TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
1441 allAmpHisOROC->SetName("Amp all OROCs");
1442 allAmpHisOROC->SetTitle("Amp all OROCs");
1443
1444 ps = new TPostScript("fArrayAmp.ps", 112);
ae28e92e 1445 if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
10757ee9 1446 for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
1447 if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat ) continue;
1448 ps->NewPage();
1449 ((TH1F*)fArrayAmp->At(i))->Draw();
1450 c1->Update(); // valgrind 3
1451 if (i > 0 && i < 36) {
1452 allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
1453 allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
1454 }
1455 }
1456 ps->NewPage();
1457 allAmpHisIROC->Draw();
1458 c1->Update(); // valgrind
1459 ps->NewPage();
1460 allAmpHisOROC->Draw();
1461 c1->Update();
1462 ps->Close();
1463 delete ps;
1464
1465 TH1F *his = 0;
1466 Double_t min = 0;
1467 Double_t max = 0;
1468 ps = new TPostScript("fArrayAmpRow.ps", 112);
ae28e92e 1469 if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
10757ee9 1470 for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
1471 his = (TH1F*)fArrayAmpRow->At(i);
1472 if (his->GetEntries() < stat) continue;
1473 ps->NewPage();
1474 min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
1475 max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
1476 his->SetAxisRange(min, max);
1477 his->Fit("pol3", "q", "", min, max);
1478 // his->Draw("error"); // don't use this line when you don't want to have empty pages in the ps-file
1479 c1->Update();
1480 }
1481 ps->Close();
1482 delete ps;
1483 delete c1;
1484 gSystem->ChangeDirectory("..");
1485}
1486
1487
f78da5ae 1488void AliTPCcalibTracks::MakeDeltaPlots(const char* pathName){
10757ee9 1489 //
1490 // creates several plots:
1491 // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1492 // the ps-files are written to the directory 'pathName', that is created if it does not exist
1493 //
1494
1495 SetStyle();
1496 gSystem->MakeDirectory(pathName);
1497 gSystem->ChangeDirectory(pathName);
1498
1499 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1500 TPostScript *ps;
1501 Double_t min = 0;
1502 Double_t max = 0;
1503
1504 ps = new TPostScript("DeltaYZ.ps", 112);
ae28e92e 1505 if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
10757ee9 1506 min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1507 max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1508 fDeltaY->SetAxisRange(min, max);
1509 ps->NewPage();
1510 fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
1511 c1->Update();
1512 ps->NewPage();
1513 max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1514 min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1515 fDeltaZ->SetAxisRange(min, max);
1516 fDeltaZ->Fit("gaus","q","",min, max);
1517 c1->Update();
1518 ps->Close();
1519 delete ps;
1520 delete c1;
1521 gSystem->ChangeDirectory("..");
1522}
1523
1524
f78da5ae 1525void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(const char* pathName){
10757ee9 1526 //
1527 // creates charge vs. driftlength plots, one TProfile for each ROC
1528 // is not correct like this, should be one TProfile for each sector and padsize
1529 //
1530
1531 SetStyle();
1532 gSystem->MakeDirectory(pathName);
1533 gSystem->ChangeDirectory(pathName);
1534
1535 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1536 TPostScript *ps;
1537 ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
ae28e92e 1538 if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
10757ee9 1539 TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
1540 TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
1541 chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
1542 chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
1543 chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
1544 chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
1545
1546 for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
1547 ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
1548 c1->Update();
1549 if (i > 0 && i < 36) {
1550 chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
1551 chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
1552 }
1553 ps->NewPage();
1554 }
1555 chargeVsDriftlengthAllIROCs->Draw();
1556 c1->Update(); // valgrind
1557 ps->NewPage();
1558 chargeVsDriftlengthAllOROCs->Draw();
1559 c1->Update();
1560 ps->Close();
1561 delete ps;
1562 delete c1;
1563 gSystem->ChangeDirectory("..");
1564}
1565
1566
f78da5ae 1567void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(const char* pathName){
10757ee9 1568 //
1569 // creates charge vs. driftlength plots, one TProfile for each ROC
1570 // under development....
1571 //
1572
1573 SetStyle();
1574 gSystem->MakeDirectory(pathName);
1575 gSystem->ChangeDirectory(pathName);
1576
1577 TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700)); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1578// TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
1579 c1->Divide(0,3);
1580 TPostScript *ps;
1581 ps = new TPostScript("chargeVsDriftlength.ps", 111);
ae28e92e 1582 if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
10757ee9 1583
1584 TProfile *chargeVsDriftlengthAllShortPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
1585 TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
1586 TProfile *chargeVsDriftlengthAllLongPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
1587 chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
1588 chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
1589 chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
1590 chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
1591 chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
1592 chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
1593
1594 for (Int_t i = 0; i < 36; i++) {
1595 c1->cd(1)->SetGridx();
1596 c1->cd(1)->SetGridy();
1597 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
1598 c1->cd(2)->SetGridx();
1599 c1->cd(2)->SetGridy();
1600 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
1601 c1->cd(3)->SetGridx();
1602 c1->cd(3)->SetGridy();
1603 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
1604 c1->Update();
1605 chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
1606 chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
1607 chargeVsDriftlengthAllLongPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
1608 ps->NewPage();
1609 }
1610 c1->cd(1)->SetGridx();
1611 c1->cd(1)->SetGridy();
1612 chargeVsDriftlengthAllShortPads->Draw();
1613 c1->cd(2)->SetGridx();
1614 c1->cd(2)->SetGridy();
1615 chargeVsDriftlengthAllMediumPads->Draw();
1616 c1->cd(3)->SetGridx();
1617 c1->cd(3)->SetGridy();
1618 chargeVsDriftlengthAllLongPads->Draw();
1619 c1->Update(); // valgrind
1620// ps->NewPage();
1621 ps->Close();
1622 delete ps;
1623 delete c1;
1624 gSystem->ChangeDirectory("..");
1625}
1626
1627
1628
f78da5ae 1629void AliTPCcalibTracks::FitResolutionNew(const char* pathName){
10757ee9 1630 //
1631 // calculates different resulution fits in Y and Z direction
1632 // the histograms are written to 'ResolutionYZ.ps'
1633 // writes calculated resolution to 'resol.txt'
1634 // all files are stored in the directory pathName
1635 //
1636
1637 SetStyle();
1638 gSystem->MakeDirectory(pathName);
1639 gSystem->ChangeDirectory(pathName);
1640
1641 TCanvas c;
1642 c.Divide(2,1);
ae28e92e 1643 if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
10757ee9 1644 TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112);
1645 TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1646 fres->SetParameter(0,0.02);
1647 fres->SetParameter(1,0.0054);
1648 fres->SetParameter(2,0.13);
1649
1650 TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
1651
1652 // create histogramw for Y-resolution
1653 TH3F * hisResY0 = (TH3F*)fResolY->At(0);
1654 hisResY0->FitSlicesZ();
1655 TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
1656 TH3F * hisResY1 = (TH3F*)fResolY->At(1);
1657 hisResY1->FitSlicesZ();
1658 TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
1659 TH3F * hisResY2 = (TH3F*)fResolY->At(2);
1660 hisResY2->FitSlicesZ();
1661 TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
1662 //
1663 ps->NewPage();
1664 c.cd(1);
1665 hisResY02->Fit(fres, "q"); // valgrind 132,072 bytes in 6 blocks are indirectly lost
1666 hisResY02->Draw("surf1");
1667 c.cd(2);
1668 MakeDiff(hisResY02,fres)->Draw("surf1");
1669 c.Update();
1670 // c.SaveAs("ResolutionYPad0.eps");
1671 ps->NewPage();
1672 c.cd(1);
1673 hisResY12->Fit(fres, "q");
1674 hisResY12->Draw("surf1");
1675 c.cd(2);
1676 MakeDiff(hisResY12,fres)->Draw("surf1");
1677 c.Update();
1678 // c.SaveAs("ResolutionYPad1.eps");
1679 ps->NewPage();
1680 c.cd(1);
1681 hisResY22->Fit(fres, "q");
1682 hisResY22->Draw("surf1");
1683 c.cd(2);
1684 MakeDiff(hisResY22,fres)->Draw("surf1");
1685 c.Update();
1686// c.SaveAs("ResolutionYPad2.eps");
1687
1688 // create histogramw for Z-resolution
1689 TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
1690 hisResZ0->FitSlicesZ();
1691 TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
1692 TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
1693 hisResZ1->FitSlicesZ();
1694 TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
1695 TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
1696 hisResZ2->FitSlicesZ();
1697 TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
1698
1699 ps->NewPage();
1700 c.cd(1);
1701 hisResZ02->Fit(fres, "q");
1702 hisResZ02->Draw("surf1");
1703 c.cd(2);
1704 MakeDiff(hisResZ02,fres)->Draw("surf1");
1705 c.Update();
1706// c.SaveAs("ResolutionZPad0.eps");
1707 ps->NewPage();
1708 c.cd(1);
1709 hisResZ12->Fit(fres, "q");
1710 hisResZ12->Draw("surf1");
1711 c.cd(2);
1712 MakeDiff(hisResZ12,fres)->Draw("surf1");
1713 c.Update();
1714// c.SaveAs("ResolutionZPad1.eps");
1715 ps->NewPage();
1716 c.cd(1);
1717 hisResZ22->Fit(fres, "q");
1718 hisResZ22->Draw("surf1");
1719 c.cd(2);
1720 MakeDiff(hisResZ22,fres)->Draw("surf1");
1721 c.Update();
1722// c.SaveAs("ResolutionZPad2.eps");
1723 ps->Close();
1724 delete ps;
1725
1726 // write calculated resoltuions to 'resol.txt'
1727 ofstream fresol("resol.txt");
1728 fresol<<"Pad 0.75 cm"<<"\n";
1729 hisResY02->Fit(fres, "q"); // valgrind
1730 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1731 hisResZ02->Fit(fres, "q");
1732 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1733 //
1734 fresol<<"Pad 1.00 cm"<<1<<"\n";
1735 hisResY12->Fit(fres, "q"); // valgrind
1736 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1737 hisResZ12->Fit(fres, "q");
1738 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1739 //
1740 fresol<<"Pad 1.50 cm"<<0<<"\n";
1741 hisResY22->Fit(fres, "q");
1742 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1743 hisResZ22->Fit(fres, "q");
1744 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1745
1746 TH1::AddDirectory(kFALSE);
1747 gSystem->ChangeDirectory("..");
1748 delete fres;
1749}
1750
1751
f78da5ae 1752void AliTPCcalibTracks::FitRMSNew(const char* pathName){
10757ee9 1753 //
1754 // calculates different resulution-rms fits in Y and Z direction
1755 // the histograms are written to 'RMS_YZ.ps'
1756 // writes calculated resolution rms to 'rms.txt'
1757 // all files are stored in the directory pathName
1758 //
1759
1760 SetStyle();
1761 gSystem->MakeDirectory(pathName);
1762 gSystem->ChangeDirectory(pathName);
1763
1764 TCanvas c; // valgrind 3 42,120 bytes in 405 blocks are still reachable 23,816 bytes in 229 blocks are still reachable
1765 c.Divide(2,1);
ae28e92e 1766 if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
10757ee9 1767 TPostScript *ps = new TPostScript("RMS_YZ.ps", 112);
1768 TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1769 frms->SetParameter(0,0.02);
1770 frms->SetParameter(1,0.0054);
1771 frms->SetParameter(2,0.13);
1772
1773 TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
1774
1775 // create histogramw for Y-RMS
1776 TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
1777 hisResY0->FitSlicesZ();
1778 TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
1779 TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
1780 hisResY1->FitSlicesZ();
1781 TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
1782 TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
1783 hisResY2->FitSlicesZ();
1784 TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
1785 //
1786 ps->NewPage();
1787 c.cd(1);
1788 hisResY02->Fit(frms, "qn0");
1789 hisResY02->Draw("surf1");
1790 c.cd(2);
1791 MakeDiff(hisResY02,frms)->Draw("surf1");
1792 c.Update();
1793// c.SaveAs("RMSYPad0.eps");
1794 ps->NewPage();
1795 c.cd(1);
1796 hisResY12->Fit(frms, "qn0"); // valgrind several blocks possibly lost
1797 hisResY12->Draw("surf1");
1798 c.cd(2);
1799 MakeDiff(hisResY12,frms)->Draw("surf1");
1800 c.Update();
1801// c.SaveAs("RMSYPad1.eps");
1802 ps->NewPage();
1803 c.cd(1);
1804 hisResY22->Fit(frms, "qn0");
1805 hisResY22->Draw("surf1");
1806 c.cd(2);
1807 MakeDiff(hisResY22,frms)->Draw("surf1");
1808 c.Update();
1809// c.SaveAs("RMSYPad2.eps");
1810
1811 // create histogramw for Z-RMS
1812 TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
1813 hisResZ0->FitSlicesZ();
1814 TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
1815 TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1);
1816 hisResZ1->FitSlicesZ();
1817 TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
1818 TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2);
1819 hisResZ2->FitSlicesZ();
1820 TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
1821 //
1822 ps->NewPage();
1823 c.cd(1);
1824 hisResZ02->Fit(frms, "qn0"); // valgrind
1825 hisResZ02->Draw("surf1");
1826 c.cd(2);
1827 MakeDiff(hisResZ02,frms)->Draw("surf1");
1828 c.Update();
1829// c.SaveAs("RMSZPad0.eps");
1830 ps->NewPage();
1831 c.cd(1);
1832 hisResZ12->Fit(frms, "qn0");
1833 hisResZ12->Draw("surf1");
1834 c.cd(2);
1835 MakeDiff(hisResZ12,frms)->Draw("surf1");
1836 c.Update();
1837// c.SaveAs("RMSZPad1.eps");
1838 ps->NewPage();
1839 c.cd(1);
1840 hisResZ22->Fit(frms, "qn0"); // valgrind 1 block possibly lost
1841 hisResZ22->Draw("surf1");
1842 c.cd(2);
1843 MakeDiff(hisResZ22,frms)->Draw("surf1");
1844 c.Update();
1845// c.SaveAs("RMSZPad2.eps");
1846
1847 // write calculated resoltuion rms to 'rms.txt'
1848 ofstream filerms("rms.txt");
1849 filerms<<"Pad 0.75 cm"<<"\n";
1850 hisResY02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
1851 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1852 hisResZ02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
1853 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1854 //
1855 filerms<<"Pad 1.00 cm"<<1<<"\n";
1856 hisResY12->Fit(frms, "qn0"); // valgrind 3,256 bytes in 22 blocks are indirectly lost
1857 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1858 hisResZ12->Fit(frms, "qn0"); // valgrind 66,036 bytes in 3 blocks are still reachable
1859 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1860 //
1861 filerms<<"Pad 1.50 cm"<<0<<"\n";
1862 hisResY22->Fit(frms, "qn0"); // valgrind 40,139 bytes in 11 blocks are still reachable 330,180 bytes in 15 blocks are possibly lost
1863 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1864 hisResZ22->Fit(frms, "qn0");
1865 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1866
1867 TH1::AddDirectory(kFALSE);
1868 gSystem->ChangeDirectory("..");
1869 ps->Close();
1870 delete ps;
1871 delete frms;
1872}
1873
1874
f78da5ae 1875void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
10757ee9 1876 //
1877 // Make tree with resolution parameters
1878 // the result is written to 'resol.root' in directory 'pathname'
1879 // file information are available in fileInfo
1880 // available variables in the tree 'Resol':
1881 // Entries: number of entries for this resolution point
1882 // nbins: number of bins that were accumulated
1883 // Dim: direction, Dim==0: y-direction, Dim==1: z-direction
1884 // Pad: padSize; short, medium and long
1885 // Length: pad length, 0.75, 1, 1.5
1886 // QMean: mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1887 // Zc: center of middle bin in drift direction
1888 // Zm: mean dirftlength for accumulated Delta-Histograms
1889 // Zs: width of driftlength bin
1890 // AngleC: center of middle bin in Angle-Direction
1891 // AngleM: mean angle for accumulated Delta-Histograms
1892 // AngleS: width of Angle-bin
1893 // Resol: sigma for gaus fit through Delta-Histograms
1894 // Sigma: error of sigma for gaus fit through Delta Histograms
1895 // MeanR: mean of the Delta-Histogram
1896 // SigmaR: rms of the Delta-Histogram
1897 // RMSm: mean of the gaus fit through RMS-Histogram
1898 // RMS: sigma of the gaus fit through RMS-Histogram
1899 // RMSe0: error of mean of gaus fit in RMS-Histogram
1900 // RMSe1: error of sigma of gaus fit in RMS-Histogram
1901 //
1902
ae28e92e 1903 if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
1904 if (GetDebugLevel() > -1) cout << " relax, the calculation will take a while..." << endl;
10757ee9 1905
1906 gSystem->MakeDirectory(pathName);
1907 gSystem->ChangeDirectory(pathName);
1908 TString kFileName = "resol.root";
1909 TTreeSRedirector fTreeResol(kFileName.Data());
1910
1911 TH3F *resArray[2][3][11];
1912 TH3F *rmsArray[2][3][11];
1913
1914 // load histograms from fArrayQDY and fArrayQDZ
1915 // into resArray and rmsArray
1916 // that is all we need here
1917 for (Int_t idim = 0; idim < 2; idim++){
1918 for (Int_t ipad = 0; ipad < 3; ipad++){
1919 for (Int_t iq = 0; iq <= 10; iq++){
1920 rmsArray[idim][ipad][iq]=0;
1921 resArray[idim][ipad][iq]=0;
1922 Int_t bin = GetBin(iq,ipad);
1923 TH3F *hresl = 0;
1924 if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
1925 if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
1926 if (!hresl) continue;
1927 resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
1928 resArray[idim][ipad][iq]->SetDirectory(0);
1929 TH3F * hreslRMS = 0;
1930 if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
1931 if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
1932 if (!hreslRMS) continue;
1933 rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
1934 rmsArray[idim][ipad][iq]->SetDirectory(0);
1935 }
1936 }
1937 }
ae28e92e 1938 if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
10757ee9 1939
1940 //--------------------------------------------------------------------------------------------
1941
1942 char name[200];
1943 Double_t qMean;
1944 Double_t zMean, angleMean, zCenter, angleCenter;
1945 Double_t zSigma, angleSigma;
1946 TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1947 TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1948 TF1 *fitFunction = new TF1("fitFunction", "gaus");
1949 Float_t entriesQ = 0;
1950 Int_t loopCounter = 1;
1951
1952 for (Int_t idim = 0; idim < 2; idim++){
1953 // Loop y-z corrdinate
1954 for (Int_t ipad = 0; ipad < 3; ipad++){
1955 // loop pad type
1956 for (Int_t iq = -1; iq < 10; iq++){
1957 // LOOP Q
ae28e92e 1958 if (GetDebugLevel() > -1)
10757ee9 1959 cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
1960 << (Int_t)((loopCounter)/66.*100) << "% done), "
1961 << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << " \r" << std::flush;
1962 loopCounter++;
1963 TH3F *hres = 0;
1964 TH3F *hrms = 0;
1965 qMean = 0;
1966 entriesQ = 0;
1967
1968 // calculate qMean
1969 if (iq == -1){
1970 // integrated spectra
1971 for (Int_t iql = 0; iql < 10; iql++){
1972 Int_t bin = GetBin(iql,ipad);
1973 TH3F *hresl = resArray[idim][ipad][iql];
1974 TH3F *hrmsl = rmsArray[idim][ipad][iql];
1975 if (!hresl) continue;
1976 if (!hrmsl) continue;
1977 entriesQ += hresl->GetEntries();
1978 qMean += hresl->GetEntries() * GetQ(bin);
1979 if (!hres) {
1980 hres = (TH3F*)hresl->Clone();
1981 hrms = (TH3F*)hrmsl->Clone();
1982 }
1983 else{
1984 hres->Add(hresl);
1985 hrms->Add(hrmsl);
1986 }
1987 }
1988 qMean /= entriesQ;
1989 qMean *= -1.; // integral mean charge
1990 }
1991 else {
1992 // loop over neighboured Q-bins
1993 // accumulate entries from neighboured Q-bins
1994 for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
1995 if (iql < 0) continue;
1996 Int_t bin = GetBin(iql,ipad);
1997 TH3F * hresl = resArray[idim][ipad][iql];
1998 TH3F * hrmsl = rmsArray[idim][ipad][iql];
1999 if (!hresl) continue;
2000 if (!hrmsl) continue;
2001 entriesQ += hresl->GetEntries();
2002 qMean += hresl->GetEntries() * GetQ(bin);
2003 if (!hres) {
2004 hres = (TH3F*) hresl->Clone();
2005 hrms = (TH3F*) hrmsl->Clone();
2006 }
2007 else{
2008 hres->Add(hresl);
2009 hrms->Add(hrmsl);
2010 }
2011 }
2012 qMean/=entriesQ;
2013 }
2014
2015 TAxis *xAxisDriftLength = hres->GetXaxis(); // driftlength / z - axis
2016 TAxis *yAxisAngle = hres->GetYaxis(); // angle axis
2017 TAxis *zAxisDelta = hres->GetZaxis(); // delta axis
2018 TAxis *zAxisRms = hrms->GetZaxis(); // rms axis
2019
2020 // loop over all angle bins
2021 for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
2022 angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
2023 // loop over all driftlength bins
2024 for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
2025 zCenter = xAxisDriftLength->GetBinCenter(ibinxDL);
2026 zSigma = xAxisDriftLength->GetBinWidth(ibinxDL);
2027 angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
2028 zMean = zCenter; // changens, when more statistic is accumulated
2029 angleMean = angleCenter; // changens, when more statistic is accumulated
2030
2031 // create 2 1D-Histograms, projectionRes and projectionRms
2032 // these histograms are delta histograms for given direction, padSize, chargeBin,
2033 // angleBin and driftLengthBin
2034 // later on they will be fitted with a gausian, its sigma is the resoltuion...
2035 sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
2036 // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2037 projectionRes->SetNameTitle(name, name);
2038 sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
2039 // TH1D * projectionRms = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2040 projectionRms->SetNameTitle(name, name);
2041
2042 projectionRes->Reset();
2043 projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2044 projectionRms->Reset();
2045 projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2046 projectionRes->SetDirectory(0);
2047 projectionRms->SetDirectory(0);
2048
2049 Double_t entries = 0;
2050 Int_t nbins = 0; // counts, how many bins were accumulated
2051 zMean = 0;
2052 angleMean = 0;
2053 entries = 0;
2054
2055 // fill projectionRes and projectionRms for given dim, ipad and iq,
2056 // as well as for given angleBin and driftlengthBin
2057 // if this gives not enough statistic, include neighbourhood
2058 // (angle and driftlength) successifely
2059 for (Int_t dbin = 0; dbin <= 8; dbin++){ // delta-bins around centered angleBin and driftlengthBin
2060 for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) { // delta-bins in angle direction
2061 for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
2062 if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; // add each bin only one time !
2063 Int_t binx2 = ibinxDL + dbinx2; // position variable in x (driftlength) direction
2064 Int_t biny2 = ibinyAngle + dbiny2; // position variable in y (angle) direction
2065 if (binx2 < 1 || biny2 < 1) continue; // don't go out of the histogram!
2066 if (binx2 >= xAxisDriftLength->GetNbins()) continue; // don't go out of the histogram!
2067 if (biny2 >= yAxisAngle->GetNbins()) continue; // don't go out of the histogram!
2068 nbins++; // count the number of accumulated bins
2069 // Fill resolution histo
2070 for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
2071 // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3); // unused variable
2072 projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
2073 entries += hres->GetBinContent(binx2, biny2, ibin3);
2074 zMean += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
2075 angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
2076 } // ibin3 loop
2077 // fill RMS histo
2078 for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
2079 projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
2080 }
2081 } //dbinx2 loop
2082 if (entries > minEntries) break; // enough statistic accumulated
2083 } // dbiny2 loop
2084 if (entries > minEntries) break; // enough statistic accumulated
2085 } // dbin loop
2086 if ( entries< minEntries) continue; // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
2087 zMean /= entries;
2088 angleMean /= entries;
2089
2090 if (entries > minEntries) {
2091 // when enough statistic is accumulated
2092 // fit Delta histograms with a gausian
2093 // of the gausian is the resolution (resol), its fit error is sigma
2094 // store also mean and RMS of the histogram
2095 Float_t xmin = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2096 Float_t xmax = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2097
2098// projectionRes->Fit("gaus", "q0", "", xmin, xmax);
2099// Float_t resol = projectionRes->GetFunction("gaus")->GetParameter(2);
2100// Float_t sigma = projectionRes->GetFunction("gaus")->GetParError(2);
2101 fitFunction->SetMaximum(xmax);
2102 fitFunction->SetMinimum(xmin);
2103 projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
2104 Float_t resol = fitFunction->GetParameter(2);
2105 Float_t sigma = fitFunction->GetParError(2);
2106
2107 Float_t meanR = projectionRes->GetMean();
2108 Float_t sigmaR = projectionRes->GetRMS();
2109 // fit also RMS histograms with a gausian
2110 // store mean and sigma of the gausian in rmsMean and rmsSigma
2111 // store also the fit errors in errorRMS and errorSigma
2112 xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2113 xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2114
2115// projectionRms->Fit("gaus","q0","",xmin,xmax);
2116// Float_t rmsMean = projectionRms->GetFunction("gaus")->GetParameter(1);
2117// Float_t rmsSigma = projectionRms->GetFunction("gaus")->GetParameter(2);
2118// Float_t errorRMS = projectionRms->GetFunction("gaus")->GetParError(1);
2119// Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
2120 projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
2121 Float_t rmsMean = fitFunction->GetParameter(1);
2122 Float_t rmsSigma = fitFunction->GetParameter(2);
2123 Float_t errorRMS = fitFunction->GetParError(1);
2124 Float_t errorSigma = fitFunction->GetParError(2);
2125
2126 Float_t length = 0.75;
2127 if (ipad == 1) length = 1;
2128 if (ipad == 2) length = 1.5;
2129
2130 fTreeResol<<"Resol"<<
2131 "Entries="<<entries<< // number of entries for this resolution point
2132 "nbins="<<nbins<< // number of bins that were accumulated
2133 "Dim="<<idim<< // direction, Dim==0: y-direction, Dim==1: z-direction
2134 "Pad="<<ipad<< // padSize; short, medium and long
2135 "Length="<<length<< // pad length, 0.75, 1, 1.5
2136 "QMean="<<qMean<< // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
2137 "Zc="<<zCenter<< // center of middle bin in drift direction
2138 "Zm="<<zMean<< // mean dirftlength for accumulated Delta-Histograms
2139 "Zs="<<zSigma<< // width of driftlength bin
2140 "AngleC="<<angleCenter<< // center of middle bin in Angle-Direction
2141 "AngleM="<<angleMean<< // mean angle for accumulated Delta-Histograms
2142 "AngleS="<<angleSigma<< // width of Angle-bin
2143 "Resol="<<resol<< // sigma for gaus fit through Delta-Histograms
2144 "Sigma="<<sigma<< // error of sigma for gaus fit through Delta Histograms
2145 "MeanR="<<meanR<< // mean of the Delta-Histogram
2146 "SigmaR="<<sigmaR<< // rms of the Delta-Histogram
2147 "RMSm="<<rmsMean<< // mean of the gaus fit through RMS-Histogram
2148 "RMSs="<<rmsSigma<< // sigma of the gaus fit through RMS-Histogram
2149 "RMSe0="<<errorRMS<< // error of mean of gaus fit in RMS-Histogram
2150 "RMSe1="<<errorSigma<< // error of sigma of gaus fit in RMS-Histogram
2151 "\n";
ae28e92e 2152 if (GetDebugLevel() > 5) {
10757ee9 2153 projectionRes->SetDirectory(fTreeResol.GetFile());
2154 projectionRes->Write(projectionRes->GetName());
2155 projectionRes->SetDirectory(0);
2156 projectionRms->SetDirectory(fTreeResol.GetFile());
2157 projectionRms->Write(projectionRms->GetName());
2158 projectionRes->SetDirectory(0);
2159 }
2160 } // if (projectionRes->GetSum() > minEntries)
2161 } // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
2162 } // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
2163
2164 } // iq-loop
2165 } // ipad-loop
2166 } // idim-loop
2167 delete projectionRes;
2168 delete projectionRms;
2169
2170// TFile resolFile(fTreeResol.GetFile());
2171 TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
2172 fileInfo.Write("fileInfo");
2173// resolFile.Close();
2174// fTreeResol.GetFile()->Close();
ae28e92e 2175 if (GetDebugLevel() > -1) cout << endl;
2176 if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
10757ee9 2177 gSystem->ChangeDirectory("..");
2178}
2179
2180
10757ee9 2181
10757ee9 2182
2183
2184Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
2185 //
2186 // function to merge several AliTPCcalibTracks objects after PROOF calculation
2187 // The object's histograms are merged via their merge functions
2188 // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
2189 //
2190
ae28e92e 2191 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl;
10757ee9 2192 if (!collectionList) return 0;
2193 if (collectionList->IsEmpty()) return -1;
2194
ae28e92e 2195 if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
2196 if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
10757ee9 2197 collectionList->Print();
2198
2199 // create a list for each data member
2200 TList* deltaYList = new TList;
2201 TList* deltaZList = new TList;
2202 TList* arrayAmpRowList = new TList;
2203 TList* rejectedTracksList = new TList;
2204 TList* hclusList = new TList;
2205 TList* clusterCutHistoList = new TList;
2206 TList* arrayAmpList = new TList;
2207 TList* arrayQDYList = new TList;
2208 TList* arrayQDZList = new TList;
2209 TList* arrayQRMSYList = new TList;
2210 TList* arrayQRMSZList = new TList;
2211 TList* arrayChargeVsDriftlengthList = new TList;
2212 TList* calPadRegionChargeVsDriftlengthList = new TList;
2213 TList* hclusterPerPadrowList = new TList;
2214 TList* hclusterPerPadrowRawList = new TList;
2215 TList* resolYList = new TList;
2216 TList* resolZList = new TList;
2217 TList* rMSYList = new TList;
2218 TList* rMSZList = new TList;
2219
2220// TList* nRowsList = new TList;
2221// TList* nSectList = new TList;
2222// TList* fileNoList = new TList;
2223
2224 TIterator *listIterator = collectionList->MakeIterator();
2225 AliTPCcalibTracks *calibTracks = 0;
ae28e92e 2226 if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
10757ee9 2227 Int_t counter = 0;
a4c5fca9 2228 while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
10757ee9 2229 // loop over all entries in the collectionList and get dataMembers into lists
2230 if (!calibTracks) continue;
2231
2232 deltaYList->Add( calibTracks->GetfDeltaY() );
2233 deltaZList->Add( calibTracks->GetfDeltaZ() );
2234 arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
2235 arrayAmpList->Add(calibTracks->GetfArrayAmp());
2236 arrayQDYList->Add(calibTracks->GetfArrayQDY());
2237 arrayQDZList->Add(calibTracks->GetfArrayQDZ());
2238 arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
2239 arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
2240 resolYList->Add(calibTracks->GetfResolY());
2241 resolZList->Add(calibTracks->GetfResolZ());
2242 rMSYList->Add(calibTracks->GetfRMSY());
2243 rMSZList->Add(calibTracks->GetfRMSZ());
2244 arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
2245 calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
2246 hclusList->Add(calibTracks->GetfHclus());
2247 rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
2248 clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
2249 hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
2250 hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
a4c5fca9 2251 //
2252 if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
2253 fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
2254 // fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
10757ee9 2255 counter++;
ae28e92e 2256 if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
af6a50bb 2257 AddHistos(calibTracks);
10757ee9 2258 }
2259
10757ee9 2260
2261 // merge data members
ae28e92e 2262 if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
10757ee9 2263 fDeltaY->Merge(deltaYList);
2264 fDeltaZ->Merge(deltaZList);
2265 fHclus->Merge(hclusList);
2266 fClusterCutHisto->Merge(clusterCutHistoList);
2267 fRejectedTracksHisto->Merge(rejectedTracksList);
2268 fHclusterPerPadrow->Merge(hclusterPerPadrowList);
2269 fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
2270
2271 TObjArray* objarray = 0;
2272 TH1* hist = 0;
2273 TList* histList = 0;
2274 TIterator *objListIterator = 0;
2275
ae28e92e 2276 if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
10757ee9 2277 // merge fArrayAmpRows
2278 for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) { // loop over data member, i<72
2279 objListIterator = arrayAmpRowList->MakeIterator();
2280 histList = new TList;
2281 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2282 // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
2283 if (!objarray) continue;
2284 hist = (TProfile*)(objarray->At(i));
2285 histList->Add(hist);
2286 }
2287 ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
2288 delete histList;
2289 delete objListIterator;
2290 }
2291
ae28e92e 2292 if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
10757ee9 2293 // merge fArrayAmps
2294 for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) { // loop over data member, i<72
d23650ed 2295 TIterator *cobjListIterator = arrayAmpList->MakeIterator();
10757ee9 2296 histList = new TList;
d23650ed 2297 while (( objarray = (TObjArray*)cobjListIterator->Next() )) {
10757ee9 2298 // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
2299 if (!objarray) continue;
2300 hist = (TH1F*)(objarray->At(i));
2301 histList->Add(hist);
2302 }
2303 ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
2304 delete histList;
d23650ed 2305 delete cobjListIterator;
10757ee9 2306 }
2307
ae28e92e 2308 if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
10757ee9 2309 // merge fArrayQDY
2310 for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
2311 objListIterator = arrayQDYList->MakeIterator();
2312 histList = new TList;
2313 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2314 // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
2315 if (!objarray) continue;
2316 hist = (TH3F*)(objarray->At(i));
2317 histList->Add(hist);
2318 }
2319 ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
2320 delete histList;
2321 delete objListIterator;
2322 }
2323
ae28e92e 2324 if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
10757ee9 2325 // merge fArrayQDZ
2326 for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2327 objListIterator = arrayQDZList->MakeIterator();
2328 histList = new TList;
2329 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2330 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2331 if (!objarray) continue;
2332 hist = (TH3F*)(objarray->At(i));
2333 histList->Add(hist);
2334 }
2335 ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
2336 delete histList;
2337 delete objListIterator;
2338 }
2339
ae28e92e 2340 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
10757ee9 2341 // merge fArrayQRMSY
2342 for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
2343 objListIterator = arrayQRMSYList->MakeIterator();
2344 histList = new TList;
2345 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2346 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2347 if (!objarray) continue;
2348 hist = (TH3F*)(objarray->At(i));
2349 histList->Add(hist);
2350 }
2351 ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
2352 delete histList;
2353 delete objListIterator;
2354 }
2355
ae28e92e 2356 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
10757ee9 2357 // merge fArrayQRMSZ
2358 for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2359 objListIterator = arrayQRMSZList->MakeIterator();
2360 histList = new TList;
2361 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2362 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2363 if (!objarray) continue;
2364 hist = (TH3F*)(objarray->At(i));
2365 histList->Add(hist);
2366 }
2367 ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
2368 delete histList;
2369 delete objListIterator;
2370 }
2371
ae28e92e 2372 if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
10757ee9 2373 // merge fArrayChargeVsDriftlength
2374 for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
2375 objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
2376 histList = new TList;
2377 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2378 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
2379 if (!objarray) continue;
2380 hist = (TProfile*)(objarray->At(i));
2381 histList->Add(hist);
2382 }
2383 ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
2384 delete histList;
2385 delete objListIterator;
2386 }
2387
ae28e92e 2388 if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
10757ee9 2389 // merge fcalPadRegionChargeVsDriftlength
2390 AliTPCCalPadRegion *cpr = 0x0;
2391
2392 /*
2393 TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
2394 while (hist = (TProfile*)regionIterator->Next()) {
2395 // loop over all calPadRegion's in destination calibTracks object
2396 objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2397 while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
2398
2399
2400 hist->Merge(...);
2401 }
2402 */
2403
2404 for (UInt_t isec = 0; isec < 36; isec++) {
2405 for (UInt_t padSize = 0; padSize < 3; padSize++){
2406 objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2407 histList = new TList;
2408 while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
2409 // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object
2410 if (!cpr) continue;
2411 hist = (TProfile*)cpr->GetObject(isec, padSize);
2412 histList->Add(hist);
2413 }
2414 ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
2415 delete histList;
2416 delete objListIterator;
2417 }
2418 }
2419
2420
2421
2422
ae28e92e 2423 if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
10757ee9 2424 // merge fResolY
2425 for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
2426 objListIterator = resolYList->MakeIterator();
2427 histList = new TList;
2428 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2429 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2430 if (!objarray) continue;
2431 hist = (TH3F*)(objarray->At(i));
2432 histList->Add(hist);
2433 }
2434 ((TH3F*)(fResolY->At(i)))->Merge(histList);
2435 delete histList;
2436 delete objListIterator;
2437 }
2438
2439 // merge fResolZ
2440 for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2441 objListIterator = resolZList->MakeIterator();
2442 histList = new TList;
2443 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2444 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2445 if (!objarray) continue;
2446 hist = (TH3F*)(objarray->At(i));
2447 histList->Add(hist);
2448 }
2449 ((TH3F*)(fResolZ->At(i)))->Merge(histList);
2450 delete histList;
2451 delete objListIterator;
2452 }
2453
2454 // merge fRMSY
2455 for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
2456 objListIterator = rMSYList->MakeIterator();
2457 histList = new TList;
2458 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2459 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2460 if (!objarray) continue;
2461 hist = (TH3F*)(objarray->At(i));
2462 histList->Add(hist);
2463 }
2464 ((TH3F*)(fRMSY->At(i)))->Merge(histList);
2465 delete histList;
2466 delete objListIterator;
2467 }
2468
2469 // merge fRMSZ
2470 for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2471 objListIterator = rMSZList->MakeIterator();
2472 histList = new TList;
2473 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2474 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2475 if (!objarray) continue;
2476 hist = (TH3F*)(objarray->At(i));
2477 histList->Add(hist);
2478 }
2479 ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
2480 delete histList;
2481 delete objListIterator;
2482 }
2483
2484 delete deltaYList;
2485 delete deltaZList;
2486 delete arrayAmpRowList;
2487 delete arrayAmpList;
2488 delete arrayQDYList;
2489 delete arrayQDZList;
2490 delete arrayQRMSYList;
2491 delete arrayQRMSZList;
2492 delete resolYList;
2493 delete resolZList;
2494 delete rMSYList;
2495 delete rMSZList;
10757ee9 2496 delete listIterator;
2497
ae28e92e 2498 if (GetDebugLevel() > 0) cout << "merging done!" << endl;
10757ee9 2499
2500 return 1;
2501}
2502
2503
2504AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
2505 //
2506 // function to test AliTPCcalibTrack::Merge:
2507 // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
2508 // this object is appended 'nCalTracks' times to a TList
2509 // A new AliTPCcalibTrack object is created which merges the list
2510 // this object is returned
2511 //
2512 /*
2513 // .L AliTPCcalibTracks.cxx+g
2514 .L libTPCcalib.so
2515 TFile f("Output.root");
2516 AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
2517 //f.Close();
2518 TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
2519 AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param");
2520 clusterParamFile.Close();
2521
2522 AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
2523 */
2524 TList *list = new TList();
2525 if (ct == 0 || clusterParam == 0) return 0;
2526 cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
2527 for (Int_t i = 0; i < nCalTracks; i++) {
2528 if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
b42cf9ed 2529 list->Add(new AliTPCcalibTracks(*ct));
10757ee9 2530 }
2531
2532 // only for check at the end
b42cf9ed 2533 AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
10757ee9 2534 Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
2535// Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
2536
2537 cout << "The list contains " << list->GetEntries() << " entries. " << endl;
2538
2539
2540 AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
2541 AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
2542 cal->Merge(list);
2543
2544 cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
2545 cal->GetfArrayAmpRow()->At(5)->Print();
2546 Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
2547
2548 cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
2549 cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " << calEntries << endl;
2550 printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n",
2551 calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
2552
2553 return cal;
2554
2555}
2556
2557
1bfaa9e9 2558void AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints){
9c171b96 2559 //
2560 // Make position corrections
2561 // for the moment Only using debug streamer
2562 // chainres - debug tree
2563 // param - parameters to be updated
2564 // maxPoints - maximal number of points using for fit
2565 // verbose - print info flag
2e5bcb67 2566 //
2567 // Current implementation - only using debug streamers
2568 //
2569
2570 /*
9c171b96 2571 //Defaults
2572 Int_t maxPoints=100000;
2573 */
2e5bcb67 2574 /*
2575 Usage:
2576 //0. Load libraries
2577 gSystem->Load("libANALYSIS");
2578 gSystem->Load("libSTAT");
2579 gSystem->Load("libTPCcalib");
2580
2581
2582 //1. Load Parameters to be modified:
2583 //e.g:
2584
162637e4 2585 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
2e5bcb67 2586 AliCDBManager::Instance()->SetRun(0)
2587 AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
2588
2589 //2. Load chain from debug streamers
2590 //
2591 //e.g
2592 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2593 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2594 AliXRDPROOFtoolkit tool;
2595 TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2596 chainres->Lookup();
2597 //3. Do fits and store results
2598 //
2599 AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
2600 TFile f("paramout.root","recreate");
2601 param->Write("clusterParam");
2602 f.Close();
2603 //4. Verification
2604 TFile f2("paramout.root");
2605 AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
2606 param2->SetInstance(param2);
2607 chainres->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line
2608
2609 */
2610
2611
2612 TStatToolkit toolkit;
2613 Double_t chi2;
2614 TVectorD fitParamY0;
2615 TVectorD fitParamY1;
2616 TVectorD fitParamZ0;
2617 TVectorD fitParamZ1;
2618 TMatrixD covMatrix;
2619 Int_t npoints;
2620
2621 chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2622 chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2623 chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2624 chainres->SetAlias("st","(sin(dt)-dt)");
2625 //
2626 chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2e5bcb67 2627 //
2628 //
2629 //
2630 TCut cutA("1");
2631 TString fstringY="";
2632 //
2633 fstringY+="(dp)++"; //1
2634 fstringY+="(dp)*di++"; //2
108953e9 2635 fstringY+="(sp)++"; //3
2636 fstringY+="(sp)*di++"; //4
2e5bcb67 2637 TString fstringZ="";
2638 fstringZ+="(dt)++"; //1
2639 fstringZ+="(dt)*di++"; //2
108953e9 2640 fstringZ+="(st)++"; //3
2641 fstringZ+="(st)*di++"; //4
2e5bcb67 2642 //
2643 // Z corrections
2644 //
2645 TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
2646 printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
798017c7 2647 param->PosZcor(0) = (TVectorD*) fitParamZ0.Clone();
2e5bcb67 2648 //
2649 TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
2650 //
2651 printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
798017c7 2652 param->PosZcor(1) = (TVectorD*) fitParamZ1.Clone();
2653 param->PosZcor(2) = (TVectorD*) fitParamZ1.Clone();
2e5bcb67 2654 //
2655 // Y corrections
2656 //
2657 TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
2658 printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
798017c7 2659 param->PosYcor(0) = (TVectorD*) fitParamY0.Clone();
2e5bcb67 2660
2661
2662 TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
2663 //
2664 printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
798017c7 2665 param->PosYcor(1) = (TVectorD*) fitParamY1.Clone();
2666 param->PosYcor(2) = (TVectorD*) fitParamY1.Clone();
2e5bcb67 2667 //
2668 //
2669 //
2670 chainres->SetAlias("fitZ0",strZ0->Data());
2671 chainres->SetAlias("fitZ1",strZ1->Data());
2672 chainres->SetAlias("fitY0",strY0->Data());
2673 chainres->SetAlias("fitY1",strY1->Data());
2674 // chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);
9c171b96 2675}
2676
10757ee9 2677
2678
af6a50bb 2679void AliTPCcalibTracks::MakeHistos(){
2680 //
2681 ////make THnSparse
2682 //
2683 //THnSparse *fHisDeltaY; // THnSparse - delta Y
2684 //THnSparse *fHisDeltaZ; // THnSparse - delta Z
2685 //THnSparse *fHisRMSY; // THnSparse - rms Y
2686 //THnSparse *fHisRMSZ; // THnSparse - rms Z
2687 //THnSparse *fHisQmax; // THnSparse - qmax
2688 //THnSparse *fHisQtot; // THnSparse - qtot
2689 // cluster performance bins
2690 // 0 - variable of interest
2691 // 1 - pad type - 0- short 1-medium 2-long pads
2692 // 2 - drift length - drift length -0-1
2693 // 3 - Qmax - Qmax - 2- 400
2694 // 4 - cog - COG position - 0-1
2695 // 5 - tan(phi) - local y angle
2696 // 6 - tan(theta) - local z angle
2697 // 7 - sector - sector number
2698 Double_t xminTrack[8], xmaxTrack[8];
2699 Int_t binsTrack[8];
2700 TString axisName[8];
2701
2702 //
2703 binsTrack[0]=100;
2704 axisName[0] ="var";
2705
2706 binsTrack[1] =3;
2707 xminTrack[1] =0; xmaxTrack[1]=3;
2708 axisName[1] ="pad type";
2709 //
2710 binsTrack[2] =10;
2711 xminTrack[2] =0; xmaxTrack[2]=1;
2712 axisName[2] ="drift length";
2713 //
2714 binsTrack[3] =10;
2715 xminTrack[3] =1; xmaxTrack[3]=400;
2716 axisName[3] ="Qmax";
2717 //
2718 binsTrack[4] =10;
2719 xminTrack[4] =0; xmaxTrack[4]=1;
2720 axisName[4] ="cog";
2721 //
2722 binsTrack[5] =10;
2723 xminTrack[5] =0; xmaxTrack[5]=2;
2724 axisName[5] ="tan(phi)";
2725 //
2726 binsTrack[6] =10;
2727 xminTrack[6] =0; xmaxTrack[6]=2;
2728 axisName[6] ="tan(theta)";
2729 //
2730 xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
2731 fHisDeltaY=new THnSparseS("#Delta_{y} (cm)","#Delta_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2732 xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
2733 fHisDeltaZ=new THnSparseS("#Delta_{z} (cm)","#Delta_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2734 xminTrack[0] =0.; xmaxTrack[0]=0.5;
2735 fHisRMSY=new THnSparseS("#RMS_{y} (cm)","#RMS_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2736 xminTrack[0] =0.; xmaxTrack[0]=0.5;
2737 fHisRMSZ=new THnSparseS("#RMS_{z} (cm)","#RMS_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2738 xminTrack[0] =0.; xmaxTrack[0]=100;
2739 fHisQmax=new THnSparseS("Qmax (ADC)","Qmax (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
2740
2741 xminTrack[0] =0.; xmaxTrack[0]=250;
2742 fHisQtot=new THnSparseS("Qtot (ADC)","Qtot (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
2743 BinLogX(fHisDeltaY,3);
2744 BinLogX(fHisDeltaZ,3);
2745 BinLogX(fHisRMSY,3);
2746 BinLogX(fHisRMSZ,3);
2747 BinLogX(fHisQmax,3);
2748 BinLogX(fHisQtot,3);
2749
2750}
2751
2752void AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
2753 //
2754 // Add histograms
2755 //
2756 if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
2757 if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
2758 if (calib->fHisRMSY) fHisRMSY->Add(calib->fHisRMSY);
2759 if (calib->fHisRMSZ) fHisRMSZ->Add(calib->fHisRMSZ);
2760}