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