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