]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug fix (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTracks.cxx
CommitLineData
10757ee9 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17///////////////////////////////////////////////////////////////////////////////
18// //
19// Class to analyse tracks for calibration //
20// to be used as a component in AliTPCSelectorTracks //
21// In the constructor you have to specify name and title //
22// to get the Object out of a file. //
23// The parameter 'clusterParam', a AliTPCClusterParam object //
24// (needed for TPC cluster error and shape parameterization) //
25// Normally you get this object out of the file 'TPCClusterParam.root' //
26// In the parameter 'cuts' the cuts are specified, that decide //
27// weather a track will be accepted for calibration or not. //
28// //
29//
30// The data flow:
31//
32/*
33 Raw Data -> Local Reconstruction -> Tracking -> Calibration -> RefData (component itself)
34 Offline/HLT Offline/HLT OCDB entries (AliTPCClusterParam)
35*/
36
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);
353 his1 = new TH1F(chname,chname,250,0,500); // valgrind
354 his1->SetXTitle("Max Amplitude (ADC)");
355 fArrayAmp->AddAt(his1,i);
356 sprintf(chname,"Amp_Sector%d",i+36);
357 his1 = new TH1F(chname,chname,200,0,600); // valgrind 3 13,408,208 bytes in 229 blocks are still reachable
358 his1->SetXTitle("Max Amplitude (ADC)");
359 fArrayAmp->AddAt(his1,i+36);
360
361 // driftlength
362 sprintf(chname, "driftlengt vs. charge, ROC %i", i);
363 prof1 = new TProfile(chname, chname, 500, 0, 250);
364 prof1->SetYTitle("Charge");
365 prof1->SetXTitle("Driftlength");
366 fArrayChargeVsDriftlength->AddAt(prof1,i);
367 sprintf(chname, "driftlengt vs. charge, ROC %i", i+36);
368 prof1 = new TProfile(chname, chname, 500, 0, 250);
369 prof1->SetYTitle("Charge");
370 prof1->SetXTitle("Driftlength");
371 fArrayChargeVsDriftlength->AddAt(prof1,i+36);
372 }
373
374 TH1::AddDirectory(kFALSE);
375
376 fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1);
377 fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1);
378
379 fResolY = new TObjArray(3);
380 fResolZ = new TObjArray(3);
381 fRMSY = new TObjArray(3);
382 fRMSZ = new TObjArray(3);
383 TH3F * his3D;
384 //
385 his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
386 fResolY->AddAt(his3D,0);
387 his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
388 fResolY->AddAt(his3D,1);
389 his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
390 fResolY->AddAt(his3D,2);
391 //
392 his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
393 fResolZ->AddAt(his3D,0);
394 his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
395 fResolZ->AddAt(his3D,1);
396 his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
397 fResolZ->AddAt(his3D,2);
398 //
399 his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
400 fRMSY->AddAt(his3D,0);
401 his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
402 fRMSY->AddAt(his3D,1);
403 his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
404 fRMSY->AddAt(his3D,2);
405 //
406 his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
407 fRMSZ->AddAt(his3D,0);
408 his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
409 fRMSZ->AddAt(his3D,1);
410 his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
411 fRMSZ->AddAt(his3D,2);
412 //
413
414 TH1::AddDirectory(kFALSE);
415
416 fArrayQDY = new TObjArray(300);
417 fArrayQDZ = new TObjArray(300);
418 fArrayQRMSY = new TObjArray(300);
419 fArrayQRMSZ = new TObjArray(300);
420 for (Int_t iq = 0; iq <= 10; iq++){
421 for (Int_t ipad = 0; ipad < 3; ipad++){
422 Int_t bin = GetBin(iq, ipad);
423 Float_t qmean = GetQ(bin);
424 char name[200];
425 sprintf(name,"ResolY Pad%d Qmiddle%f",ipad, qmean);
426 his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
427 fArrayQDY->AddAt(his3D, bin);
428 sprintf(name,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
429 his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
430 fArrayQDZ->AddAt(his3D, bin);
431 sprintf(name,"RMSY Pad%d Qmiddle%f",ipad, qmean);
432 his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
433 fArrayQRMSY->AddAt(his3D, bin);
434 sprintf(name,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
435 his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
436 fArrayQRMSZ->AddAt(his3D, bin);
437 }
438 }
439
440 fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region");
441 TProfile *tempProf;
442 for (UInt_t padSize = 0; padSize < 3; padSize++) {
443 for (UInt_t isector = 0; isector < 36; isector++) {
444 if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector);
445 if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium pads", isector);
446 if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long pads", isector);
447 tempProf = new TProfile(chname, chname, 500, 0, 250);
448 tempProf->SetYTitle("Charge");
449 tempProf->SetXTitle("Driftlength");
450 fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize);
451 }
452 }
453
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;
827 Double_t chi2 = 0;
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);
834 chi2 += chi2Y(0, 0);
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);
841 chi2 += chi2Z(0, 0);
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
853 TTreeSRedirector *cstream = GetDebugStreamer();
854 if (cstream){
855 (*cstream)<<"Cut8"<<
856 "chi2="<<chi2<<
857 "\n";
858 }
10757ee9 859 }
860
861 // current padrow has no kink
862
863 // get fit parameters from pol2 fit:
864 Double_t paramY[4], paramZ[4];
9f0beaf7 865 paramY[0] = fFitterParY.GetParameter(0);
866 paramY[1] = fFitterParY.GetParameter(1);
867 paramY[2] = fFitterParY.GetParameter(2);
868 paramZ[0] = fFitterParZ.GetParameter(0);
869 paramZ[1] = fFitterParZ.GetParameter(1);
870 paramZ[2] = fFitterParZ.GetParameter(2);
10757ee9 871
872 Double_t tracky = paramY[0];
873 Double_t trackz = paramZ[0];
874 Float_t deltay = tracky - cluster0->GetY();
875 Float_t deltaz = trackz - cluster0->GetZ();
876 Float_t angley = paramY[1] - paramY[0] / xref;
877 Float_t anglez = paramZ[1];
878
879 Float_t max = cluster0->GetMax();
880 UInt_t isegment = cluster0->GetDetector() % 36;
881 Int_t padSize = 0; // short pads
882 if (cluster0->GetDetector() >= 36) {
883 padSize = 1; // medium pads
884 if (cluster0->GetRow() > 63) padSize = 2; // long pads
885 }
886
887 // =========================================
888 // wirte collected information to histograms
889 // =========================================
890
891 TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector);
892 if ( irow >= kFirstLargePad) max /= kLargePadSize;
893 profAmpRow->Fill( (Double_t)cluster0->GetRow(), max );
894 TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector);
895 hisAmp->Fill(max);
896
897 // remove the following two lines one day:
898 TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
899 profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
900
901 TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
902 profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
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"<<
931 "padSize="<<padSize<<
932 "angley="<<angley<<
933 "anglez="<<anglez<<
934 "zdr="<<zdrift<<
935 "dy="<<deltay<<
936 "dz="<<deltaz<<
937 "sy="<<sy<<
938 "sz="<<sz<<
939 "\n";
940 }
941
10757ee9 942 if (useForResol){
943 fDeltaY->Fill(deltay);
944 fDeltaZ->Fill(deltaz);
945 his3 = (TH3F*)fResolY->At(padSize);
946 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
947 his3 = (TH3F*)fResolZ->At(padSize);
948 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
949 his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
950 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
951 his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
952 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
953 }
954
955 //=============================================================================================
956
957 if (useForResol && nclFound > 2 * kMinRatio * kDelta
ae28e92e 958 && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
2acad464 959 if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
10757ee9 960 FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
961 } // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
962
963 } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
964} // FillResolutionHistoLocal(...)
965
966
967
968void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t angley, Float_t anglez, Int_t nclFound, Int_t kDelta) {
969 //
970 // - debug part of FillResolutionHistoLocal -
971 // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
ae28e92e 972 // called only for GetStreamLevel() > 4
10757ee9 973 // fill resolution trees
974 //
975
976 Int_t sector = cluster0->GetDetector();
977 Float_t xref = cluster0->GetX();
978 Int_t padSize = 0; // short pads
979 if (cluster0->GetDetector() >= 36) {
980 padSize = 1; // medium pads
981 if (cluster0->GetRow() > 63) padSize = 2; // long pads
982 }
983
9f0beaf7 984 static TLinearFitter fitY0(3, "pol2");
985 static TLinearFitter fitZ0(3, "pol2");
986 static TLinearFitter fitY2(5, "hyp4");
987 static TLinearFitter fitZ2(5, "hyp4");
988 static TLinearFitter fitY2Q(5, "hyp4");
989 static TLinearFitter fitZ2Q(5, "hyp4");
990 static TLinearFitter fitY2S(5, "hyp4");
991 static TLinearFitter fitZ2S(5, "hyp4");
992 fitY0.ClearPoints();
993 fitZ0.ClearPoints();
994 fitY2.ClearPoints();
995 fitZ2.ClearPoints();
996 fitY2Q.ClearPoints();
997 fitZ2Q.ClearPoints();
998 fitY2S.ClearPoints();
999 fitZ2S.ClearPoints();
1000
1001 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
1002 // loop over irow +- kDelta rows (neighboured rows)
1003 //
1004 //
1005 if (idelta == 0) continue;
1006 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
1007 AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
1008 if (!cluster) continue;
1009 if (cluster->GetType() < 0) continue;
1010 if (cluster->GetDetector() != sector) continue;
1011 Double_t x = cluster->GetX() - xref;
1012 Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
1013 Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
1014 //
1015 Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1016 Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1017 Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
10757ee9 1018 TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
9f0beaf7 1019 Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
10757ee9 1020 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
9f0beaf7 1021 Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1022 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1023 TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
1024 Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1025 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1026 TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
1027 sigmaYS = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
1028 sigmaZS = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
1029 //
1030 if (kDelta != 0){
1031 fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
1032 fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
1033 }
1034 Double_t xxx[4];
1035 xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
1036 xxx[1] = x;
1037 xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
1038 xxx[3] = x * x;
1039 fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
1040 fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
1041 fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
1042 fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
1043 fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
1044 fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
1045 //
1046 } // neigbouhood-loop
10757ee9 1047 //
9f0beaf7 1048 fitY0.Eval();
1049 fitZ0.Eval();
1050 fitY2.Eval();
1051 fitZ2.Eval();
1052 fitY2Q.Eval();
1053 fitZ2Q.Eval();
1054 fitY2S.Eval();
1055 fitZ2S.Eval();
1056 Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
1057 Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
1058 Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
1059 Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
1060 Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
1061 Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
1062 Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
1063 Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
1064 //
1065 static TVectorD parY0(3);
1066 static TMatrixD matY0(3, 3);
1067 static TVectorD parZ0(3);
1068 static TMatrixD matZ0(3, 3);
1069 fitY0.GetParameters(parY0);
1070 fitY0.GetCovarianceMatrix(matY0);
1071 fitZ0.GetParameters(parZ0);
1072 fitZ0.GetCovarianceMatrix(matZ0);
1073 //
1074 static TVectorD parY2(5);
1075 static TMatrixD matY2(5,5);
1076 static TVectorD parZ2(5);
1077 static TMatrixD matZ2(5,5);
1078 fitY2.GetParameters(parY2);
1079 fitY2.GetCovarianceMatrix(matY2);
1080 fitZ2.GetParameters(parZ2);
1081 fitZ2.GetCovarianceMatrix(matZ2);
1082 //
1083 static TVectorD parY2Q(5);
1084 static TMatrixD matY2Q(5,5);
1085 static TVectorD parZ2Q(5);
1086 static TMatrixD matZ2Q(5,5);
1087 fitY2Q.GetParameters(parY2Q);
1088 fitY2Q.GetCovarianceMatrix(matY2Q);
1089 fitZ2Q.GetParameters(parZ2Q);
1090 fitZ2Q.GetCovarianceMatrix(matZ2Q);
1091 static TVectorD parY2S(5);
1092 static TMatrixD matY2S(5,5);
1093 static TVectorD parZ2S(5);
1094 static TMatrixD matZ2S(5,5);
1095 fitY2S.GetParameters(parY2S);
1096 fitY2S.GetCovarianceMatrix(matY2S);
1097 fitZ2S.GetParameters(parZ2S);
1098 fitZ2S.GetCovarianceMatrix(matZ2S);
1099 Float_t sigmaY0 = TMath::Sqrt(matY0(0,0));
1100 Float_t sigmaZ0 = TMath::Sqrt(matZ0(0,0));
1101 Float_t sigmaDY0 = TMath::Sqrt(matY0(1,1));
1102 Float_t sigmaDZ0 = TMath::Sqrt(matZ0(1,1));
1103 Float_t sigmaY2 = TMath::Sqrt(matY2(1,1));
1104 Float_t sigmaZ2 = TMath::Sqrt(matZ2(1,1));
1105 Float_t sigmaDY2 = TMath::Sqrt(matY2(3,3));
1106 Float_t sigmaDZ2 = TMath::Sqrt(matZ2(3,3));
1107 Float_t sigmaY2Q = TMath::Sqrt(matY2Q(1,1));
1108 Float_t sigmaZ2Q = TMath::Sqrt(matZ2Q(1,1));
1109 Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
1110 Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
1111 Float_t sigmaY2S = TMath::Sqrt(matY2S(1,1));
1112 Float_t sigmaZ2S = TMath::Sqrt(matZ2S(1,1));
1113 Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
1114 Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
1115
1116 // Error parameters
1117 Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
1118 Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
1119 //
1120 Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1121 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1122 Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1123 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1124 Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1125 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1126 Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1127 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1128
1129
1130 // RMS parameters
1131 Float_t meanRMSY = 0;
1132 Float_t meanRMSZ = 0;
1133 Int_t nclRMS = 0;
1134 for (Int_t idelta = -2; idelta <= 2; idelta++){
1135 // loop over neighbourhood
1136 if (idelta+irow < 0 || idelta+irow > 159) continue;
1137 // if (idelta+irow>159) continue;
1138 AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
1139 if (!cluster) continue;
1140 meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
1141 meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
1142 nclRMS++;
1143 }
1144 meanRMSY /= nclRMS;
1145 meanRMSZ /= nclRMS;
1146
1147 Float_t rmsY = TMath::Sqrt(cluster0->GetSigmaY2());
1148 Float_t rmsZ = TMath::Sqrt(cluster0->GetSigmaZ2());
1149 Float_t rmsYT = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1150 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1151 Float_t rmsZT = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1152 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1153 Float_t rmsYT0 = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1154 TMath::Abs(angley));
1155 Float_t rmsZT0 = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
10757ee9 1156 TMath::Abs(anglez));
9f0beaf7 1157 Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1158 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1159 Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1160 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1161 Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1162 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1163 rmsY,meanRMSY);
1164 Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1165 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1166 rmsZ,meanRMSZ);
1167 //
1168 // cluster debug
1169 TTreeSRedirector *cstream = GetDebugStreamer();
1170 if (cstream){
1171 (*cstream)<<"ResolCl"<< // valgrind 3 40,000 bytes in 1 blocks are possibly
1172 "Sector="<<sector<<
1173 "Cl.="<<cluster0<<
1174 "CSigmaY0="<<csigmaY0<< // cluster errorY
1175 "CSigmaYQ="<<csigmaYQ<< // cluster errorY - q scaled
1176 "CSigmaYS="<<csigmaYS<< // cluster errorY - q scaled
1177 "CSigmaZ0="<<csigmaZ0<< //
1178 "CSigmaZQ="<<csigmaZQ<<
1179 "CSigmaZS="<<csigmaZS<<
1180 "shapeYF="<<rmsYFactor<<
1181 "shapeZF="<<rmsZFactor<<
1182 "rmsY="<<rmsY<<
1183 "rmsZ="<<rmsZ<<
1184 "rmsYM="<<meanRMSY<<
1185 "rmsZM="<<meanRMSZ<<
1186 "rmsYT="<<rmsYT<<
1187 "rmsZT="<<rmsZT<<
1188 "rmsYT0="<<rmsYT0<<
1189 "rmsZT0="<<rmsZT0<<
1190 "rmsYS="<<rmsYSigma<<
1191 "rmsZS="<<rmsZSigma<<
1192 "padSize="<<padSize<<
1193 "Ncl="<<nclFound<<
1194 "PY0.="<<&parY0<<
1195 "PZ0.="<<&parZ0<<
1196 "SigmaY0="<<sigmaY0<<
1197 "SigmaZ0="<<sigmaZ0<<
1198 "angley="<<angley<<
1199 "anglez="<<anglez<<
1200 "\n";
1201
1202 // tracklet dubug
1203 (*cstream)<<"ResolTr"<<
1204 "padSize="<<padSize<<
1205 "IPad="<<padSize<<
1206 "Sector="<<sector<<
1207 "Ncl="<<nclFound<<
1208 "chi2Y0="<<chi2Y0<<
1209 "chi2Z0="<<chi2Z0<<
1210 "chi2Y2="<<chi2Y2<<
1211 "chi2Z2="<<chi2Z2<<
1212 "chi2Y2Q="<<chi2Y2Q<<
1213 "chi2Z2Q="<<chi2Z2Q<<
1214 "chi2Y2S="<<chi2Y2S<<
1215 "chi2Z2S="<<chi2Z2S<<
1216 "PY0.="<<&parY0<<
1217 "PZ0.="<<&parZ0<<
1218 "PY2.="<<&parY2<<
1219 "PZ2.="<<&parZ2<<
1220 "PY2Q.="<<&parY2Q<<
1221 "PZ2Q.="<<&parZ2Q<<
1222 "PY2S.="<<&parY2S<<
1223 "PZ2S.="<<&parZ2S<<
1224 "SigmaY0="<<sigmaY0<<
1225 "SigmaZ0="<<sigmaZ0<<
1226 "SigmaDY0="<<sigmaDY0<<
1227 "SigmaDZ0="<<sigmaDZ0<<
1228 "SigmaY2="<<sigmaY2<<
1229 "SigmaZ2="<<sigmaZ2<<
1230 "SigmaDY2="<<sigmaDY2<<
1231 "SigmaDZ2="<<sigmaDZ2<<
1232 "SigmaY2Q="<<sigmaY2Q<<
1233 "SigmaZ2Q="<<sigmaZ2Q<<
1234 "SigmaDY2Q="<<sigmaDY2Q<<
1235 "SigmaDZ2Q="<<sigmaDZ2Q<<
1236 "SigmaY2S="<<sigmaY2S<<
1237 "SigmaZ2S="<<sigmaZ2S<<
1238 "SigmaDY2S="<<sigmaDY2S<<
1239 "SigmaDZ2S="<<sigmaDZ2S<<
1240 "angley="<<angley<<
1241 "anglez="<<anglez<<
1242 "\n";
1243 }
10757ee9 1244}
1245
1246
1247
1248
1249
1250TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
1251 //
1252 // creates a new histogram which contains the difference between
1253 // the histogram hfit and the function func
1254 //
1255 TH2D * result = (TH2D*)hfit->Clone(); // valgrind 3 40,139 bytes in 11 blocks are still reachable
1256 result->SetTitle(Form("%s fit residuals",result->GetTitle()));
1257 result->SetName(Form("%s fit residuals",result->GetName()));
1258 TAxis *xaxis = hfit->GetXaxis();
1259 TAxis *yaxis = hfit->GetYaxis();
1260 Double_t x[2];
1261 for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
1262 x[1] = yaxis->GetBinCenter(biny);
1263 for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
1264 x[0] = xaxis->GetBinCenter(binx);
1265 Int_t bin = hfit->GetBin(binx, biny);
1266 Double_t val = hfit->GetBinContent(bin);
1267// result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
1268 result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
1269 }
1270 }
1271 return result;
1272}
1273
1274
1275void AliTPCcalibTracks::SetStyle() const {
1276 //
1277 // set style, can be called by all draw functions
1278 //
1279 gROOT->SetStyle("Plain");
1280 gStyle->SetFillColor(10);
1281 gStyle->SetPadColor(10);
1282 gStyle->SetCanvasColor(10);
1283 gStyle->SetStatColor(10);
1284 gStyle->SetPalette(1,0);
1285 gStyle->SetNumberContours(60);
1286}
1287
1288
1289void AliTPCcalibTracks::Draw(Option_t* opt){
1290 //
1291 // draw-function of AliTPCcalibTracks
1292 // will draws some exemplaric pictures
1293 //
1294
ae28e92e 1295 if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
10757ee9 1296 SetStyle();
1297 Double_t min = 0;
1298 Double_t max = 0;
1299 TCanvas *c1 = new TCanvas();
1300 c1->Divide(0, 3);
1301 TVirtualPad *upperThird = c1->GetPad(1);
1302 TVirtualPad *middleThird = c1->GetPad(2);
1303 TVirtualPad *lowerThird = c1->GetPad(3);
1304 upperThird->Divide(2,0);
1305 TVirtualPad *upleft = upperThird->GetPad(1);
1306 TVirtualPad *upright = upperThird->GetPad(2);
1307 middleThird->Divide(2,0);
1308 TVirtualPad *middleLeft = middleThird->GetPad(1);
1309 TVirtualPad *middleRight = middleThird->GetPad(2);
1310 lowerThird->Divide(2,0);
1311 TVirtualPad *downLeft = lowerThird->GetPad(1);
1312 TVirtualPad *downRight = lowerThird->GetPad(2);
1313
1314
1315 upleft->cd(0);
1316 min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1317 max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1318 fDeltaY->SetAxisRange(min, max);
1319 fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
1320 c1->Update();
1321
1322 upright->cd(0);
1323 max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1324 min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1325 fDeltaZ->SetAxisRange(min, max);
1326 fDeltaZ->Fit("gaus","q","",min, max);
1327 c1->Update();
1328
1329 middleLeft->cd();
1330 fHclus->Draw(opt);
1331
1332 middleRight->cd();
1333 fRejectedTracksHisto->Draw(opt);
1334 TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
1335 TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
1336 TText *t2 = pt->AddText("2: kMinClusters");
1337 TText *t3 = pt->AddText("3: kMinRatio");
1338 TText *t4 = pt->AddText("4: kMax1pt");
1339 t1 = t1; t2 = t2; t3 = t3; t4 = t4; // avoid compiler warnings
1340 pt->SetToolTipText("Legend for failed cuts");
1341 pt->Draw();
1342
1343 downLeft->cd();
1344 fHclusterPerPadrowRaw->Draw(opt);
1345
1346 downRight->cd();
1347 fHclusterPerPadrow->Draw(opt);
1348}
1349
1350
1351void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){
1352 //
1353 // all functions are called, that produce pictures
1354 // the histograms are written to the directory 'pathName'
1355 // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
1356 // 'stat' is also the number of minEntries for MakeResPlotsQTree
1357 //
1358
ae28e92e 1359 if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
10757ee9 1360 MakeAmpPlots(stat, pathName);
1361 MakeDeltaPlots(pathName);
1362 FitResolutionNew(pathName);
1363 FitRMSNew(pathName);
1364 MakeChargeVsDriftLengthPlots(pathName);
1365// MakeResPlotsQ(1, 1);
1366 MakeResPlotsQTree(stat, pathName);
1367}
1368
1369
1370void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){
1371 //
1372 // creates several plots:
1373 // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
1374 // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
1375 // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
1376 // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1377 // Empty histograms (sectors without data) are not written to file
1378 // the ps-files are written to the directory 'pathName', that is created if it does not exist
1379 // 'stat': only histograms with more than 'stat' entries are written to file.
1380 //
1381
1382 SetStyle();
1383 gSystem->MakeDirectory(pathName);
1384 gSystem->ChangeDirectory(pathName);
1385
1386 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1387 TPostScript *ps;
1388 // histograms with accumulated amplitude for all IROCs and OROCs
1389 TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
1390 allAmpHisIROC->SetName("Amp all IROCs");
1391 allAmpHisIROC->SetTitle("Amp all IROCs");
1392 TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
1393 allAmpHisOROC->SetName("Amp all OROCs");
1394 allAmpHisOROC->SetTitle("Amp all OROCs");
1395
1396 ps = new TPostScript("fArrayAmp.ps", 112);
ae28e92e 1397 if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
10757ee9 1398 for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
1399 if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat ) continue;
1400 ps->NewPage();
1401 ((TH1F*)fArrayAmp->At(i))->Draw();
1402 c1->Update(); // valgrind 3
1403 if (i > 0 && i < 36) {
1404 allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
1405 allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
1406 }
1407 }
1408 ps->NewPage();
1409 allAmpHisIROC->Draw();
1410 c1->Update(); // valgrind
1411 ps->NewPage();
1412 allAmpHisOROC->Draw();
1413 c1->Update();
1414 ps->Close();
1415 delete ps;
1416
1417 TH1F *his = 0;
1418 Double_t min = 0;
1419 Double_t max = 0;
1420 ps = new TPostScript("fArrayAmpRow.ps", 112);
ae28e92e 1421 if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
10757ee9 1422 for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
1423 his = (TH1F*)fArrayAmpRow->At(i);
1424 if (his->GetEntries() < stat) continue;
1425 ps->NewPage();
1426 min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
1427 max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
1428 his->SetAxisRange(min, max);
1429 his->Fit("pol3", "q", "", min, max);
1430 // his->Draw("error"); // don't use this line when you don't want to have empty pages in the ps-file
1431 c1->Update();
1432 }
1433 ps->Close();
1434 delete ps;
1435 delete c1;
1436 gSystem->ChangeDirectory("..");
1437}
1438
1439
1440void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){
1441 //
1442 // creates several plots:
1443 // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1444 // the ps-files are written to the directory 'pathName', that is created if it does not exist
1445 //
1446
1447 SetStyle();
1448 gSystem->MakeDirectory(pathName);
1449 gSystem->ChangeDirectory(pathName);
1450
1451 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1452 TPostScript *ps;
1453 Double_t min = 0;
1454 Double_t max = 0;
1455
1456 ps = new TPostScript("DeltaYZ.ps", 112);
ae28e92e 1457 if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
10757ee9 1458 min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1459 max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1460 fDeltaY->SetAxisRange(min, max);
1461 ps->NewPage();
1462 fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
1463 c1->Update();
1464 ps->NewPage();
1465 max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1466 min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1467 fDeltaZ->SetAxisRange(min, max);
1468 fDeltaZ->Fit("gaus","q","",min, max);
1469 c1->Update();
1470 ps->Close();
1471 delete ps;
1472 delete c1;
1473 gSystem->ChangeDirectory("..");
1474}
1475
1476
1477void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){
1478 //
1479 // creates charge vs. driftlength plots, one TProfile for each ROC
1480 // is not correct like this, should be one TProfile for each sector and padsize
1481 //
1482
1483 SetStyle();
1484 gSystem->MakeDirectory(pathName);
1485 gSystem->ChangeDirectory(pathName);
1486
1487 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1488 TPostScript *ps;
1489 ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
ae28e92e 1490 if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
10757ee9 1491 TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
1492 TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
1493 chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
1494 chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
1495 chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
1496 chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
1497
1498 for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
1499 ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
1500 c1->Update();
1501 if (i > 0 && i < 36) {
1502 chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
1503 chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
1504 }
1505 ps->NewPage();
1506 }
1507 chargeVsDriftlengthAllIROCs->Draw();
1508 c1->Update(); // valgrind
1509 ps->NewPage();
1510 chargeVsDriftlengthAllOROCs->Draw();
1511 c1->Update();
1512 ps->Close();
1513 delete ps;
1514 delete c1;
1515 gSystem->ChangeDirectory("..");
1516}
1517
1518
1519void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){
1520 //
1521 // creates charge vs. driftlength plots, one TProfile for each ROC
1522 // under development....
1523 //
1524
1525 SetStyle();
1526 gSystem->MakeDirectory(pathName);
1527 gSystem->ChangeDirectory(pathName);
1528
1529 TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700)); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1530// TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
1531 c1->Divide(0,3);
1532 TPostScript *ps;
1533 ps = new TPostScript("chargeVsDriftlength.ps", 111);
ae28e92e 1534 if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
10757ee9 1535
1536 TProfile *chargeVsDriftlengthAllShortPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
1537 TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
1538 TProfile *chargeVsDriftlengthAllLongPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
1539 chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
1540 chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
1541 chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
1542 chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
1543 chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
1544 chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
1545
1546 for (Int_t i = 0; i < 36; i++) {
1547 c1->cd(1)->SetGridx();
1548 c1->cd(1)->SetGridy();
1549 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
1550 c1->cd(2)->SetGridx();
1551 c1->cd(2)->SetGridy();
1552 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
1553 c1->cd(3)->SetGridx();
1554 c1->cd(3)->SetGridy();
1555 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
1556 c1->Update();
1557 chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
1558 chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
1559 chargeVsDriftlengthAllLongPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
1560 ps->NewPage();
1561 }
1562 c1->cd(1)->SetGridx();
1563 c1->cd(1)->SetGridy();
1564 chargeVsDriftlengthAllShortPads->Draw();
1565 c1->cd(2)->SetGridx();
1566 c1->cd(2)->SetGridy();
1567 chargeVsDriftlengthAllMediumPads->Draw();
1568 c1->cd(3)->SetGridx();
1569 c1->cd(3)->SetGridy();
1570 chargeVsDriftlengthAllLongPads->Draw();
1571 c1->Update(); // valgrind
1572// ps->NewPage();
1573 ps->Close();
1574 delete ps;
1575 delete c1;
1576 gSystem->ChangeDirectory("..");
1577}
1578
1579
1580
1581void AliTPCcalibTracks::FitResolutionNew(char* pathName){
1582 //
1583 // calculates different resulution fits in Y and Z direction
1584 // the histograms are written to 'ResolutionYZ.ps'
1585 // writes calculated resolution to 'resol.txt'
1586 // all files are stored in the directory pathName
1587 //
1588
1589 SetStyle();
1590 gSystem->MakeDirectory(pathName);
1591 gSystem->ChangeDirectory(pathName);
1592
1593 TCanvas c;
1594 c.Divide(2,1);
ae28e92e 1595 if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
10757ee9 1596 TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112);
1597 TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1598 fres->SetParameter(0,0.02);
1599 fres->SetParameter(1,0.0054);
1600 fres->SetParameter(2,0.13);
1601
1602 TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
1603
1604 // create histogramw for Y-resolution
1605 TH3F * hisResY0 = (TH3F*)fResolY->At(0);
1606 hisResY0->FitSlicesZ();
1607 TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
1608 TH3F * hisResY1 = (TH3F*)fResolY->At(1);
1609 hisResY1->FitSlicesZ();
1610 TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
1611 TH3F * hisResY2 = (TH3F*)fResolY->At(2);
1612 hisResY2->FitSlicesZ();
1613 TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
1614 //
1615 ps->NewPage();
1616 c.cd(1);
1617 hisResY02->Fit(fres, "q"); // valgrind 132,072 bytes in 6 blocks are indirectly lost
1618 hisResY02->Draw("surf1");
1619 c.cd(2);
1620 MakeDiff(hisResY02,fres)->Draw("surf1");
1621 c.Update();
1622 // c.SaveAs("ResolutionYPad0.eps");
1623 ps->NewPage();
1624 c.cd(1);
1625 hisResY12->Fit(fres, "q");
1626 hisResY12->Draw("surf1");
1627 c.cd(2);
1628 MakeDiff(hisResY12,fres)->Draw("surf1");
1629 c.Update();
1630 // c.SaveAs("ResolutionYPad1.eps");
1631 ps->NewPage();
1632 c.cd(1);
1633 hisResY22->Fit(fres, "q");
1634 hisResY22->Draw("surf1");
1635 c.cd(2);
1636 MakeDiff(hisResY22,fres)->Draw("surf1");
1637 c.Update();
1638// c.SaveAs("ResolutionYPad2.eps");
1639
1640 // create histogramw for Z-resolution
1641 TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
1642 hisResZ0->FitSlicesZ();
1643 TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
1644 TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
1645 hisResZ1->FitSlicesZ();
1646 TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
1647 TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
1648 hisResZ2->FitSlicesZ();
1649 TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
1650
1651 ps->NewPage();
1652 c.cd(1);
1653 hisResZ02->Fit(fres, "q");
1654 hisResZ02->Draw("surf1");
1655 c.cd(2);
1656 MakeDiff(hisResZ02,fres)->Draw("surf1");
1657 c.Update();
1658// c.SaveAs("ResolutionZPad0.eps");
1659 ps->NewPage();
1660 c.cd(1);
1661 hisResZ12->Fit(fres, "q");
1662 hisResZ12->Draw("surf1");
1663 c.cd(2);
1664 MakeDiff(hisResZ12,fres)->Draw("surf1");
1665 c.Update();
1666// c.SaveAs("ResolutionZPad1.eps");
1667 ps->NewPage();
1668 c.cd(1);
1669 hisResZ22->Fit(fres, "q");
1670 hisResZ22->Draw("surf1");
1671 c.cd(2);
1672 MakeDiff(hisResZ22,fres)->Draw("surf1");
1673 c.Update();
1674// c.SaveAs("ResolutionZPad2.eps");
1675 ps->Close();
1676 delete ps;
1677
1678 // write calculated resoltuions to 'resol.txt'
1679 ofstream fresol("resol.txt");
1680 fresol<<"Pad 0.75 cm"<<"\n";
1681 hisResY02->Fit(fres, "q"); // valgrind
1682 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1683 hisResZ02->Fit(fres, "q");
1684 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1685 //
1686 fresol<<"Pad 1.00 cm"<<1<<"\n";
1687 hisResY12->Fit(fres, "q"); // valgrind
1688 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1689 hisResZ12->Fit(fres, "q");
1690 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1691 //
1692 fresol<<"Pad 1.50 cm"<<0<<"\n";
1693 hisResY22->Fit(fres, "q");
1694 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1695 hisResZ22->Fit(fres, "q");
1696 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1697
1698 TH1::AddDirectory(kFALSE);
1699 gSystem->ChangeDirectory("..");
1700 delete fres;
1701}
1702
1703
1704void AliTPCcalibTracks::FitRMSNew(char* pathName){
1705 //
1706 // calculates different resulution-rms fits in Y and Z direction
1707 // the histograms are written to 'RMS_YZ.ps'
1708 // writes calculated resolution rms to 'rms.txt'
1709 // all files are stored in the directory pathName
1710 //
1711
1712 SetStyle();
1713 gSystem->MakeDirectory(pathName);
1714 gSystem->ChangeDirectory(pathName);
1715
1716 TCanvas c; // valgrind 3 42,120 bytes in 405 blocks are still reachable 23,816 bytes in 229 blocks are still reachable
1717 c.Divide(2,1);
ae28e92e 1718 if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
10757ee9 1719 TPostScript *ps = new TPostScript("RMS_YZ.ps", 112);
1720 TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1721 frms->SetParameter(0,0.02);
1722 frms->SetParameter(1,0.0054);
1723 frms->SetParameter(2,0.13);
1724
1725 TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
1726
1727 // create histogramw for Y-RMS
1728 TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
1729 hisResY0->FitSlicesZ();
1730 TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
1731 TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
1732 hisResY1->FitSlicesZ();
1733 TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
1734 TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
1735 hisResY2->FitSlicesZ();
1736 TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
1737 //
1738 ps->NewPage();
1739 c.cd(1);
1740 hisResY02->Fit(frms, "qn0");
1741 hisResY02->Draw("surf1");
1742 c.cd(2);
1743 MakeDiff(hisResY02,frms)->Draw("surf1");
1744 c.Update();
1745// c.SaveAs("RMSYPad0.eps");
1746 ps->NewPage();
1747 c.cd(1);
1748 hisResY12->Fit(frms, "qn0"); // valgrind several blocks possibly lost
1749 hisResY12->Draw("surf1");
1750 c.cd(2);
1751 MakeDiff(hisResY12,frms)->Draw("surf1");
1752 c.Update();
1753// c.SaveAs("RMSYPad1.eps");
1754 ps->NewPage();
1755 c.cd(1);
1756 hisResY22->Fit(frms, "qn0");
1757 hisResY22->Draw("surf1");
1758 c.cd(2);
1759 MakeDiff(hisResY22,frms)->Draw("surf1");
1760 c.Update();
1761// c.SaveAs("RMSYPad2.eps");
1762
1763 // create histogramw for Z-RMS
1764 TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
1765 hisResZ0->FitSlicesZ();
1766 TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
1767 TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1);
1768 hisResZ1->FitSlicesZ();
1769 TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
1770 TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2);
1771 hisResZ2->FitSlicesZ();
1772 TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
1773 //
1774 ps->NewPage();
1775 c.cd(1);
1776 hisResZ02->Fit(frms, "qn0"); // valgrind
1777 hisResZ02->Draw("surf1");
1778 c.cd(2);
1779 MakeDiff(hisResZ02,frms)->Draw("surf1");
1780 c.Update();
1781// c.SaveAs("RMSZPad0.eps");
1782 ps->NewPage();
1783 c.cd(1);
1784 hisResZ12->Fit(frms, "qn0");
1785 hisResZ12->Draw("surf1");
1786 c.cd(2);
1787 MakeDiff(hisResZ12,frms)->Draw("surf1");
1788 c.Update();
1789// c.SaveAs("RMSZPad1.eps");
1790 ps->NewPage();
1791 c.cd(1);
1792 hisResZ22->Fit(frms, "qn0"); // valgrind 1 block possibly lost
1793 hisResZ22->Draw("surf1");
1794 c.cd(2);
1795 MakeDiff(hisResZ22,frms)->Draw("surf1");
1796 c.Update();
1797// c.SaveAs("RMSZPad2.eps");
1798
1799 // write calculated resoltuion rms to 'rms.txt'
1800 ofstream filerms("rms.txt");
1801 filerms<<"Pad 0.75 cm"<<"\n";
1802 hisResY02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
1803 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1804 hisResZ02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
1805 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1806 //
1807 filerms<<"Pad 1.00 cm"<<1<<"\n";
1808 hisResY12->Fit(frms, "qn0"); // valgrind 3,256 bytes in 22 blocks are indirectly lost
1809 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1810 hisResZ12->Fit(frms, "qn0"); // valgrind 66,036 bytes in 3 blocks are still reachable
1811 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1812 //
1813 filerms<<"Pad 1.50 cm"<<0<<"\n";
1814 hisResY22->Fit(frms, "qn0"); // valgrind 40,139 bytes in 11 blocks are still reachable 330,180 bytes in 15 blocks are possibly lost
1815 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1816 hisResZ22->Fit(frms, "qn0");
1817 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1818
1819 TH1::AddDirectory(kFALSE);
1820 gSystem->ChangeDirectory("..");
1821 ps->Close();
1822 delete ps;
1823 delete frms;
1824}
1825
1826
1827void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){
1828 //
1829 // Make tree with resolution parameters
1830 // the result is written to 'resol.root' in directory 'pathname'
1831 // file information are available in fileInfo
1832 // available variables in the tree 'Resol':
1833 // Entries: number of entries for this resolution point
1834 // nbins: number of bins that were accumulated
1835 // Dim: direction, Dim==0: y-direction, Dim==1: z-direction
1836 // Pad: padSize; short, medium and long
1837 // Length: pad length, 0.75, 1, 1.5
1838 // QMean: mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1839 // Zc: center of middle bin in drift direction
1840 // Zm: mean dirftlength for accumulated Delta-Histograms
1841 // Zs: width of driftlength bin
1842 // AngleC: center of middle bin in Angle-Direction
1843 // AngleM: mean angle for accumulated Delta-Histograms
1844 // AngleS: width of Angle-bin
1845 // Resol: sigma for gaus fit through Delta-Histograms
1846 // Sigma: error of sigma for gaus fit through Delta Histograms
1847 // MeanR: mean of the Delta-Histogram
1848 // SigmaR: rms of the Delta-Histogram
1849 // RMSm: mean of the gaus fit through RMS-Histogram
1850 // RMS: sigma of the gaus fit through RMS-Histogram
1851 // RMSe0: error of mean of gaus fit in RMS-Histogram
1852 // RMSe1: error of sigma of gaus fit in RMS-Histogram
1853 //
1854
ae28e92e 1855 if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
1856 if (GetDebugLevel() > -1) cout << " relax, the calculation will take a while..." << endl;
10757ee9 1857
1858 gSystem->MakeDirectory(pathName);
1859 gSystem->ChangeDirectory(pathName);
1860 TString kFileName = "resol.root";
1861 TTreeSRedirector fTreeResol(kFileName.Data());
1862
1863 TH3F *resArray[2][3][11];
1864 TH3F *rmsArray[2][3][11];
1865
1866 // load histograms from fArrayQDY and fArrayQDZ
1867 // into resArray and rmsArray
1868 // that is all we need here
1869 for (Int_t idim = 0; idim < 2; idim++){
1870 for (Int_t ipad = 0; ipad < 3; ipad++){
1871 for (Int_t iq = 0; iq <= 10; iq++){
1872 rmsArray[idim][ipad][iq]=0;
1873 resArray[idim][ipad][iq]=0;
1874 Int_t bin = GetBin(iq,ipad);
1875 TH3F *hresl = 0;
1876 if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
1877 if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
1878 if (!hresl) continue;
1879 resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
1880 resArray[idim][ipad][iq]->SetDirectory(0);
1881 TH3F * hreslRMS = 0;
1882 if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
1883 if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
1884 if (!hreslRMS) continue;
1885 rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
1886 rmsArray[idim][ipad][iq]->SetDirectory(0);
1887 }
1888 }
1889 }
ae28e92e 1890 if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
10757ee9 1891
1892 //--------------------------------------------------------------------------------------------
1893
1894 char name[200];
1895 Double_t qMean;
1896 Double_t zMean, angleMean, zCenter, angleCenter;
1897 Double_t zSigma, angleSigma;
1898 TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1899 TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1900 TF1 *fitFunction = new TF1("fitFunction", "gaus");
1901 Float_t entriesQ = 0;
1902 Int_t loopCounter = 1;
1903
1904 for (Int_t idim = 0; idim < 2; idim++){
1905 // Loop y-z corrdinate
1906 for (Int_t ipad = 0; ipad < 3; ipad++){
1907 // loop pad type
1908 for (Int_t iq = -1; iq < 10; iq++){
1909 // LOOP Q
ae28e92e 1910 if (GetDebugLevel() > -1)
10757ee9 1911 cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
1912 << (Int_t)((loopCounter)/66.*100) << "% done), "
1913 << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << " \r" << std::flush;
1914 loopCounter++;
1915 TH3F *hres = 0;
1916 TH3F *hrms = 0;
1917 qMean = 0;
1918 entriesQ = 0;
1919
1920 // calculate qMean
1921 if (iq == -1){
1922 // integrated spectra
1923 for (Int_t iql = 0; iql < 10; iql++){
1924 Int_t bin = GetBin(iql,ipad);
1925 TH3F *hresl = resArray[idim][ipad][iql];
1926 TH3F *hrmsl = rmsArray[idim][ipad][iql];
1927 if (!hresl) continue;
1928 if (!hrmsl) continue;
1929 entriesQ += hresl->GetEntries();
1930 qMean += hresl->GetEntries() * GetQ(bin);
1931 if (!hres) {
1932 hres = (TH3F*)hresl->Clone();
1933 hrms = (TH3F*)hrmsl->Clone();
1934 }
1935 else{
1936 hres->Add(hresl);
1937 hrms->Add(hrmsl);
1938 }
1939 }
1940 qMean /= entriesQ;
1941 qMean *= -1.; // integral mean charge
1942 }
1943 else {
1944 // loop over neighboured Q-bins
1945 // accumulate entries from neighboured Q-bins
1946 for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
1947 if (iql < 0) continue;
1948 Int_t bin = GetBin(iql,ipad);
1949 TH3F * hresl = resArray[idim][ipad][iql];
1950 TH3F * hrmsl = rmsArray[idim][ipad][iql];
1951 if (!hresl) continue;
1952 if (!hrmsl) continue;
1953 entriesQ += hresl->GetEntries();
1954 qMean += hresl->GetEntries() * GetQ(bin);
1955 if (!hres) {
1956 hres = (TH3F*) hresl->Clone();
1957 hrms = (TH3F*) hrmsl->Clone();
1958 }
1959 else{
1960 hres->Add(hresl);
1961 hrms->Add(hrmsl);
1962 }
1963 }
1964 qMean/=entriesQ;
1965 }
1966
1967 TAxis *xAxisDriftLength = hres->GetXaxis(); // driftlength / z - axis
1968 TAxis *yAxisAngle = hres->GetYaxis(); // angle axis
1969 TAxis *zAxisDelta = hres->GetZaxis(); // delta axis
1970 TAxis *zAxisRms = hrms->GetZaxis(); // rms axis
1971
1972 // loop over all angle bins
1973 for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1974 angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1975 // loop over all driftlength bins
1976 for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1977 zCenter = xAxisDriftLength->GetBinCenter(ibinxDL);
1978 zSigma = xAxisDriftLength->GetBinWidth(ibinxDL);
1979 angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
1980 zMean = zCenter; // changens, when more statistic is accumulated
1981 angleMean = angleCenter; // changens, when more statistic is accumulated
1982
1983 // create 2 1D-Histograms, projectionRes and projectionRms
1984 // these histograms are delta histograms for given direction, padSize, chargeBin,
1985 // angleBin and driftLengthBin
1986 // later on they will be fitted with a gausian, its sigma is the resoltuion...
1987 sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
1988 // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1989 projectionRes->SetNameTitle(name, name);
1990 sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
1991 // TH1D * projectionRms = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1992 projectionRms->SetNameTitle(name, name);
1993
1994 projectionRes->Reset();
1995 projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1996 projectionRms->Reset();
1997 projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1998 projectionRes->SetDirectory(0);
1999 projectionRms->SetDirectory(0);
2000
2001 Double_t entries = 0;
2002 Int_t nbins = 0; // counts, how many bins were accumulated
2003 zMean = 0;
2004 angleMean = 0;
2005 entries = 0;
2006
2007 // fill projectionRes and projectionRms for given dim, ipad and iq,
2008 // as well as for given angleBin and driftlengthBin
2009 // if this gives not enough statistic, include neighbourhood
2010 // (angle and driftlength) successifely
2011 for (Int_t dbin = 0; dbin <= 8; dbin++){ // delta-bins around centered angleBin and driftlengthBin
2012 for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) { // delta-bins in angle direction
2013 for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
2014 if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; // add each bin only one time !
2015 Int_t binx2 = ibinxDL + dbinx2; // position variable in x (driftlength) direction
2016 Int_t biny2 = ibinyAngle + dbiny2; // position variable in y (angle) direction
2017 if (binx2 < 1 || biny2 < 1) continue; // don't go out of the histogram!
2018 if (binx2 >= xAxisDriftLength->GetNbins()) continue; // don't go out of the histogram!
2019 if (biny2 >= yAxisAngle->GetNbins()) continue; // don't go out of the histogram!
2020 nbins++; // count the number of accumulated bins
2021 // Fill resolution histo
2022 for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
2023 // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3); // unused variable
2024 projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
2025 entries += hres->GetBinContent(binx2, biny2, ibin3);
2026 zMean += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
2027 angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
2028 } // ibin3 loop
2029 // fill RMS histo
2030 for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
2031 projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
2032 }
2033 } //dbinx2 loop
2034 if (entries > minEntries) break; // enough statistic accumulated
2035 } // dbiny2 loop
2036 if (entries > minEntries) break; // enough statistic accumulated
2037 } // dbin loop
2038 if ( entries< minEntries) continue; // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
2039 zMean /= entries;
2040 angleMean /= entries;
2041
2042 if (entries > minEntries) {
2043 // when enough statistic is accumulated
2044 // fit Delta histograms with a gausian
2045 // of the gausian is the resolution (resol), its fit error is sigma
2046 // store also mean and RMS of the histogram
2047 Float_t xmin = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2048 Float_t xmax = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2049
2050// projectionRes->Fit("gaus", "q0", "", xmin, xmax);
2051// Float_t resol = projectionRes->GetFunction("gaus")->GetParameter(2);
2052// Float_t sigma = projectionRes->GetFunction("gaus")->GetParError(2);
2053 fitFunction->SetMaximum(xmax);
2054 fitFunction->SetMinimum(xmin);
2055 projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
2056 Float_t resol = fitFunction->GetParameter(2);
2057 Float_t sigma = fitFunction->GetParError(2);
2058
2059 Float_t meanR = projectionRes->GetMean();
2060 Float_t sigmaR = projectionRes->GetRMS();
2061 // fit also RMS histograms with a gausian
2062 // store mean and sigma of the gausian in rmsMean and rmsSigma
2063 // store also the fit errors in errorRMS and errorSigma
2064 xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2065 xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2066
2067// projectionRms->Fit("gaus","q0","",xmin,xmax);
2068// Float_t rmsMean = projectionRms->GetFunction("gaus")->GetParameter(1);
2069// Float_t rmsSigma = projectionRms->GetFunction("gaus")->GetParameter(2);
2070// Float_t errorRMS = projectionRms->GetFunction("gaus")->GetParError(1);
2071// Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
2072 projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
2073 Float_t rmsMean = fitFunction->GetParameter(1);
2074 Float_t rmsSigma = fitFunction->GetParameter(2);
2075 Float_t errorRMS = fitFunction->GetParError(1);
2076 Float_t errorSigma = fitFunction->GetParError(2);
2077
2078 Float_t length = 0.75;
2079 if (ipad == 1) length = 1;
2080 if (ipad == 2) length = 1.5;
2081
2082 fTreeResol<<"Resol"<<
2083 "Entries="<<entries<< // number of entries for this resolution point
2084 "nbins="<<nbins<< // number of bins that were accumulated
2085 "Dim="<<idim<< // direction, Dim==0: y-direction, Dim==1: z-direction
2086 "Pad="<<ipad<< // padSize; short, medium and long
2087 "Length="<<length<< // pad length, 0.75, 1, 1.5
2088 "QMean="<<qMean<< // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
2089 "Zc="<<zCenter<< // center of middle bin in drift direction
2090 "Zm="<<zMean<< // mean dirftlength for accumulated Delta-Histograms
2091 "Zs="<<zSigma<< // width of driftlength bin
2092 "AngleC="<<angleCenter<< // center of middle bin in Angle-Direction
2093 "AngleM="<<angleMean<< // mean angle for accumulated Delta-Histograms
2094 "AngleS="<<angleSigma<< // width of Angle-bin
2095 "Resol="<<resol<< // sigma for gaus fit through Delta-Histograms
2096 "Sigma="<<sigma<< // error of sigma for gaus fit through Delta Histograms
2097 "MeanR="<<meanR<< // mean of the Delta-Histogram
2098 "SigmaR="<<sigmaR<< // rms of the Delta-Histogram
2099 "RMSm="<<rmsMean<< // mean of the gaus fit through RMS-Histogram
2100 "RMSs="<<rmsSigma<< // sigma of the gaus fit through RMS-Histogram
2101 "RMSe0="<<errorRMS<< // error of mean of gaus fit in RMS-Histogram
2102 "RMSe1="<<errorSigma<< // error of sigma of gaus fit in RMS-Histogram
2103 "\n";
ae28e92e 2104 if (GetDebugLevel() > 5) {
10757ee9 2105 projectionRes->SetDirectory(fTreeResol.GetFile());
2106 projectionRes->Write(projectionRes->GetName());
2107 projectionRes->SetDirectory(0);
2108 projectionRms->SetDirectory(fTreeResol.GetFile());
2109 projectionRms->Write(projectionRms->GetName());
2110 projectionRes->SetDirectory(0);
2111 }
2112 } // if (projectionRes->GetSum() > minEntries)
2113 } // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
2114 } // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
2115
2116 } // iq-loop
2117 } // ipad-loop
2118 } // idim-loop
2119 delete projectionRes;
2120 delete projectionRms;
2121
2122// TFile resolFile(fTreeResol.GetFile());
2123 TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
2124 fileInfo.Write("fileInfo");
2125// resolFile.Close();
2126// fTreeResol.GetFile()->Close();
ae28e92e 2127 if (GetDebugLevel() > -1) cout << endl;
2128 if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
10757ee9 2129 gSystem->ChangeDirectory("..");
2130}
2131
2132
10757ee9 2133
10757ee9 2134
2135
2136Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
2137 //
2138 // function to merge several AliTPCcalibTracks objects after PROOF calculation
2139 // The object's histograms are merged via their merge functions
2140 // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
2141 //
2142
ae28e92e 2143 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl;
10757ee9 2144 if (!collectionList) return 0;
2145 if (collectionList->IsEmpty()) return -1;
2146
ae28e92e 2147 if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
2148 if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
10757ee9 2149 collectionList->Print();
2150
2151 // create a list for each data member
2152 TList* deltaYList = new TList;
2153 TList* deltaZList = new TList;
2154 TList* arrayAmpRowList = new TList;
2155 TList* rejectedTracksList = new TList;
2156 TList* hclusList = new TList;
2157 TList* clusterCutHistoList = new TList;
2158 TList* arrayAmpList = new TList;
2159 TList* arrayQDYList = new TList;
2160 TList* arrayQDZList = new TList;
2161 TList* arrayQRMSYList = new TList;
2162 TList* arrayQRMSZList = new TList;
2163 TList* arrayChargeVsDriftlengthList = new TList;
2164 TList* calPadRegionChargeVsDriftlengthList = new TList;
2165 TList* hclusterPerPadrowList = new TList;
2166 TList* hclusterPerPadrowRawList = new TList;
2167 TList* resolYList = new TList;
2168 TList* resolZList = new TList;
2169 TList* rMSYList = new TList;
2170 TList* rMSZList = new TList;
2171
2172// TList* nRowsList = new TList;
2173// TList* nSectList = new TList;
2174// TList* fileNoList = new TList;
2175
2176 TIterator *listIterator = collectionList->MakeIterator();
2177 AliTPCcalibTracks *calibTracks = 0;
ae28e92e 2178 if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
10757ee9 2179 Int_t counter = 0;
2180 while ( (calibTracks = (AliTPCcalibTracks*)listIterator->Next()) ){
2181 // loop over all entries in the collectionList and get dataMembers into lists
2182 if (!calibTracks) continue;
2183
2184 deltaYList->Add( calibTracks->GetfDeltaY() );
2185 deltaZList->Add( calibTracks->GetfDeltaZ() );
2186 arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
2187 arrayAmpList->Add(calibTracks->GetfArrayAmp());
2188 arrayQDYList->Add(calibTracks->GetfArrayQDY());
2189 arrayQDZList->Add(calibTracks->GetfArrayQDZ());
2190 arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
2191 arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
2192 resolYList->Add(calibTracks->GetfResolY());
2193 resolZList->Add(calibTracks->GetfResolZ());
2194 rMSYList->Add(calibTracks->GetfRMSY());
2195 rMSZList->Add(calibTracks->GetfRMSZ());
2196 arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
2197 calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
2198 hclusList->Add(calibTracks->GetfHclus());
2199 rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
2200 clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
2201 hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
2202 hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
2203 fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
2204 fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
2205 counter++;
ae28e92e 2206 if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
10757ee9 2207 }
2208
10757ee9 2209
2210 // merge data members
ae28e92e 2211 if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
10757ee9 2212 fDeltaY->Merge(deltaYList);
2213 fDeltaZ->Merge(deltaZList);
2214 fHclus->Merge(hclusList);
2215 fClusterCutHisto->Merge(clusterCutHistoList);
2216 fRejectedTracksHisto->Merge(rejectedTracksList);
2217 fHclusterPerPadrow->Merge(hclusterPerPadrowList);
2218 fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
2219
2220 TObjArray* objarray = 0;
2221 TH1* hist = 0;
2222 TList* histList = 0;
2223 TIterator *objListIterator = 0;
2224
ae28e92e 2225 if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
10757ee9 2226 // merge fArrayAmpRows
2227 for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) { // loop over data member, i<72
2228 objListIterator = arrayAmpRowList->MakeIterator();
2229 histList = new TList;
2230 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2231 // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
2232 if (!objarray) continue;
2233 hist = (TProfile*)(objarray->At(i));
2234 histList->Add(hist);
2235 }
2236 ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
2237 delete histList;
2238 delete objListIterator;
2239 }
2240
ae28e92e 2241 if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
10757ee9 2242 // merge fArrayAmps
2243 for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) { // loop over data member, i<72
2244 TIterator *objListIterator = arrayAmpList->MakeIterator();
2245 histList = new TList;
2246 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2247 // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
2248 if (!objarray) continue;
2249 hist = (TH1F*)(objarray->At(i));
2250 histList->Add(hist);
2251 }
2252 ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
2253 delete histList;
2254 delete objListIterator;
2255 }
2256
ae28e92e 2257 if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
10757ee9 2258 // merge fArrayQDY
2259 for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
2260 objListIterator = arrayQDYList->MakeIterator();
2261 histList = new TList;
2262 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2263 // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
2264 if (!objarray) continue;
2265 hist = (TH3F*)(objarray->At(i));
2266 histList->Add(hist);
2267 }
2268 ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
2269 delete histList;
2270 delete objListIterator;
2271 }
2272
ae28e92e 2273 if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
10757ee9 2274 // merge fArrayQDZ
2275 for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2276 objListIterator = arrayQDZList->MakeIterator();
2277 histList = new TList;
2278 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2279 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2280 if (!objarray) continue;
2281 hist = (TH3F*)(objarray->At(i));
2282 histList->Add(hist);
2283 }
2284 ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
2285 delete histList;
2286 delete objListIterator;
2287 }
2288
ae28e92e 2289 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
10757ee9 2290 // merge fArrayQRMSY
2291 for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
2292 objListIterator = arrayQRMSYList->MakeIterator();
2293 histList = new TList;
2294 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2295 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2296 if (!objarray) continue;
2297 hist = (TH3F*)(objarray->At(i));
2298 histList->Add(hist);
2299 }
2300 ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
2301 delete histList;
2302 delete objListIterator;
2303 }
2304
ae28e92e 2305 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
10757ee9 2306 // merge fArrayQRMSZ
2307 for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2308 objListIterator = arrayQRMSZList->MakeIterator();
2309 histList = new TList;
2310 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2311 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2312 if (!objarray) continue;
2313 hist = (TH3F*)(objarray->At(i));
2314 histList->Add(hist);
2315 }
2316 ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
2317 delete histList;
2318 delete objListIterator;
2319 }
2320
ae28e92e 2321 if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
10757ee9 2322 // merge fArrayChargeVsDriftlength
2323 for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
2324 objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
2325 histList = new TList;
2326 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2327 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
2328 if (!objarray) continue;
2329 hist = (TProfile*)(objarray->At(i));
2330 histList->Add(hist);
2331 }
2332 ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
2333 delete histList;
2334 delete objListIterator;
2335 }
2336
ae28e92e 2337 if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
10757ee9 2338 // merge fcalPadRegionChargeVsDriftlength
2339 AliTPCCalPadRegion *cpr = 0x0;
2340
2341 /*
2342 TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
2343 while (hist = (TProfile*)regionIterator->Next()) {
2344 // loop over all calPadRegion's in destination calibTracks object
2345 objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2346 while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
2347
2348
2349 hist->Merge(...);
2350 }
2351 */
2352
2353 for (UInt_t isec = 0; isec < 36; isec++) {
2354 for (UInt_t padSize = 0; padSize < 3; padSize++){
2355 objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2356 histList = new TList;
2357 while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
2358 // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object
2359 if (!cpr) continue;
2360 hist = (TProfile*)cpr->GetObject(isec, padSize);
2361 histList->Add(hist);
2362 }
2363 ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
2364 delete histList;
2365 delete objListIterator;
2366 }
2367 }
2368
2369
2370
2371
ae28e92e 2372 if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
10757ee9 2373 // merge fResolY
2374 for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
2375 objListIterator = resolYList->MakeIterator();
2376 histList = new TList;
2377 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2378 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2379 if (!objarray) continue;
2380 hist = (TH3F*)(objarray->At(i));
2381 histList->Add(hist);
2382 }
2383 ((TH3F*)(fResolY->At(i)))->Merge(histList);
2384 delete histList;
2385 delete objListIterator;
2386 }
2387
2388 // merge fResolZ
2389 for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2390 objListIterator = resolZList->MakeIterator();
2391 histList = new TList;
2392 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2393 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2394 if (!objarray) continue;
2395 hist = (TH3F*)(objarray->At(i));
2396 histList->Add(hist);
2397 }
2398 ((TH3F*)(fResolZ->At(i)))->Merge(histList);
2399 delete histList;
2400 delete objListIterator;
2401 }
2402
2403 // merge fRMSY
2404 for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
2405 objListIterator = rMSYList->MakeIterator();
2406 histList = new TList;
2407 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2408 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2409 if (!objarray) continue;
2410 hist = (TH3F*)(objarray->At(i));
2411 histList->Add(hist);
2412 }
2413 ((TH3F*)(fRMSY->At(i)))->Merge(histList);
2414 delete histList;
2415 delete objListIterator;
2416 }
2417
2418 // merge fRMSZ
2419 for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2420 objListIterator = rMSZList->MakeIterator();
2421 histList = new TList;
2422 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2423 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2424 if (!objarray) continue;
2425 hist = (TH3F*)(objarray->At(i));
2426 histList->Add(hist);
2427 }
2428 ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
2429 delete histList;
2430 delete objListIterator;
2431 }
2432
2433 delete deltaYList;
2434 delete deltaZList;
2435 delete arrayAmpRowList;
2436 delete arrayAmpList;
2437 delete arrayQDYList;
2438 delete arrayQDZList;
2439 delete arrayQRMSYList;
2440 delete arrayQRMSZList;
2441 delete resolYList;
2442 delete resolZList;
2443 delete rMSYList;
2444 delete rMSZList;
10757ee9 2445 delete listIterator;
2446
ae28e92e 2447 if (GetDebugLevel() > 0) cout << "merging done!" << endl;
10757ee9 2448
2449 return 1;
2450}
2451
2452
2453AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
2454 //
2455 // function to test AliTPCcalibTrack::Merge:
2456 // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
2457 // this object is appended 'nCalTracks' times to a TList
2458 // A new AliTPCcalibTrack object is created which merges the list
2459 // this object is returned
2460 //
2461 /*
2462 // .L AliTPCcalibTracks.cxx+g
2463 .L libTPCcalib.so
2464 TFile f("Output.root");
2465 AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
2466 //f.Close();
2467 TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
2468 AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param");
2469 clusterParamFile.Close();
2470
2471 AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
2472 */
2473 TList *list = new TList();
2474 if (ct == 0 || clusterParam == 0) return 0;
2475 cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
2476 for (Int_t i = 0; i < nCalTracks; i++) {
2477 if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
b42cf9ed 2478 list->Add(new AliTPCcalibTracks(*ct));
10757ee9 2479 }
2480
2481 // only for check at the end
b42cf9ed 2482 AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
10757ee9 2483 Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
2484// Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
2485
2486 cout << "The list contains " << list->GetEntries() << " entries. " << endl;
2487
2488
2489 AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
2490 AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
2491 cal->Merge(list);
2492
2493 cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
2494 cal->GetfArrayAmpRow()->At(5)->Print();
2495 Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
2496
2497 cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
2498 cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " << calEntries << endl;
2499 printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n",
2500 calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
2501
2502 return cal;
2503
2504}
2505
2506
9c171b96 2507void AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints, Int_t verbose){
2508 //
2509 // Make position corrections
2510 // for the moment Only using debug streamer
2511 // chainres - debug tree
2512 // param - parameters to be updated
2513 // maxPoints - maximal number of points using for fit
2514 // verbose - print info flag
2e5bcb67 2515 //
2516 // Current implementation - only using debug streamers
2517 //
2518
2519 /*
9c171b96 2520 //Defaults
2521 Int_t maxPoints=100000;
2522 */
2e5bcb67 2523 /*
2524 Usage:
2525 //0. Load libraries
2526 gSystem->Load("libANALYSIS");
2527 gSystem->Load("libSTAT");
2528 gSystem->Load("libTPCcalib");
2529
2530
2531 //1. Load Parameters to be modified:
2532 //e.g:
2533
2534 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
2535 AliCDBManager::Instance()->SetRun(0)
2536 AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
2537
2538 //2. Load chain from debug streamers
2539 //
2540 //e.g
2541 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2542 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2543 AliXRDPROOFtoolkit tool;
2544 TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2545 chainres->Lookup();
2546 //3. Do fits and store results
2547 //
2548 AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
2549 TFile f("paramout.root","recreate");
2550 param->Write("clusterParam");
2551 f.Close();
2552 //4. Verification
2553 TFile f2("paramout.root");
2554 AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
2555 param2->SetInstance(param2);
2556 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
2557
2558 */
2559
2560
2561 TStatToolkit toolkit;
2562 Double_t chi2;
2563 TVectorD fitParamY0;
2564 TVectorD fitParamY1;
2565 TVectorD fitParamZ0;
2566 TVectorD fitParamZ1;
2567 TMatrixD covMatrix;
2568 Int_t npoints;
2569
2570 chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2571 chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2572 chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2573 chainres->SetAlias("st","(sin(dt)-dt)");
2574 //
2575 chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2576 chainres->SetAlias("dq","sqrt(15./(5+Cl.fMax))");
2577 chainres->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
2578 chainres->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
2579 //
2580 //
2581 //
2582 TCut cutA("1");
2583 TString fstringY="";
2584 //
2585 fstringY+="(dp)++"; //1
2586 fstringY+="(dp)*di++"; //2
2587 fstringY+="(dp)*dq++"; //3
2588 fstringY+="(dp)*sy++"; //4
2589 //
2590 fstringY+="(sp)++"; //5
2591 fstringY+="(sp)*di++"; //6
2592 fstringY+="(sp)*dq++"; //7
2593 fstringY+="(sp)*sy++"; //8
2594 //
2595
2596 TString fstringZ="";
2597 fstringZ+="(dt)++"; //1
2598 fstringZ+="(dt)*di++"; //2
2599 fstringZ+="(dt)*dq++"; //3
2600 fstringZ+="(dt)*sz++"; //4
2601 //
2602 fstringZ+="(st)++"; //5
2603 fstringZ+="(st)*di++"; //6
2604 fstringZ+="(st)*dq++"; //7
2605 fstringZ+="(st)*sz++"; //8
2606 //
2607 // Z corrections
2608 //
2609 TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
2610 printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2611 param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
2612 //
2613 TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
2614 //
2615 printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2616 param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
2617 param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
2618 //
2619 // Y corrections
2620 //
2621 TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
2622 printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2623 param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
2624
2625
2626 TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
2627 //
2628 printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2629 param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
2630 param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
2631 //
2632 //
2633 //
2634 chainres->SetAlias("fitZ0",strZ0->Data());
2635 chainres->SetAlias("fitZ1",strZ1->Data());
2636 chainres->SetAlias("fitY0",strY0->Data());
2637 chainres->SetAlias("fitY1",strY1->Data());
2638 // chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);
9c171b96 2639}
2640
2641/*
2e5bcb67 2642 gSystem->Load("libSTAT.so");
9c171b96 2643 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2644 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2645 AliXRDPROOFtoolkit tool;
2646 TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2647 chainres->Lookup();
2648 //
2649 //
2650 //
2651
2652
9c171b96 2653
2654 TStatToolkit toolkit;
2655 Double_t chi2;
2656 TVectorD fitParamY0;
2657 TVectorD fitParamY1;
2658 TVectorD fitParamZ0;
2659 TVectorD fitParamZ1;
2660 TMatrixD covMatrix;
2661 Int_t npoints;
2662
2663 chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2664 chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2665 chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2666 chainres->SetAlias("st","(sin(dt)-dt)");
2667 //
2668 chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2669 chainres->SetAlias("dq","sqrt(15./(5+Cl.fMax))");
2670 chainres->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
2671 chainres->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
2672 //
2673 //
2674 //
2675 TCut cutA("1")
2676
2677 TString fstringY="";
2678 //
2679 fstringY+="(dp)++"; //1
2680 fstringY+="(dp)*di++"; //2
2681 fstringY+="(dp)*dq++"; //3
2682 fstringY+="(dp)*sy++"; //4
2683 //
2684 fstringY+="(sp)++"; //5
2685 fstringY+="(sp)*di++"; //6
2686 fstringY+="(sp)*dq++"; //7
2687 fstringY+="(sp)*sy++"; //8
2688 //
2689
2690 TString fstringZ="";
2691 fstringZ+="(dt)++"; //1
2692 fstringZ+="(dt)*di++"; //2
2693 fstringZ+="(dt)*dq++"; //3
2694 fstringZ+="(dt)*sz++"; //4
2695 //
2696 fstringZ+="(st)++"; //5
2697 fstringZ+="(st)*di++"; //6
2698 fstringZ+="(st)*dq++"; //7
2699 fstringZ+="(st)*sz++"; //8
2700
2701
2e5bcb67 2702 TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,100000);
2703 TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,100000);
9c171b96 2704
2e5bcb67 2705 TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,100000);
2706 TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,100000);
9c171b96 2707
2e5bcb67 2708 chainres->SetAlias("fitZ0",strZ0.Data());
2709 chainres->SetAlias("fitZ1",strZ1.Data());
2710 chainres->SetAlias("fitY0",strY0.Data());
2711 chainres->SetAlias("fitY1",strY1.Data());
9c171b96 2712 chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000)
2713
2714
2715
2716*/
2717
10757ee9 2718
2719