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