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