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