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