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