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