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