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