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