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