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