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