]> 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
Getter for subcorrections
[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// //
60e31b96 19// Class for calibration of the cluster properties: //
20// Cluster position resolution (RMS) and short term distortions (within pad or within time bin)
21//
22// Algorithm:
23// 1. Creation of the residual/properties histograms
24// 2. Filling of the histograms.
25// 2.a Traklet fitting
26// 2.b Resudual filling
27// 3. Postprocessing: Creation of the tree with the mean values and RMS at different bins
28// 4. : Paramaterization fitting.
29// In old AliRoot the RMS paramterization was used to characterize cluster resolution
30// Mean value /bias was neglected
31// 5. To be implemented
32// a.) lookup table for the distortion parmaterization - THn
33// b.) Usage in the reconstruction
34//
35//
36// 1. Creation of the histograms: MakeHistos()
37//
38// 6 n dimensional histograms are filled during the calibration:
39// cluster performance bins
40// 0 - variable of interest
41// 1 - pad type - 0- IROC Short 1- OCROC medium 2 - OROC long pads
42// 2 - drift length - drift length -0-1
43// 3 - Qmax - Qmax - 2- 400
44// 4 - cog - COG position - 0-1
45// 5 - tan(phi) - local angle
46// Histograms:
47// THnSparse *fHisDeltaY; // THnSparse - delta Y between the cluster and tracklet interpolation (+-X(5?) rows)
48// THnSparse *fHisDeltaZ; // THnSparse - delta Z
49// THnSparse *fHisRMSY; // THnSparse - rms Y
50// THnSparse *fHisRMSZ; // THnSparse - rms Z
51// THnSparse *fHisQmax; // THnSparse - qmax
52// THnSparse *fHisQtot; // THnSparse - qtot
53//
54//
55//
56
57
58 //
10757ee9 59// The parameter 'clusterParam', a AliTPCClusterParam object //
60// (needed for TPC cluster error and shape parameterization) //
60e31b96 61//
10757ee9 62// Normally you get this object out of the file 'TPCClusterParam.root' //
63// In the parameter 'cuts' the cuts are specified, that decide //
64// weather a track will be accepted for calibration or not. //
65// //
66//
67// The data flow:
68//
69/*
70 Raw Data -> Local Reconstruction -> Tracking -> Calibration -> RefData (component itself)
71 Offline/HLT Offline/HLT OCDB entries (AliTPCClusterParam)
72*/
73
60e31b96 74/*
75 Example usage - creation of the residual trees:
76 TFile f("CalibObjects.root");
77 AliTPCcalibTracks *calibTracks = ( AliTPCcalibTracks *)TPCCalib->FindObject("calibTracks");
78 TH2 * his2 = calibTracks->fHisDeltaY->Projection(0,4);
79 his2->SetName("his2");
80 his2->FitSlicesY();
81
82
83 TTreeSRedirector *pcstream = new TTreeSRedirector("clusterLookup.root");
84 AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaY, pcstream, 0);
85 AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaZ, pcstream, 1);
86 delete pcstream;
87 TFile fl("clusterLookup.root");
88 TCut cutNI="cogA==0&&angleA==0&&qmaxA==0&&zA==0&&ipadA==0"
89
90*/
9f0beaf7 91
10757ee9 92 //
93///////////////////////////////////////////////////////////////////////////////
94
95//
96// ROOT includes
97//
98#include <iostream>
99#include <fstream>
100using namespace std;
101
102#include <TROOT.h>
103#include <TChain.h>
104#include <TFile.h>
105#include <TH3F.h>
5b00528f 106#include <TProfile.h>
107
10757ee9 108//
109//#include <TPDGCode.h>
110#include <TStyle.h>
111#include "TLinearFitter.h"
112//#include "TMatrixD.h"
113#include "TTreeStream.h"
114#include "TF1.h"
115#include <TCanvas.h>
116#include <TGraph2DErrors.h>
117#include "TPostScript.h"
118#include "TCint.h"
119
120#include <TH2D.h>
121#include <TF2.h>
122#include <TSystem.h>
123#include <TCollection.h>
124#include <iostream>
125#include <TLinearFitter.h>
2e5bcb67 126#include <TString.h>
10757ee9 127
128//
129// AliROOT includes
130//
131#include "AliMagF.h"
132#include "AliTracker.h"
133#include "AliESD.h"
134#include "AliESDtrack.h"
135#include "AliESDfriend.h"
136#include "AliESDfriendTrack.h"
137#include "AliTPCseed.h"
138#include "AliTPCclusterMI.h"
139#include "AliTPCROC.h"
140
141#include "AliTPCParamSR.h"
142#include "AliTrackPointArray.h"
143#include "AliTPCcalibTracks.h"
144#include "AliTPCClusterParam.h"
145#include "AliTPCcalibTracksCuts.h"
10757ee9 146#include "AliTPCCalPad.h"
147#include "AliTPCCalROC.h"
148#include "TText.h"
149#include "TPaveText.h"
150#include "TSystem.h"
2e5bcb67 151#include "TStatToolkit.h"
152#include "TCut.h"
af6a50bb 153#include "THnSparse.h"
46e89793 154#include "AliRieman.h"
42d30447 155#include "TRandom.h"
10757ee9 156
10757ee9 157
158
159ClassImp(AliTPCcalibTracks)
160
161
162AliTPCcalibTracks::AliTPCcalibTracks():
c32da879 163 AliTPCcalibBase(),
2c632057 164 fClusterParam(0),
2c632057 165 fROC(0),
af6a50bb 166 fHisDeltaY(0), // THnSparse - delta Y
167 fHisDeltaZ(0), // THnSparse - delta Z
168 fHisRMSY(0), // THnSparse - rms Y
169 fHisRMSZ(0), // THnSparse - rms Z
170 fHisQmax(0), // THnSparse - qmax
171 fHisQtot(0), // THnSparse - qtot
42d30447 172 fPtDownscaleRatio(20), // pt downscaling ratio (use subsample of data)
173 fQDownscaleRatio(20), // Q downscaling ratio (use subsample of dta)
2c632057 174 fArrayQDY(0),
175 fArrayQDZ(0),
176 fArrayQRMSY(0),
177 fArrayQRMSZ(0),
2c632057 178 fResolY(0),
179 fResolZ(0),
180 fRMSY(0),
181 fRMSZ(0),
182 fCuts(0),
2c632057 183 fRejectedTracksHisto(0),
2c632057 184 fClusterCutHisto(0),
185 fCalPadClusterPerPad(0),
9f0beaf7 186 fCalPadClusterPerPadRaw(0)
10757ee9 187{
188 //
189 // AliTPCcalibTracks default constructor
190 //
ae28e92e 191 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;
10757ee9 192}
193
194
b42cf9ed 195
196AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
ae28e92e 197 AliTPCcalibBase(calibTracks),
2c632057 198 fClusterParam(0),
2c632057 199 fROC(0),
af6a50bb 200 fHisDeltaY(0), // THnSparse - delta Y
201 fHisDeltaZ(0), // THnSparse - delta Z
202 fHisRMSY(0), // THnSparse - rms Y
203 fHisRMSZ(0), // THnSparse - rms Z
204 fHisQmax(0), // THnSparse - qmax
205 fHisQtot(0), // THnSparse - qtot
42d30447 206 fPtDownscaleRatio(20), // pt downscaling ratio (use subsample of data)
207 fQDownscaleRatio(20), // Q downscaling ratio (use subsample of dta)
2c632057 208 fArrayQDY(0),
209 fArrayQDZ(0),
210 fArrayQRMSY(0),
211 fArrayQRMSZ(0),
2c632057 212 fResolY(0),
213 fResolZ(0),
214 fRMSY(0),
215 fRMSZ(0),
216 fCuts(0),
2c632057 217 fRejectedTracksHisto(0),
2c632057 218 fClusterCutHisto(0),
219 fCalPadClusterPerPad(0),
9f0beaf7 220 fCalPadClusterPerPadRaw(0)
2c632057 221{
10757ee9 222 //
223 // AliTPCcalibTracks copy constructor
224 //
ae28e92e 225 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
10757ee9 226
227 Bool_t dirStatus = TH1::AddDirectoryStatus();
228 TH1::AddDirectory(kFALSE);
229
230 Int_t length = -1;
10757ee9 231
b42cf9ed 232 (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
10757ee9 233 fArrayQDY= new TObjArray(length);
234 fArrayQDZ= new TObjArray(length);
235 fArrayQRMSY= new TObjArray(length);
236 fArrayQRMSZ= new TObjArray(length);
237 for (Int_t i = 0; i < length; i++) {
b42cf9ed 238 fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
239 fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
240 fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
241 fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
10757ee9 242 }
243
b42cf9ed 244 (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
10757ee9 245 fResolY = new TObjArray(length);
246 fResolZ = new TObjArray(length);
247 fRMSY = new TObjArray(length);
248 fRMSZ = new TObjArray(length);
249 for (Int_t i = 0; i < length; i++) {
b42cf9ed 250 fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
251 fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
252 fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
253 fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
10757ee9 254 }
255
10757ee9 256
b42cf9ed 257 fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
258 fRejectedTracksHisto = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
b42cf9ed 259 fCalPadClusterPerPad = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
260 fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
261
262 fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(),
263 calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
b42cf9ed 264 SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
10757ee9 265 TH1::AddDirectory(dirStatus); // set status back to original status
266// cout << "+++++ end of copy constructor +++++" << endl; // TO BE REMOVED
267}
268
269
b42cf9ed 270AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
271 //
272 // assgnment operator
273 //
274 if (this != &calibTracks) {
275 new (this) AliTPCcalibTracks(calibTracks);
276 }
277 return *this;
278
279}
280
281
10757ee9 282AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts, Int_t logLevel) :
c32da879 283 AliTPCcalibBase(),
2c632057 284 fClusterParam(0),
2c632057 285 fROC(0),
af6a50bb 286 fHisDeltaY(0), // THnSparse - delta Y
287 fHisDeltaZ(0), // THnSparse - delta Z
288 fHisRMSY(0), // THnSparse - rms Y
289 fHisRMSZ(0), // THnSparse - rms Z
290 fHisQmax(0), // THnSparse - qmax
291 fHisQtot(0), // THnSparse - qtot
42d30447 292 fPtDownscaleRatio(20), // pt downscaling ratio (use subsample of data)
293 fQDownscaleRatio(20), // Q downscaling ratio (use subsample of dta)
2c632057 294 fArrayQDY(0),
295 fArrayQDZ(0),
296 fArrayQRMSY(0),
297 fArrayQRMSZ(0),
2c632057 298 fResolY(0),
299 fResolZ(0),
300 fRMSY(0),
301 fRMSZ(0),
302 fCuts(0),
2c632057 303 fRejectedTracksHisto(0),
2c632057 304 fClusterCutHisto(0),
305 fCalPadClusterPerPad(0),
9f0beaf7 306 fCalPadClusterPerPadRaw(0)
2c632057 307 {
10757ee9 308 //
309 // AliTPCcalibTracks constructor
310 // specify 'name' and 'title' of your object
311 // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
312 // In the parameter 'cuts' the cuts are specified, that decide
313 // weather a track will be accepted for calibration or not.
ae28e92e 314 //
315 // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
10757ee9 316 //
317 // All histograms are instatiated in this constructor.
318 //
c32da879 319 this->SetName(name);
320 this->SetTitle(title);
321
ae28e92e 322 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
10757ee9 323 G__SetCatchException(0);
324
325 fClusterParam = clusterParam;
326 if (fClusterParam){
327 fClusterParam->SetInstance(fClusterParam);
328 }
329 else {
330 Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
331 }
332 fCuts = cuts;
ae28e92e 333 SetDebugLevel(logLevel);
af6a50bb 334 MakeHistos();
10757ee9 335
336 TH1::AddDirectory(kFALSE);
337
8b9b3187 338 fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
10757ee9 339 fCalPadClusterPerPad = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
340 fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
341 fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
342
10757ee9 343
344 TH1::AddDirectory(kFALSE);
345
10757ee9 346
347 fResolY = new TObjArray(3);
348 fResolZ = new TObjArray(3);
349 fRMSY = new TObjArray(3);
350 fRMSZ = new TObjArray(3);
351 TH3F * his3D;
352 //
353 his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
354 fResolY->AddAt(his3D,0);
355 his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
356 fResolY->AddAt(his3D,1);
357 his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
358 fResolY->AddAt(his3D,2);
359 //
360 his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
361 fResolZ->AddAt(his3D,0);
362 his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
363 fResolZ->AddAt(his3D,1);
364 his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
365 fResolZ->AddAt(his3D,2);
366 //
367 his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
368 fRMSY->AddAt(his3D,0);
369 his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
370 fRMSY->AddAt(his3D,1);
371 his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
372 fRMSY->AddAt(his3D,2);
373 //
374 his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
375 fRMSZ->AddAt(his3D,0);
376 his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
377 fRMSZ->AddAt(his3D,1);
378 his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
379 fRMSZ->AddAt(his3D,2);
380 //
381
382 TH1::AddDirectory(kFALSE);
383
384 fArrayQDY = new TObjArray(300);
385 fArrayQDZ = new TObjArray(300);
386 fArrayQRMSY = new TObjArray(300);
387 fArrayQRMSZ = new TObjArray(300);
388 for (Int_t iq = 0; iq <= 10; iq++){
389 for (Int_t ipad = 0; ipad < 3; ipad++){
390 Int_t bin = GetBin(iq, ipad);
391 Float_t qmean = GetQ(bin);
d23650ed 392 char hname[200];
f69a0aa7 393 snprintf(hname,100,"ResolY Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 394 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
10757ee9 395 fArrayQDY->AddAt(his3D, bin);
f69a0aa7 396 snprintf(hname,100,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 397 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
10757ee9 398 fArrayQDZ->AddAt(his3D, bin);
f69a0aa7 399 snprintf(hname,100,"RMSY Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 400 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
10757ee9 401 fArrayQRMSY->AddAt(his3D, bin);
f69a0aa7 402 snprintf(hname,100,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 403 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
10757ee9 404 fArrayQRMSZ->AddAt(his3D, bin);
405 }
406 }
407
10757ee9 408
10757ee9 409
ae28e92e 410 if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl;
10757ee9 411 cout << "end of main constructor" << endl; // TO BE REMOVED
412}
413
414
415AliTPCcalibTracks::~AliTPCcalibTracks() {
416 //
417 // AliTPCcalibTracks destructor
418 //
419
ae28e92e 420 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
10757ee9 421 Int_t length = 0;
10757ee9 422
10757ee9 423
f69a0aa7 424 if (fResolY) {
425 length = fResolY->GetEntriesFast();
426 for (Int_t i = 0; i < length; i++){
427 delete fResolY->At(i);
428 delete fResolZ->At(i);
429 delete fRMSY->At(i);
430 delete fRMSZ->At(i);
431 }
432 delete fResolY;
433 delete fResolZ;
434 delete fRMSY;
435 delete fRMSZ;
10757ee9 436 }
10757ee9 437
f69a0aa7 438 if (fArrayQDY) {
439 length = fArrayQDY->GetEntriesFast();
440 for (Int_t i = 0; i < length; i++){
441 delete fArrayQDY->At(i);
442 delete fArrayQDZ->At(i);
443 delete fArrayQRMSY->At(i);
444 delete fArrayQRMSZ->At(i);
445 }
10757ee9 446 }
447
10757ee9 448
9f0beaf7 449
10757ee9 450 delete fArrayQDY;
451 delete fArrayQDZ;
452 delete fArrayQRMSY;
453 delete fArrayQRMSZ;
10757ee9 454
10757ee9 455 delete fRejectedTracksHisto;
456 delete fClusterCutHisto;
10757ee9 457 if (fCalPadClusterPerPad) delete fCalPadClusterPerPad;
458 if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
af6a50bb 459 delete fHisDeltaY; // THnSparse - delta Y
460 delete fHisDeltaZ; // THnSparse - delta Z
461 delete fHisRMSY; // THnSparse - rms Y
462 delete fHisRMSZ; // THnSparse - rms Z
463 delete fHisQmax; // THnSparse - qmax
464 delete fHisQtot; // THnSparse - qtot
465
10757ee9 466}
467
468
10757ee9 469
c32da879 470void AliTPCcalibTracks::Process(AliTPCseed *track){
10757ee9 471 //
472 // To be called in the selector
473 // first AcceptTrack is evaluated, then calls all the following analyse functions:
474 // FillResolutionHistoLocal(track)
af6a50bb 475
10757ee9 476 //
42d30447 477 Double_t scalept= TMath::Min(1./TMath::Abs(track->GetParameter()[4]),2.)/0.5;
478 Bool_t isSelected = (TMath::Exp(scalept)>fPtDownscaleRatio*gRandom->Rndm());
479 if (isSelected) return;
480
ae28e92e 481 if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
10757ee9 482 Int_t accpetStatus = AcceptTrack(track);
483 if (accpetStatus == 0) {
484 FillResolutionHistoLocal(track);
10757ee9 485 }
486 else fRejectedTracksHisto->Fill(accpetStatus);
487}
488
489
490
491Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
492 //
493 // calculate bins for given q and pad type
494 // used in TObjArray
495 //
496 Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );
497 res *= 3;
498 res += pad;
499 return res;
500}
501
502
503Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
504 //
505 // calculate bins for given iq and pad type
506 // used in TObjArray
507 //
508 return iq * 3 + pad;;
509}
510
511
512Float_t AliTPCcalibTracks::GetQ(Int_t bin){
513 //
514 // returns to bin belonging charge
515 // (bin / 3 + 3)^2
516 //
517 Int_t bin0 = bin / 3;
518 bin0 += 3;
519 return bin0 * bin0;
520}
521
522
523Float_t AliTPCcalibTracks::GetPad(Int_t bin){
524 //
525 // returns to bin belonging pad
526 // bin % 3
527 //
528 return bin % 3;
529}
530
531
532
533Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
534 //
535 // Function, that decides wheather a given track is accepted for
536 // the analysis or not.
537 // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
538 // Returns 0 if a track is accepted or an integer different from 0
539 // to indicate the failed cut
540 //
541 const Int_t kMinClusters = fCuts->GetMinClusters();
542 const Float_t kMinRatio = fCuts->GetMinRatio();
543 const Float_t kMax1pt = fCuts->GetMax1pt();
544 const Float_t kEdgeYXCutNoise = fCuts->GetEdgeYXCutNoise();
545 const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
546
547 //
548 // edge induced noise tracks - NEXT RELEASE will be removed during tracking
549 if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
550 if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
551 if (track->GetNumberOfClusters() < kMinClusters) return 2;
552 Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
553 if (ratio < kMinRatio) return 3;
554// Float_t mpt = track->Get1Pt(); // Get1Pt() doesn't exist any more
555 Float_t mpt = track->GetSigned1Pt();
556 if (TMath::Abs(mpt) > kMax1pt) return 4;
557 //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
558 //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
559 //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
560
430a3403 561 if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");
10757ee9 562 return 0;
563}
564
565
566void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
567 //
568 // fill resolution histograms - localy - tracklet in the neighborhood
569 // write debug information to 'TPCSelectorDebug.root'
570 //
571 // _ the main function, called during track analysis _
572 //
573 // loop over all padrows along the track
574 // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
575 //
576 // loop again over all padrows along the track
577 // fit tracklet (clusters in given padrow +- kDelta padrows)
578 // with polynom of 2nd order and two polynoms of 1st order
579 // take both polynoms of 1st order, calculate difference of their parameters
580 // add covariance matrixes and calculate chi2 of this difference
581 // if this chi2 is bigger than a given threshold, assume that the current cluster is
582 // a kink an goto next padrow
583 // if not:
46e89793 584 // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
10757ee9 585 //
586 // write debug information to 'TPCSelectorDebug.root'
587 // only for every kDeltaWriteDebugStream'th padrow to reduce data volume
588 // and to avoid redundant data
589 //
590
7d14c1c1 591 /* {//SG
592 static long n=0;
593 if( n==0 ){
594 if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){
595 cout<<"Map found "<<endl;
596 }else{
597 cout<<"Can't find the map "<<endl;
598 }
599 }
600 if( n%100 ==0 )cout<<n<<" Tracks processed"<<endl;
601 n++;
602 }*/
9f0beaf7 603 static TLinearFitter fFitterParY(3,"pol2"); //
604 static TLinearFitter fFitterParZ(3,"pol2"); //
46e89793 605 static TLinearFitter fFitterParYWeight(3,"pol2"); //
606 static TLinearFitter fFitterParZWeight(3,"pol2"); //
607 fFitterParY.StoreData(kTRUE);
608 fFitterParZ.StoreData(kTRUE);
609 fFitterParYWeight.StoreData(kTRUE);
610 fFitterParZWeight.StoreData(kTRUE);
ae28e92e 611 if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
46e89793 612 const Int_t kDelta = 10; // delta rows to fit
613 const Float_t kMinRatio = 0.75; // minimal ratio
614 const Float_t kChi2Cut = 10.; // cut chi2 - left right
615 const Float_t kSigmaCut = 3.; //sigma cut
616 const Float_t kdEdxCut=300;
617 const Float_t kPtCut=0.300;
618
619 if (track->GetTPCsignal()>kdEdxCut) return;
620 if (TMath::Abs(AliTracker::GetBz())>0.1 &&TMath::Abs(track->Pt())<kPtCut) return;
621
622 // estimate mean error
623 Int_t nTrackletsAll = 0; // number of tracklets for given track
624 Float_t csigmaY = 0; // mean sigma for tracklet refit in Y direction
625 Float_t csigmaZ = 0; // mean sigma for tracklet refit in Z direction
626 Int_t nClusters = 0; // working variable, number of clusters per tracklet
627 Int_t sectorG = -1; // working variable, sector of tracklet, has to stay constant for one tracklet
628 Double_t refX=0;
629 // ---------------------------------------------------------------------
630 for (Int_t irow = 0; irow < 159; irow++){
631 // loop over all rows along the track
632 // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
633 // calculate mean chi^2 for this track-fit in Y and Z direction
634 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
635 if (!cluster0) continue; // no cluster found
636 Int_t sector = cluster0->GetDetector();
637
638 Int_t ipad = TMath::Nint(cluster0->GetPad());
639 Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
640 fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
641
642 if (sector != sectorG){
643 // track leaves sector before it crossed enough rows to fit / initialization
644 nClusters = 0;
9f0beaf7 645 fFitterParY.ClearPoints();
646 fFitterParZ.ClearPoints();
46e89793 647 sectorG = sector;
648 }
649 else {
650 nClusters++;
651 if (refX<1) refX=cluster0->GetX()+kDelta*0.5;
652 Double_t x = cluster0->GetX()-refX;
653 fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
654 fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
655 //
656 if ( nClusters >= kDelta + 3 ){
657 // if more than 13 (kDelta+3) clusters were added to the fitters
658 // fit the tracklet, increase trackletCounter
659 fFitterParY.Eval();
660 fFitterParZ.Eval();
661 nTrackletsAll++;
662 csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
663 csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
664 nClusters = -1;
665 fFitterParY.ClearPoints();
666 fFitterParZ.ClearPoints();
667 refX=0;
10757ee9 668 }
46e89793 669 }
670 } // for (Int_t irow = 0; irow < 159; irow++)
671 // mean chi^2 for all tracklet fits in Y and in Z direction:
672 csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
673 csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
3d34659a 674 if (csigmaY<=TMath::KUncertainty() || csigmaZ<= TMath::KUncertainty()) return;
46e89793 675 // ---------------------------------------------------------------------
676 //
677 //
10757ee9 678
46e89793 679 for (Int_t irow = kDelta; irow < 159-kDelta; irow++){
680 // loop again over all rows along the track
681 // do analysis
682 //
683 Int_t nclFound = 0; // number of clusters in the neighborhood
684 Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster
685 Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster
686 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
687 if (!cluster0) continue;
688 Int_t sector = cluster0->GetDetector();
689 Float_t xref = cluster0->GetX();
690
691 // Make Fit
692 fFitterParY.ClearPoints();
693 fFitterParZ.ClearPoints();
694 fFitterParYWeight.ClearPoints();
695 fFitterParZWeight.ClearPoints();
696 // fit tracklet (clusters in given padrow +- kDelta padrows)
697 // with polynom of 2nd order and two polynoms of 1st order
698 // take both polynoms of 1st order, calculate difference of their parameters
699 // add covariance matrixes and calculate chi2 of this difference
700 // if this chi2 is bigger than a given threshold, assume that the current cluster is
701 // a kink an goto next padrow
7d14c1c1 702
703 // first fit the track angle for wave correction
704
705 AliRieman riemanFitAngle( 2*kDelta ); //SG
706
707 if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){
708 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
709 // loop over irow +- kDelta rows (neighboured rows)
710 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
711 AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
712 if (!currentCluster) continue;
713 if (currentCluster->GetType() < 0) continue;
714 if (currentCluster->GetDetector() != sector) continue;
715 riemanFitAngle.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ);
716 }
717 if( riemanFitAngle.GetN()>3 ) riemanFitAngle.Update();
718 }
719
720 // do fit
721
46e89793 722 AliRieman riemanFit(2*kDelta);
723 AliRieman riemanFitW(2*kDelta);
724 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
725 // loop over irow +- kDelta rows (neighboured rows)
726 //
727 //
728 if (idelta == 0) continue; // don't use center cluster
729 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
730 AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
731 if (!currentCluster) continue;
732 if (currentCluster->GetType() < 0) continue;
733 if (currentCluster->GetDetector() != sector) continue;
734 nclFound++;
735 if (idelta < 0){
736 ncl0++;
9f0beaf7 737 }
46e89793 738 if (idelta > 0){
739 ncl1++;
10757ee9 740 }
7d14c1c1 741 //SG!!!
742 double dY = 0;
743 if( fClusterParam ){
744 Int_t padSize = 0; // short pads
745 if (currentCluster->GetDetector() >= 36) {
746 padSize = 1; // medium pads
747 if (currentCluster->GetRow() > 63) padSize = 2; // long pads
748 }
749 dY = fClusterParam->GetWaveCorrection( padSize,
750 currentCluster->GetZ(),
751 currentCluster->GetMax(),
752 currentCluster->GetPad(),
753 riemanFitAngle.GetDYat(currentCluster->GetX())
754 );
755 }
3d34659a 756 riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), csigmaY,csigmaZ);
757 riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), TMath::Abs(csigmaY*TMath::Sqrt(1+TMath::Abs(idelta))),TMath::Abs(csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta))));
46e89793 758 } // loop over neighbourhood for fitter filling
759 if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
760 if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
761 if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow
762 riemanFit.Update();
763 riemanFitW.Update();
764 Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5));
765 Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5));
766 Double_t paramYR[3], paramZR[3];
767 Double_t paramYRW[3], paramZRW[3];
768 //
769 paramYR[0] = riemanFit.GetYat(xref);
770 paramZR[0] = riemanFit.GetZat(xref);
771 paramYRW[0] = riemanFitW.GetYat(xref);
772 paramZRW[0] = riemanFitW.GetZat(xref);
773 //
774 paramYR[1] = riemanFit.GetDYat(xref);
775 paramZR[1] = riemanFit.GetDZat(xref);
776 paramYRW[1] = riemanFitW.GetDYat(xref);
777 paramZRW[1] = riemanFitW.GetDZat(xref);
778 //
779 Int_t reject=0;
780 if (chi2R>kChi2Cut) reject+=1;
781 if (chi2RW>kChi2Cut) reject+=2;
782 if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4;
783 if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8;
784 if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16;
785 if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32;
786 //
787 TTreeSRedirector *cstream = GetDebugStreamer();
788 // get fit parameters from pol2 fit:
789 Double_t tracky = paramYR[0];
790 Double_t trackz = paramZR[0];
791 Float_t deltay = cluster0->GetY()-tracky;
792 Float_t deltaz = cluster0->GetZ()-trackz;
793 Float_t angley = paramYR[1];
794 Float_t anglez = paramZR[1];
795
796 Int_t padSize = 0; // short pads
797 if (cluster0->GetDetector() >= 36) {
10757ee9 798 padSize = 1; // medium pads
799 if (cluster0->GetRow() > 63) padSize = 2; // long pads
46e89793 800 }
801 Int_t ipad = TMath::Nint(cluster0->GetPad());
802 if (cstream){
803 Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
804 (*cstream)<<"Resol0"<<
805 "run="<<fRun<< // run number
806 "event="<<fEvent<< // event number
807 "time="<<fTime<< // time stamp of event
808 "trigger="<<fTrigger<< // trigger
809 "mag="<<fMagF<< // magnetic field
810 "padSize="<<padSize<<
811 //
812 "reject="<<reject<<
813 "cl.="<<cluster0<< // cluster info
814 "xref="<<xref<< // cluster reference X
815 //rieman fit
816 "yr="<<paramYR[0]<< // track position y - rieman fit
817 "zr="<<paramZR[0]<< // track position z - rieman fit
818 "yrW="<<paramYRW[0]<< // track position y - rieman fit - weight
819 "zrW="<<paramZRW[0]<< // track position z - rieman fit - weight
820 "dyr="<<paramYR[1]<< // track position y - rieman fit
821 "dzr="<<paramZR[1]<< // track position z - rieman fit
822 "dyrW="<<paramYRW[1]<< // track position y - rieman fit - weight
823 "dzrW="<<paramZRW[1]<< // track position z - rieman fit - weight
824 //
825 "angley="<<angley<< // angle par fit
826 "anglez="<<anglez<< // angle par fit
827 "zdr="<<zdrift<< //
828 "dy="<<deltay<<
829 "dz="<<deltaz<<
830 "sy="<<csigmaY<<
831 "sz="<<csigmaZ<<
832 "chi2R="<<chi2R<<
833 "chi2RW="<<chi2RW<<
834 "\n";
835 }
836 if (reject) continue;
837
838
839 // =========================================
840 // wirte collected information to histograms
841 // =========================================
842
843 Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
844 fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
845
846
847 TH3F * his3 = 0;
848 his3 = (TH3F*)fRMSY->At(padSize);
7e22eccb 849 if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
46e89793 850 his3 = (TH3F*)fRMSZ->At(padSize);
7e22eccb 851 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
46e89793 852
853 his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
7e22eccb 854 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
46e89793 855 his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
7e22eccb 856 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
46e89793 857
858
859 his3 = (TH3F*)fResolY->At(padSize);
860 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
861 his3 = (TH3F*)fResolZ->At(padSize);
862 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
863 his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
864 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
865 his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
866 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
867 //=============================================================================================
868 //
869 // Fill THN histograms
870 //
42d30447 871 Double_t scaleQ= TMath::Min(Double_t(cluster0->GetMax()),200.)/50.;
872 Bool_t isSelected = (TMath::Exp(scaleQ)>fQDownscaleRatio*gRandom->Rndm());
873
46e89793 874 Double_t xvar[9];
60e31b96 875 xvar[1]=padSize; // pad type
876 xvar[2]=cluster0->GetZ(); //
46e89793 877 xvar[3]=cluster0->GetMax();
878
879 xvar[0]=deltay;
7d14c1c1 880 double cog = cluster0->GetPad()-Int_t(cluster0->GetPad());// distance to the center of the pad
881 xvar[4] = cog;
882 xvar[5]=angley;
883
884 if( TMath::Abs(cog-0.5)<1.e-8 ) xvar[4]=-1; // fill one-pad clusters in the underflow bin
885 fHisDeltaY->Fill(xvar);
886
887 xvar[4]=cog;
46e89793 888 xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
889 fHisRMSY->Fill(xvar);
890
891 xvar[0]=deltaz;
60e31b96 892 xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin()); // distance to the center of the time bin
46e89793 893 xvar[5]=anglez;
894 fHisDeltaZ->Fill(xvar);
895 xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
896 fHisRMSZ->Fill(xvar);
897
898 } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
899} // FillResolutionHistoLocal(...)
900
10757ee9 901
902
903
904
905
10757ee9 906
907
908void AliTPCcalibTracks::SetStyle() const {
909 //
910 // set style, can be called by all draw functions
911 //
912 gROOT->SetStyle("Plain");
913 gStyle->SetFillColor(10);
914 gStyle->SetPadColor(10);
915 gStyle->SetCanvasColor(10);
916 gStyle->SetStatColor(10);
917 gStyle->SetPalette(1,0);
918 gStyle->SetNumberContours(60);
919}
920
921
10757ee9 922
f78da5ae 923void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
10757ee9 924 //
925 // all functions are called, that produce pictures
926 // the histograms are written to the directory 'pathName'
927 // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
928 // 'stat' is also the number of minEntries for MakeResPlotsQTree
929 //
930
ae28e92e 931 if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
46e89793 932 MakeResPlotsQTree(stat, pathName);
10757ee9 933}
934
935
10757ee9 936
937
f78da5ae 938void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
10757ee9 939 //
940 // Make tree with resolution parameters
941 // the result is written to 'resol.root' in directory 'pathname'
942 // file information are available in fileInfo
943 // available variables in the tree 'Resol':
944 // Entries: number of entries for this resolution point
945 // nbins: number of bins that were accumulated
946 // Dim: direction, Dim==0: y-direction, Dim==1: z-direction
947 // Pad: padSize; short, medium and long
948 // Length: pad length, 0.75, 1, 1.5
949 // QMean: mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
950 // Zc: center of middle bin in drift direction
951 // Zm: mean dirftlength for accumulated Delta-Histograms
952 // Zs: width of driftlength bin
953 // AngleC: center of middle bin in Angle-Direction
954 // AngleM: mean angle for accumulated Delta-Histograms
955 // AngleS: width of Angle-bin
956 // Resol: sigma for gaus fit through Delta-Histograms
957 // Sigma: error of sigma for gaus fit through Delta Histograms
958 // MeanR: mean of the Delta-Histogram
959 // SigmaR: rms of the Delta-Histogram
960 // RMSm: mean of the gaus fit through RMS-Histogram
961 // RMS: sigma of the gaus fit through RMS-Histogram
962 // RMSe0: error of mean of gaus fit in RMS-Histogram
963 // RMSe1: error of sigma of gaus fit in RMS-Histogram
964 //
965
ae28e92e 966 if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
967 if (GetDebugLevel() > -1) cout << " relax, the calculation will take a while..." << endl;
10757ee9 968
969 gSystem->MakeDirectory(pathName);
970 gSystem->ChangeDirectory(pathName);
971 TString kFileName = "resol.root";
972 TTreeSRedirector fTreeResol(kFileName.Data());
973
974 TH3F *resArray[2][3][11];
975 TH3F *rmsArray[2][3][11];
976
977 // load histograms from fArrayQDY and fArrayQDZ
978 // into resArray and rmsArray
979 // that is all we need here
980 for (Int_t idim = 0; idim < 2; idim++){
981 for (Int_t ipad = 0; ipad < 3; ipad++){
982 for (Int_t iq = 0; iq <= 10; iq++){
983 rmsArray[idim][ipad][iq]=0;
984 resArray[idim][ipad][iq]=0;
985 Int_t bin = GetBin(iq,ipad);
986 TH3F *hresl = 0;
987 if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
988 if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
989 if (!hresl) continue;
990 resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
991 resArray[idim][ipad][iq]->SetDirectory(0);
992 TH3F * hreslRMS = 0;
993 if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
994 if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
995 if (!hreslRMS) continue;
996 rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
997 rmsArray[idim][ipad][iq]->SetDirectory(0);
998 }
999 }
1000 }
ae28e92e 1001 if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
10757ee9 1002
1003 //--------------------------------------------------------------------------------------------
1004
1005 char name[200];
1006 Double_t qMean;
1007 Double_t zMean, angleMean, zCenter, angleCenter;
1008 Double_t zSigma, angleSigma;
1009 TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1010 TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1011 TF1 *fitFunction = new TF1("fitFunction", "gaus");
1012 Float_t entriesQ = 0;
1013 Int_t loopCounter = 1;
1014
1015 for (Int_t idim = 0; idim < 2; idim++){
1016 // Loop y-z corrdinate
1017 for (Int_t ipad = 0; ipad < 3; ipad++){
1018 // loop pad type
1019 for (Int_t iq = -1; iq < 10; iq++){
1020 // LOOP Q
ae28e92e 1021 if (GetDebugLevel() > -1)
10757ee9 1022 cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
1023 << (Int_t)((loopCounter)/66.*100) << "% done), "
1024 << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << " \r" << std::flush;
1025 loopCounter++;
1026 TH3F *hres = 0;
1027 TH3F *hrms = 0;
1028 qMean = 0;
1029 entriesQ = 0;
1030
1031 // calculate qMean
1032 if (iq == -1){
1033 // integrated spectra
1034 for (Int_t iql = 0; iql < 10; iql++){
1035 Int_t bin = GetBin(iql,ipad);
1036 TH3F *hresl = resArray[idim][ipad][iql];
1037 TH3F *hrmsl = rmsArray[idim][ipad][iql];
1038 if (!hresl) continue;
1039 if (!hrmsl) continue;
1040 entriesQ += hresl->GetEntries();
1041 qMean += hresl->GetEntries() * GetQ(bin);
1042 if (!hres) {
1043 hres = (TH3F*)hresl->Clone();
1044 hrms = (TH3F*)hrmsl->Clone();
1045 }
1046 else{
1047 hres->Add(hresl);
1048 hrms->Add(hrmsl);
1049 }
1050 }
1051 qMean /= entriesQ;
1052 qMean *= -1.; // integral mean charge
1053 }
1054 else {
1055 // loop over neighboured Q-bins
1056 // accumulate entries from neighboured Q-bins
1057 for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
1058 if (iql < 0) continue;
1059 Int_t bin = GetBin(iql,ipad);
1060 TH3F * hresl = resArray[idim][ipad][iql];
1061 TH3F * hrmsl = rmsArray[idim][ipad][iql];
1062 if (!hresl) continue;
1063 if (!hrmsl) continue;
1064 entriesQ += hresl->GetEntries();
1065 qMean += hresl->GetEntries() * GetQ(bin);
1066 if (!hres) {
1067 hres = (TH3F*) hresl->Clone();
1068 hrms = (TH3F*) hrmsl->Clone();
1069 }
1070 else{
1071 hres->Add(hresl);
1072 hrms->Add(hrmsl);
1073 }
1074 }
1075 qMean/=entriesQ;
1076 }
f69a0aa7 1077 if (!hres) continue;
1078 if (!hrms) continue;
1079
10757ee9 1080 TAxis *xAxisDriftLength = hres->GetXaxis(); // driftlength / z - axis
1081 TAxis *yAxisAngle = hres->GetYaxis(); // angle axis
1082 TAxis *zAxisDelta = hres->GetZaxis(); // delta axis
1083 TAxis *zAxisRms = hrms->GetZaxis(); // rms axis
1084
1085 // loop over all angle bins
1086 for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1087 angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1088 // loop over all driftlength bins
1089 for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1090 zCenter = xAxisDriftLength->GetBinCenter(ibinxDL);
1091 zSigma = xAxisDriftLength->GetBinWidth(ibinxDL);
1092 angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
1093 zMean = zCenter; // changens, when more statistic is accumulated
1094 angleMean = angleCenter; // changens, when more statistic is accumulated
1095
1096 // create 2 1D-Histograms, projectionRes and projectionRms
1097 // these histograms are delta histograms for given direction, padSize, chargeBin,
1098 // angleBin and driftLengthBin
1099 // later on they will be fitted with a gausian, its sigma is the resoltuion...
4aa37f93 1100 snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
10757ee9 1101 // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1102 projectionRes->SetNameTitle(name, name);
4aa37f93 1103 snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
10757ee9 1104 // TH1D * projectionRms = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1105 projectionRms->SetNameTitle(name, name);
1106
1107 projectionRes->Reset();
1108 projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1109 projectionRms->Reset();
1110 projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1111 projectionRes->SetDirectory(0);
1112 projectionRms->SetDirectory(0);
1113
1114 Double_t entries = 0;
1115 Int_t nbins = 0; // counts, how many bins were accumulated
1116 zMean = 0;
1117 angleMean = 0;
1118 entries = 0;
1119
1120 // fill projectionRes and projectionRms for given dim, ipad and iq,
1121 // as well as for given angleBin and driftlengthBin
1122 // if this gives not enough statistic, include neighbourhood
1123 // (angle and driftlength) successifely
1124 for (Int_t dbin = 0; dbin <= 8; dbin++){ // delta-bins around centered angleBin and driftlengthBin
1125 for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) { // delta-bins in angle direction
1126 for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
1127 if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; // add each bin only one time !
1128 Int_t binx2 = ibinxDL + dbinx2; // position variable in x (driftlength) direction
1129 Int_t biny2 = ibinyAngle + dbiny2; // position variable in y (angle) direction
1130 if (binx2 < 1 || biny2 < 1) continue; // don't go out of the histogram!
1131 if (binx2 >= xAxisDriftLength->GetNbins()) continue; // don't go out of the histogram!
1132 if (biny2 >= yAxisAngle->GetNbins()) continue; // don't go out of the histogram!
1133 nbins++; // count the number of accumulated bins
1134 // Fill resolution histo
1135 for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
1136 // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3); // unused variable
1137 projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
1138 entries += hres->GetBinContent(binx2, biny2, ibin3);
1139 zMean += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
1140 angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
1141 } // ibin3 loop
1142 // fill RMS histo
1143 for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
1144 projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
1145 }
1146 } //dbinx2 loop
1147 if (entries > minEntries) break; // enough statistic accumulated
1148 } // dbiny2 loop
1149 if (entries > minEntries) break; // enough statistic accumulated
1150 } // dbin loop
1151 if ( entries< minEntries) continue; // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
1152 zMean /= entries;
1153 angleMean /= entries;
1154
1155 if (entries > minEntries) {
1156 // when enough statistic is accumulated
1157 // fit Delta histograms with a gausian
1158 // of the gausian is the resolution (resol), its fit error is sigma
1159 // store also mean and RMS of the histogram
1160 Float_t xmin = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1161 Float_t xmax = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1162
1163// projectionRes->Fit("gaus", "q0", "", xmin, xmax);
1164// Float_t resol = projectionRes->GetFunction("gaus")->GetParameter(2);
1165// Float_t sigma = projectionRes->GetFunction("gaus")->GetParError(2);
1166 fitFunction->SetMaximum(xmax);
1167 fitFunction->SetMinimum(xmin);
1168 projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
1169 Float_t resol = fitFunction->GetParameter(2);
1170 Float_t sigma = fitFunction->GetParError(2);
1171
1172 Float_t meanR = projectionRes->GetMean();
1173 Float_t sigmaR = projectionRes->GetRMS();
1174 // fit also RMS histograms with a gausian
1175 // store mean and sigma of the gausian in rmsMean and rmsSigma
1176 // store also the fit errors in errorRMS and errorSigma
1177 xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1178 xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1179
1180// projectionRms->Fit("gaus","q0","",xmin,xmax);
1181// Float_t rmsMean = projectionRms->GetFunction("gaus")->GetParameter(1);
1182// Float_t rmsSigma = projectionRms->GetFunction("gaus")->GetParameter(2);
1183// Float_t errorRMS = projectionRms->GetFunction("gaus")->GetParError(1);
1184// Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
1185 projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
1186 Float_t rmsMean = fitFunction->GetParameter(1);
1187 Float_t rmsSigma = fitFunction->GetParameter(2);
1188 Float_t errorRMS = fitFunction->GetParError(1);
1189 Float_t errorSigma = fitFunction->GetParError(2);
1190
1191 Float_t length = 0.75;
1192 if (ipad == 1) length = 1;
1193 if (ipad == 2) length = 1.5;
1194
1195 fTreeResol<<"Resol"<<
1196 "Entries="<<entries<< // number of entries for this resolution point
1197 "nbins="<<nbins<< // number of bins that were accumulated
1198 "Dim="<<idim<< // direction, Dim==0: y-direction, Dim==1: z-direction
1199 "Pad="<<ipad<< // padSize; short, medium and long
1200 "Length="<<length<< // pad length, 0.75, 1, 1.5
1201 "QMean="<<qMean<< // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1202 "Zc="<<zCenter<< // center of middle bin in drift direction
1203 "Zm="<<zMean<< // mean dirftlength for accumulated Delta-Histograms
1204 "Zs="<<zSigma<< // width of driftlength bin
1205 "AngleC="<<angleCenter<< // center of middle bin in Angle-Direction
1206 "AngleM="<<angleMean<< // mean angle for accumulated Delta-Histograms
1207 "AngleS="<<angleSigma<< // width of Angle-bin
1208 "Resol="<<resol<< // sigma for gaus fit through Delta-Histograms
1209 "Sigma="<<sigma<< // error of sigma for gaus fit through Delta Histograms
1210 "MeanR="<<meanR<< // mean of the Delta-Histogram
1211 "SigmaR="<<sigmaR<< // rms of the Delta-Histogram
1212 "RMSm="<<rmsMean<< // mean of the gaus fit through RMS-Histogram
1213 "RMSs="<<rmsSigma<< // sigma of the gaus fit through RMS-Histogram
1214 "RMSe0="<<errorRMS<< // error of mean of gaus fit in RMS-Histogram
1215 "RMSe1="<<errorSigma<< // error of sigma of gaus fit in RMS-Histogram
1216 "\n";
ae28e92e 1217 if (GetDebugLevel() > 5) {
10757ee9 1218 projectionRes->SetDirectory(fTreeResol.GetFile());
1219 projectionRes->Write(projectionRes->GetName());
1220 projectionRes->SetDirectory(0);
1221 projectionRms->SetDirectory(fTreeResol.GetFile());
1222 projectionRms->Write(projectionRms->GetName());
1223 projectionRes->SetDirectory(0);
1224 }
1225 } // if (projectionRes->GetSum() > minEntries)
1226 } // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
1227 } // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
1228
1229 } // iq-loop
1230 } // ipad-loop
1231 } // idim-loop
1232 delete projectionRes;
1233 delete projectionRms;
1234
1235// TFile resolFile(fTreeResol.GetFile());
1236 TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
1237 fileInfo.Write("fileInfo");
1238// resolFile.Close();
1239// fTreeResol.GetFile()->Close();
ae28e92e 1240 if (GetDebugLevel() > -1) cout << endl;
1241 if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
10757ee9 1242 gSystem->ChangeDirectory("..");
1243}
1244
1245
10757ee9 1246
10757ee9 1247
1248
1249Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
1250 //
1251 // function to merge several AliTPCcalibTracks objects after PROOF calculation
1252 // The object's histograms are merged via their merge functions
1253 // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
1254 //
1255
ae28e92e 1256 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl;
10757ee9 1257 if (!collectionList) return 0;
1258 if (collectionList->IsEmpty()) return -1;
1259
ae28e92e 1260 if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
1261 if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
10757ee9 1262 collectionList->Print();
1263
1264 // create a list for each data member
1265 TList* deltaYList = new TList;
1266 TList* deltaZList = new TList;
1267 TList* arrayAmpRowList = new TList;
1268 TList* rejectedTracksList = new TList;
10757ee9 1269 TList* clusterCutHistoList = new TList;
1270 TList* arrayAmpList = new TList;
1271 TList* arrayQDYList = new TList;
1272 TList* arrayQDZList = new TList;
1273 TList* arrayQRMSYList = new TList;
1274 TList* arrayQRMSZList = new TList;
10757ee9 1275 TList* resolYList = new TList;
1276 TList* resolZList = new TList;
1277 TList* rMSYList = new TList;
1278 TList* rMSZList = new TList;
1279
1280// TList* nRowsList = new TList;
1281// TList* nSectList = new TList;
1282// TList* fileNoList = new TList;
1283
1284 TIterator *listIterator = collectionList->MakeIterator();
1285 AliTPCcalibTracks *calibTracks = 0;
ae28e92e 1286 if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
10757ee9 1287 Int_t counter = 0;
a4c5fca9 1288 while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
10757ee9 1289 // loop over all entries in the collectionList and get dataMembers into lists
10757ee9 1290
10757ee9 1291 arrayQDYList->Add(calibTracks->GetfArrayQDY());
1292 arrayQDZList->Add(calibTracks->GetfArrayQDZ());
1293 arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
1294 arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
1295 resolYList->Add(calibTracks->GetfResolY());
1296 resolZList->Add(calibTracks->GetfResolZ());
1297 rMSYList->Add(calibTracks->GetfRMSY());
1298 rMSZList->Add(calibTracks->GetfRMSZ());
10757ee9 1299 rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
1300 clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
a4c5fca9 1301 //
1302 if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
1303 fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
1304 // fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
10757ee9 1305 counter++;
ae28e92e 1306 if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
af6a50bb 1307 AddHistos(calibTracks);
10757ee9 1308 }
1309
10757ee9 1310
1311 // merge data members
ae28e92e 1312 if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
10757ee9 1313 fClusterCutHisto->Merge(clusterCutHistoList);
1314 fRejectedTracksHisto->Merge(rejectedTracksList);
10757ee9 1315
1316 TObjArray* objarray = 0;
1317 TH1* hist = 0;
1318 TList* histList = 0;
1319 TIterator *objListIterator = 0;
1320
46e89793 1321
ae28e92e 1322 if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
10757ee9 1323 // merge fArrayQDY
1324 for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
1325 objListIterator = arrayQDYList->MakeIterator();
1326 histList = new TList;
1327 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1328 // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1329 hist = (TH3F*)(objarray->At(i));
1330 histList->Add(hist);
1331 }
1332 ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
1333 delete histList;
1334 delete objListIterator;
1335 }
1336
ae28e92e 1337 if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
10757ee9 1338 // merge fArrayQDZ
1339 for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1340 objListIterator = arrayQDZList->MakeIterator();
1341 histList = new TList;
1342 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1343 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1344 hist = (TH3F*)(objarray->At(i));
1345 histList->Add(hist);
1346 }
1347 ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
1348 delete histList;
1349 delete objListIterator;
1350 }
1351
ae28e92e 1352 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
10757ee9 1353 // merge fArrayQRMSY
1354 for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
1355 objListIterator = arrayQRMSYList->MakeIterator();
1356 histList = new TList;
1357 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1358 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
03533457 1359 // if (!objarray) continue; // removed for coverity -> JMT
10757ee9 1360 hist = (TH3F*)(objarray->At(i));
1361 histList->Add(hist);
1362 }
1363 ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
1364 delete histList;
1365 delete objListIterator;
1366 }
1367
ae28e92e 1368 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
10757ee9 1369 // merge fArrayQRMSZ
1370 for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1371 objListIterator = arrayQRMSZList->MakeIterator();
1372 histList = new TList;
1373 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1374 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1375 hist = (TH3F*)(objarray->At(i));
1376 histList->Add(hist);
1377 }
1378 ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
1379 delete histList;
1380 delete objListIterator;
1381 }
1382
10757ee9 1383
10757ee9 1384
1385
1386
1387
ae28e92e 1388 if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
10757ee9 1389 // merge fResolY
1390 for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
1391 objListIterator = resolYList->MakeIterator();
1392 histList = new TList;
1393 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1394 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1395 hist = (TH3F*)(objarray->At(i));
1396 histList->Add(hist);
1397 }
1398 ((TH3F*)(fResolY->At(i)))->Merge(histList);
1399 delete histList;
1400 delete objListIterator;
1401 }
1402
1403 // merge fResolZ
1404 for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1405 objListIterator = resolZList->MakeIterator();
1406 histList = new TList;
1407 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1408 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1409 hist = (TH3F*)(objarray->At(i));
1410 histList->Add(hist);
1411 }
1412 ((TH3F*)(fResolZ->At(i)))->Merge(histList);
1413 delete histList;
1414 delete objListIterator;
1415 }
1416
1417 // merge fRMSY
1418 for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
1419 objListIterator = rMSYList->MakeIterator();
1420 histList = new TList;
1421 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1422 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1423 hist = (TH3F*)(objarray->At(i));
1424 histList->Add(hist);
1425 }
1426 ((TH3F*)(fRMSY->At(i)))->Merge(histList);
1427 delete histList;
1428 delete objListIterator;
1429 }
1430
1431 // merge fRMSZ
1432 for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1433 objListIterator = rMSZList->MakeIterator();
1434 histList = new TList;
1435 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1436 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1437 hist = (TH3F*)(objarray->At(i));
1438 histList->Add(hist);
1439 }
1440 ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
1441 delete histList;
1442 delete objListIterator;
1443 }
1444
1445 delete deltaYList;
1446 delete deltaZList;
1447 delete arrayAmpRowList;
1448 delete arrayAmpList;
1449 delete arrayQDYList;
1450 delete arrayQDZList;
1451 delete arrayQRMSYList;
1452 delete arrayQRMSZList;
1453 delete resolYList;
1454 delete resolZList;
1455 delete rMSYList;
1456 delete rMSZList;
10757ee9 1457 delete listIterator;
1458
ae28e92e 1459 if (GetDebugLevel() > 0) cout << "merging done!" << endl;
10757ee9 1460
1461 return 1;
1462}
1463
1464
10757ee9 1465
1466
af6a50bb 1467void AliTPCcalibTracks::MakeHistos(){
1468 //
1469 ////make THnSparse
1470 //
1471 //THnSparse *fHisDeltaY; // THnSparse - delta Y
1472 //THnSparse *fHisDeltaZ; // THnSparse - delta Z
1473 //THnSparse *fHisRMSY; // THnSparse - rms Y
1474 //THnSparse *fHisRMSZ; // THnSparse - rms Z
1475 //THnSparse *fHisQmax; // THnSparse - qmax
1476 //THnSparse *fHisQtot; // THnSparse - qtot
1477 // cluster performance bins
1478 // 0 - variable of interest
1479 // 1 - pad type - 0- short 1-medium 2-long pads
1480 // 2 - drift length - drift length -0-1
1481 // 3 - Qmax - Qmax - 2- 400
1482 // 4 - cog - COG position - 0-1
1483 // 5 - tan(phi) - local y angle
1484 // 6 - tan(theta) - local z angle
1485 // 7 - sector - sector number
1486 Double_t xminTrack[8], xmaxTrack[8];
1487 Int_t binsTrack[8];
1488 TString axisName[8];
1489
1490 //
46e89793 1491 binsTrack[0]=200;
af6a50bb 1492 axisName[0] ="var";
1493
1494 binsTrack[1] =3;
1495 xminTrack[1] =0; xmaxTrack[1]=3;
1496 axisName[1] ="pad type";
1497 //
46e89793 1498 binsTrack[2] =20;
1499 xminTrack[2] =-250; xmaxTrack[2]=250;
1500 axisName[2] ="z";
af6a50bb 1501 //
1502 binsTrack[3] =10;
1503 xminTrack[3] =1; xmaxTrack[3]=400;
1504 axisName[3] ="Qmax";
1505 //
46e89793 1506 binsTrack[4] =20;
af6a50bb 1507 xminTrack[4] =0; xmaxTrack[4]=1;
1508 axisName[4] ="cog";
1509 //
46e89793 1510 binsTrack[5] =15;
1511 xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
1512 axisName[5] ="tan(angle)";
af6a50bb 1513 //
af6a50bb 1514 //
46e89793 1515 xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1516 fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1517 xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1518 fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1519 xminTrack[0] =0.; xmaxTrack[0]=0.5;
46e89793 1520 fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1521 xminTrack[0] =0.; xmaxTrack[0]=0.5;
46e89793 1522 fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1523 xminTrack[0] =0.; xmaxTrack[0]=100;
46e89793 1524 fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1525
1526 xminTrack[0] =0.; xmaxTrack[0]=250;
46e89793 1527 fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1528
1529
1530 for (Int_t ivar=0;ivar<6;ivar++){
1531 fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1532 fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1533 fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1534 fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1535 fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1536 fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1537 fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
1538 fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1539 fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1540 fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1541 fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1542 fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1543 }
1544
1545
af6a50bb 1546 BinLogX(fHisDeltaY,3);
1547 BinLogX(fHisDeltaZ,3);
1548 BinLogX(fHisRMSY,3);
1549 BinLogX(fHisRMSZ,3);
1550 BinLogX(fHisQmax,3);
1551 BinLogX(fHisQtot,3);
1552
1553}
1554
1555void AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
1556 //
1557 // Add histograms
1558 //
1559 if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
1560 if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
1561 if (calib->fHisRMSY) fHisRMSY->Add(calib->fHisRMSY);
1562 if (calib->fHisRMSZ) fHisRMSZ->Add(calib->fHisRMSZ);
1563}
46e89793 1564
1565
1566
1567void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
1568 //
1569 // Dump summary info
1570 //
1571 // 0.OBJ: TAxis var
1572 // 1.OBJ: TAxis pad type
1573 // 2.OBJ: TAxis z
1574 // 3.OBJ: TAxis Qmax
1575 // 4.OBJ: TAxis cog
1576 // 5.OBJ: TAxis tan(angle)
1577 //
1578 if (ptype>3) return;
1579 Int_t idim[6]={0,1,2,3,4,5};
1580 TString hname[4]={"dy","dz","rmsy","rmsz"};
1581 //
1582 Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
1583 Int_t first5=hisInput->GetAxis(5)->GetFirst();
1584 Int_t last5 =hisInput->GetAxis(5)->GetLast();
1585 //
1586 for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){ // axis 5 - local angle
1587 Bool_t bin5All=kTRUE;
1588 if (ibin5>=first5){
1589 hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
1590 if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
1591 bin5All=kFALSE;
1592 }
1593 Double_t x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
1594 THnSparse * his5= hisInput->Projection(5,idim); //projected histogram according selection 5
1595 Int_t nbins4=his5->GetAxis(4)->GetNbins();
1596 Int_t first4=his5->GetAxis(4)->GetFirst();
1597 Int_t last4 =his5->GetAxis(4)->GetLast();
1598
1599 for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){ // axis 4 - cog
1600 Bool_t bin4All=kTRUE;
1601 if (ibin4>=first4){
1602 his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
1603 if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
1604 bin4All=kFALSE;
1605 }
1606 Double_t x4= his5->GetAxis(4)->GetBinCenter(ibin4);
1607 THnSparse * his4= his5->Projection(4,idim); //projected histogram according selection 4
1608 //
1609 Int_t nbins3=his4->GetAxis(3)->GetNbins();
1610 Int_t first3=his4->GetAxis(3)->GetFirst();
1611 Int_t last3 =his4->GetAxis(3)->GetLast();
1612 //
1613 for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - Qmax
1614 Bool_t bin3All=kTRUE;
1615 if (ibin3>=first3){
1616 his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
1617 bin3All=kFALSE;
1618 }
1619 Double_t x3= his4->GetAxis(3)->GetBinCenter(ibin3);
1620 THnSparse * his3= his4->Projection(3,idim); //projected histogram according selection 3
1621 //
1622 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1623 Int_t first2 = his3->GetAxis(2)->GetFirst();
1624 Int_t last2 = his3->GetAxis(2)->GetLast();
1625 //
1626 for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){ // axis 2 - z
1627 Bool_t bin2All=kTRUE;
1628 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1629 if (ibin2>=first2){
1630 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
1631 if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
1632 bin2All=kFALSE;
1633 }
1634 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
1635 //
1636 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
1637 Int_t first1 = his2->GetAxis(1)->GetFirst();
1638 Int_t last1 = his2->GetAxis(1)->GetLast();
1639 for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){ //axis 1 - pad type
1640 Bool_t bin1All=kTRUE;
1641 if (ibin1>=first1){
1642 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1643 bin1All=kFALSE;
1644 }
1645 Double_t x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
1646 TH1 * hisDelta = his2->Projection(0);
1647 Double_t entries = hisDelta->GetEntries();
1648 Double_t mean=0, rms=0;
1649 if (entries>10){
1650 mean = hisDelta->GetMean();
1651 rms = hisDelta->GetRMS();
1652 hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
1653 mean = hisDelta->GetMean();
1654 rms = hisDelta->GetRMS();
1655 }
1656 //
1657 (*pcstream)<<hname[ptype].Data()<<
1658 // flag - intgrated values for given bin
60e31b96 1659 "angleA="<<bin5All<<
46e89793 1660 "cogA="<<bin4All<<
1661 "qmaxA="<<bin3All<<
1662 "zA="<<bin2All<<
1663 "ipadA="<<bin1All<<
1664 // center of bin value
60e31b96 1665 "angle="<<x5<< // local angle
1666 "cog="<<x4<< // distance cluster to center
1667 "qmax="<<x3<< // qmax
1668 "z="<<x2<< // z of the cluster
1669 "ipad="<<x1<< // type of the region
46e89793 1670 // mean values
1671 //
1672 "entries="<<entries<<
1673 "mean="<<mean<<
1674 "rms="<<rms<<
1675 "ptype="<<ptype<< //
1676 "\n";
1677 delete hisDelta;
1678 printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);
1679 }//loop z
1680 delete his2;
1681 } // loop Q max
1682 delete his3;
1683 } // loop COG
1684 delete his4;
1685 }//loop local angle
1686 delete his5;
1687 }
1688}
7d14c1c1 1689
1690
1691int AliTPCcalibTracks::GetTHnStat(const THnBase *H, THnBase *&Mean, THnBase *&Sigma, THnBase *&Entr )
1692{
1693 // H is THnSparse:
1694 //
1695 // OBJ: TAxis var var
1696 // OBJ: TAxis axis 1
1697 // OBJ: TAxis axis 2
1698 // ...
1699
1700 // Output:
1701
1702 // Mean, Sigma and Entr is THnF of mean, RMS and entries of var:
1703 //
1704 // OBJ: TAxis axis 1
1705 // OBJ: TAxis axis 2
1706 // ..
1707
1708 Mean = 0;
1709 Sigma = 0;
1710 Entr = 0;
1711 if( !H ) return -1;
1712
1713 Bool_t ok = kTRUE;
1714
1715 int nDim = H->GetNdimensions();
1716 Long64_t nFilledBins = H->GetNbins();
1717 Long64_t nStatBins = 0;
1718
1719 Int_t *nBins = new Int_t [nDim];
1720 Double_t *xMin = new Double_t [nDim];
1721 Double_t *xMax = new Double_t [nDim];
1722 Int_t *idx = new Int_t [nDim];
1723
1724 TString nameMean = H->GetName();
1725 TString nameSigma = H->GetName();
1726 TString nameEntr = H->GetName();
1727 nameMean+="_Mean";
1728 nameSigma+="_Sigma";
1729 nameEntr+="_Entr";
1730
1731 ok = ok && ( nBins && xMin && xMax && idx );
1732
1733 if( ok ){
1734 for( int i=0; i<nDim; i++){
1735 xMin[i] = H->GetAxis(i)->GetXmin();
1736 xMax[i] = H->GetAxis(i)->GetXmax();
1737 nBins[i] = H->GetAxis(i)->GetNbins();
1738 }
1739
1740 Mean = new THnSparseF( nameMean.Data(), nameMean.Data(), nDim-1, nBins+1, xMin+1, xMax+1 );
1741 Sigma = new THnSparseF( nameSigma.Data(), nameSigma.Data(), nDim-1, nBins+1, xMin+1, xMax+1 );
1742 Entr = new THnSparseF( nameEntr.Data(), nameEntr.Data(), nDim-1, nBins+1, xMin+1, xMax+1 );
1743 }
1744
1745 ok = ok && ( Mean && Sigma && Entr );
1746
1747 if( ok ){
1748 for( int i=0; i<nDim-1; i++){
1749 const TAxis *ax = H->GetAxis(i+1);
1750 Mean->GetAxis(i)->SetName(ax->GetName());
1751 Mean->GetAxis(i)->SetTitle(ax->GetTitle());
1752 Sigma->GetAxis(i)->SetName(ax->GetName());
1753 Sigma->GetAxis(i)->SetTitle(ax->GetTitle());
1754 Entr->GetAxis(i)->SetName(ax->GetName());
1755 Entr->GetAxis(i)->SetTitle(ax->GetTitle());
1756 if( ax->GetXbins()->fN>0 ){
1757 Mean->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1758 Sigma->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1759 Entr->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1760 }
1761 }
1762
1763 // Allocate bins
1764
1765 for( Long64_t binS=0; binS<nFilledBins; binS++){
1766 H->GetBinContent(binS,idx);
1767 Mean->GetBin(idx+1,kTRUE);
1768 Sigma->GetBin(idx+1,kTRUE);
1769 Entr->GetBin(idx+1,kTRUE);
1770 }
1771
1772 nStatBins = Mean->GetNbins();
1773
1774 }
1775
1776
1777 Long64_t *hMap = new Long64_t[nFilledBins];
1778 Double_t *hVal = new Double_t[nFilledBins];
1779 Double_t *hEntr = new Double_t[nFilledBins];
1780 Double_t *meanD = new Double_t[nStatBins];
1781 Double_t *sigmaD = new Double_t[nStatBins];
1782
1783 ok = ok && ( hMap && hVal && hEntr && meanD && sigmaD );
1784
1785 // Map bins to Stat bins
1786
1787 if( ok ){
1788 for( Long64_t binS=0; binS<nFilledBins; binS++){
1789 double val = H->GetBinContent(binS,idx);
1790 Long64_t bin = Mean->GetBin(idx+1,kFALSE);
1791 ok = ok && ( bin>=0 && bin<nStatBins && bin==Sigma->GetBin(idx+1,kFALSE) && bin== Entr->GetBin(idx+1,kFALSE) );
1792 if( !ok ) break;
1793 if( val < 0 ){
1794 cout << "AliTPCcalibTracks: GetSparseStat : Unexpected content of the input THn"<<endl;
1795 bin = -1;
1796 }
1797 if( idx[0]<1 || idx[0]>nBins[0] ) bin = -1;
1798 hMap[binS] = bin;
1799 hVal[binS] = idx[0];
1800 hEntr[binS] = val;
1801 }
1802 }
1803
1804 // do 2 iteration with cut at 4 Sigma
1805
1806 for( int iter=0; ok && iter<2; iter++ ){
1807
1808 // clean up entries
1809
1810 for( Long64_t bin=0; bin < nStatBins; bin++){
1811 Entr->SetBinContent(bin, 0);
1812 meanD[bin]=0;
1813 sigmaD[bin]=0;
1814 }
1815
1816 // get Entries and Mean
1817
1818 for( Long64_t binS=0; binS<nFilledBins; binS++){
1819 Long64_t bin = hMap[binS];
1820 Double_t val = hVal[binS];
1821 Double_t entr = hEntr[binS];
1822 if( bin < 0 ) continue;
1823 if( iter!=0 ){ // cut
1824 double s2 = Sigma->GetBinContent(bin);
1825 double d = val - Mean->GetBinContent(bin);
1826 if( d*d>s2*16 ) continue;
1827 }
1828 meanD[bin]+= entr*val;
1829 Entr->AddBinContent(bin,entr);
1830 }
1831
1832 for( Long64_t bin=0; bin<nStatBins; bin++){
1833 double val = meanD[bin];
1834 double sum = Entr->GetBinContent(bin);
1835 if( sum>=1 ){
1836 Mean->SetBinContent(bin, val/sum );
1837 }
1838 else Mean->SetBinContent(bin, 0);
1839 Entr->SetBinContent(bin, 0);
1840 }
1841
1842 // get RMS
1843
1844 for( Long64_t binS=0; binS<nFilledBins; binS++){
1845 Long64_t bin = hMap[binS];
1846 Double_t val = hVal[binS];
1847 Double_t entr = hEntr[binS];
1848 if( bin < 0 ) continue;
1849 double d2 = val - Mean->GetBinContent(bin);
1850 d2 *= d2;
1851 if( iter!=0 ){ // cut
1852 double s2 = Sigma->GetBinContent(bin);
1853 if( d2>s2*16 ) continue;
1854 }
1855 sigmaD[bin] += entr*d2;
1856 Entr->AddBinContent(bin,entr);
1857 }
1858
1859 for( Long64_t bin=0; bin<nStatBins; bin++ ){
1860 double val = sigmaD[bin];
1861 double sum = Entr->GetBinContent(bin);
1862 if( sum>=1 && val>=0 ){
1863 Sigma->SetBinContent(bin, val/sum );
1864 }
1865 else Sigma->SetBinContent(bin, 0);
1866 }
1867 }
1868
1869 // scale the Mean and the Sigma
1870
1871 if( ok ){
1872 double v0 = H->GetAxis(0)->GetBinCenter(1);
1873 double vscale = H->GetAxis(0)->GetBinWidth(1);
1874
1875 for( Long64_t bin=0; bin<nStatBins; bin++){
1876 double m = Mean->GetBinContent(bin);
1877 double s2 = Sigma->GetBinContent(bin);
1878 double sum = Entr->GetBinContent(bin);
1879 if( sum>0 && s2>=0 ){
1880 Mean->SetBinContent(bin, v0 + (m-1)*vscale );
1881 Sigma->SetBinContent(bin, sqrt(s2)*vscale );
1882 }
1883 else{
1884 Mean->SetBinContent(bin, 0);
1885 Sigma->SetBinContent(bin, 0);
1886 Entr->SetBinContent(bin, 0);
1887 }
1888 }
1889 }
1890
d4fa6a95 1891 delete[] nBins;
1892 delete[] xMin;
1893 delete[] xMax;
1894 delete[] idx;
1895 delete[] hMap;
1896 delete[] hVal;
1897 delete[] hEntr;
1898 delete[] meanD;
1899 delete[] sigmaD;
7d14c1c1 1900
1901 if( !ok ){
1902 cout << "AliTPCcalibTracks: GetSparseStat : No memory or internal Error "<<endl;
1903 delete Mean;
1904 delete Sigma;
1905 delete Entr;
1906 Mean =0;
1907 Sigma = 0;
1908 Entr = 0;
1909 return -1;
1910 }
1911
1912 return 0;
1913}
1914
1915
1916
1917int AliTPCcalibTracks::CreateWaveCorrection(const THnBase *DeltaY, THnBase *&MeanY, THnBase *&SigmaY, THnBase *&EntrY,
1918 Bool_t MirrorZ, Bool_t MirrorPad, Bool_t MirrorAngle, Int_t MinStat )
1919{
1920 // DeltaY is THnSparse:
1921 //
1922 // OBJ: TAxis 0 var var
1923 // OBJ: TAxis 1 pad type pad type
1924 // OBJ: TAxis 2 z z
1925 // OBJ: TAxis 3 Qmax Qmax
1926 // OBJ: TAxis 4 cog cog
1927 // OBJ: TAxis 5 tan(angle) tan(angle)
1928 // ...
1929
1930 MeanY = 0;
1931 SigmaY = 0;
1932 EntrY = 0;
1933
1934 int nDim = DeltaY->GetNdimensions();
d4fa6a95 1935 if( nDim<6 ){
1936 cout << "AliTPCcalibTracks: CreateWaveCorrection: Unknown input"<<endl;
1937 return -1;
7d14c1c1 1938 }
1939
d4fa6a95 1940 Int_t *nBins = new Int_t[nDim];
1941 Int_t *nBinsNew = new Int_t[nDim];
1942 Double_t *xMin = new Double_t[nDim];
1943 Double_t *xMax = new Double_t[nDim];
1944 Int_t *idx = new Int_t[nDim];
1945 THnSparseD *mergedDeltaY = 0;
7d14c1c1 1946
d4fa6a95 1947 int ret = 0;
7d14c1c1 1948
d4fa6a95 1949 if( !nBins || !nBinsNew || !xMin || !xMax || !idx ){
1950 ret = -1;
1951 cout << "AliTPCcalibTracks: CreateWaveCorrection: Out of memory"<<endl;
7d14c1c1 1952 }
d4fa6a95 1953
1954 if( ret==0 ){
1955
1956 for( int i=0; i<nDim; i++){
1957 xMin[i] = DeltaY->GetAxis(i)->GetXmin();
1958 xMax[i] = DeltaY->GetAxis(i)->GetXmax();
1959 nBins[i] = DeltaY->GetAxis(i)->GetNbins();
1960 nBinsNew[i] = nBins[i];
1961 }
7d14c1c1 1962
d4fa6a95 1963 // Merge cog axis
1964 if( MirrorPad ){
1965 Int_t centralBin = DeltaY->GetAxis(4)->FindFixBin(0.5);
1966 xMin[4] = DeltaY->GetAxis(4)->GetBinLowEdge(centralBin);
1967 nBinsNew[4] = nBinsNew[4]-centralBin+1;
7d14c1c1 1968 }
7d14c1c1 1969
d4fa6a95 1970 // Merge Z axis
1971 if( MirrorZ ){
1972 Int_t centralBin = DeltaY->GetAxis(2)->FindFixBin(0.0);
1973 xMin[2] = DeltaY->GetAxis(2)->GetBinLowEdge(centralBin)-0.0;
1974 nBinsNew[2] = nBinsNew[2]-centralBin+1;
1975 }
7d14c1c1 1976
d4fa6a95 1977 // Merge Angle axis
1978 if( MirrorAngle ){
1979 Int_t centralBin = DeltaY->GetAxis(5)->FindFixBin(0.0);
1980 xMin[5] = DeltaY->GetAxis(5)->GetBinLowEdge(centralBin)-0.0;
1981 nBinsNew[5] = nBinsNew[5]-centralBin+1;
7d14c1c1 1982 }
1983
d4fa6a95 1984 // Merge Sparse
1985
1986 mergedDeltaY = new THnSparseD("mergedDeltaY", "mergedDeltaY",nDim, nBinsNew, xMin, xMax );
1987
1988 if( !mergedDeltaY ){
1989 cout << "AliTPCcalibTracks: CreateWaveCorrection: Can not copy a Sparse"<<endl;
1990 ret = -1;
7d14c1c1 1991 }
d4fa6a95 1992 }
1993
1994 if( ret == 0 ){
1995
1996 for( int i=0; i<nDim; i++){
1997 const TAxis *ax = DeltaY->GetAxis(i);
1998 mergedDeltaY->GetAxis(i)->SetName(ax->GetName());
1999 mergedDeltaY->GetAxis(i)->SetTitle(ax->GetTitle());
2000 if( ax->GetXbins()->fN>0 ){
2001 mergedDeltaY->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
2002 }
7d14c1c1 2003 }
d4fa6a95 2004
2005 for( Long64_t binS=0; binS<DeltaY->GetNbins(); binS++){
2006 Double_t stat = DeltaY->GetBinContent(binS,idx);
2007 if( stat < 1 ) continue;
2008 bool swap0=0;
2009
2010 if( MirrorPad && idx[4]>0){ // underflow reserved for contains one-pad clusters, don't mirror
2011 Double_t v = DeltaY->GetAxis(4)->GetBinCenter(idx[4]);
2012 if( v < 0.5 ) swap0 = !swap0;
2013 idx[4] = mergedDeltaY->GetAxis(4)->FindFixBin( 0.5 + TMath::Abs(0.5 - v) );
2014 }
2015
2016 if( MirrorZ ){
2017 Double_t v = DeltaY->GetAxis(2)->GetBinCenter(idx[2]);
2018 if( v < 0.0 ) swap0 = !swap0;
2019 if( idx[2]<=0 ) idx[2] = nBinsNew[2]+1;
2020 else idx[2] = mergedDeltaY->GetAxis(2)->FindFixBin( TMath::Abs(v) );
2021 }
2022
2023 if( MirrorAngle ){
2024 Double_t v = DeltaY->GetAxis(5)->GetBinCenter(idx[5]);
2025 if( idx[5]<=0 ) idx[5] = nBinsNew[5]+1;
2026 else idx[5] = mergedDeltaY->GetAxis(5)->FindFixBin( TMath::Abs(v) );
2027 }
7d14c1c1 2028
d4fa6a95 2029 if( swap0 ){
2030 if( idx[0]<=0 ) idx[0] = nBinsNew[0]+1;
2031 else if( idx[0] >= nBins[0]+1 ) idx[0] = 0;
2032 else {
2033 Double_t v = DeltaY->GetAxis(0)->GetBinCenter(idx[0]);
2034 idx[0] = mergedDeltaY->GetAxis(0)->FindFixBin(-v);
2035 }
7d14c1c1 2036 }
d4fa6a95 2037
2038 Long64_t bin = mergedDeltaY->GetBin(idx,kTRUE);
2039 if( bin<0 ){
2040 cout << "AliTPCcalibTracks: CreateWaveCorrection : wrong bining"<<endl;
2041 continue;
2042 }
2043 mergedDeltaY->AddBinContent(bin,stat);
7d14c1c1 2044 }
2045
d4fa6a95 2046 ret = GetTHnStat(mergedDeltaY, MeanY, SigmaY, EntrY );
7d14c1c1 2047 }
2048
d4fa6a95 2049 if( ret==0 ){
2050
2051 MeanY->SetName("TPCWaveCorrectionMap");
2052 MeanY->SetTitle("TPCWaveCorrectionMap");
2053 SigmaY->SetName("TPCResolutionMap");
2054 SigmaY->SetTitle("TPCResolutionMap");
2055 EntrY->SetName("TPCWaveCorrectionEntr");
2056 EntrY->SetTitle("TPCWaveCorrectionEntr");
7d14c1c1 2057
d4fa6a95 2058 for( Long64_t bin=0; bin<MeanY->GetNbins(); bin++){
2059 Double_t stat = EntrY->GetBinContent(bin,idx);
2060
2061 // Normalize : Set no correction for one pad clusters
2062
2063 if( idx[3]<=0 ) MeanY->SetBinContent(bin,0);
2064
2065 // Suppress bins with low statistic
2066
2067 if( stat<MinStat ){
2068 EntrY->SetBinContent(bin,0);
2069 MeanY->SetBinContent(bin,0);
2070 SigmaY->SetBinContent(bin,-1);
2071 }
2072 }
2073
2074 }
7d14c1c1 2075
d4fa6a95 2076 delete[] nBins;
2077 delete[] nBinsNew;
2078 delete[] xMin;
2079 delete[] xMax;
2080 delete[] idx;
2081 delete mergedDeltaY;
7d14c1c1 2082
d4fa6a95 2083 if( ret!=0 ){
2084 delete MeanY;
2085 delete SigmaY;
2086 delete EntrY;
2087 MeanY = 0;
2088 SigmaY = 0;
2089 EntrY = 0;
7d14c1c1 2090 }
7d14c1c1 2091
d4fa6a95 2092 return ret;
7d14c1c1 2093}
2094
2095int AliTPCcalibTracks::UpdateClusterParam( AliTPCClusterParam *cParam, Bool_t MirrorZ, Bool_t MirrorPad, Bool_t MirrorAngle, Int_t MinStat )
2096{
2097 if( !cParam || !fHisDeltaY) return -1;
2098
2099 THnBase *meanY = 0;
2100 THnBase *sigmaY = 0;
2101 THnBase *entrY = 0;
2102 int ret = CreateWaveCorrection(fHisDeltaY, meanY, sigmaY, entrY, MirrorZ, MirrorAngle, MirrorPad, MinStat );
2103 if( ret<0 ) return ret;
2104 cParam->SetWaveCorrectionMap(meanY);
2105 cParam->SetResolutionYMap(sigmaY);
2106 delete meanY;
2107 delete sigmaY;
2108 delete entrY;
2109 return 0;
2110}