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