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