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