]> 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 display of calibration type in CalibViewer (Jens)
[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));
663 // ---------------------------------------------------------------------
664 //
665 //
10757ee9 666
46e89793 667 for (Int_t irow = kDelta; irow < 159-kDelta; irow++){
668 // loop again over all rows along the track
669 // do analysis
670 //
671 Int_t nclFound = 0; // number of clusters in the neighborhood
672 Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster
673 Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster
674 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
675 if (!cluster0) continue;
676 Int_t sector = cluster0->GetDetector();
677 Float_t xref = cluster0->GetX();
678
679 // Make Fit
680 fFitterParY.ClearPoints();
681 fFitterParZ.ClearPoints();
682 fFitterParYWeight.ClearPoints();
683 fFitterParZWeight.ClearPoints();
684 // fit tracklet (clusters in given padrow +- kDelta padrows)
685 // with polynom of 2nd order and two polynoms of 1st order
686 // take both polynoms of 1st order, calculate difference of their parameters
687 // add covariance matrixes and calculate chi2 of this difference
688 // if this chi2 is bigger than a given threshold, assume that the current cluster is
689 // a kink an goto next padrow
7d14c1c1 690
691 // first fit the track angle for wave correction
692
693 AliRieman riemanFitAngle( 2*kDelta ); //SG
694
695 if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){
696 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
697 // loop over irow +- kDelta rows (neighboured rows)
698 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
699 AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
700 if (!currentCluster) continue;
701 if (currentCluster->GetType() < 0) continue;
702 if (currentCluster->GetDetector() != sector) continue;
703 riemanFitAngle.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ);
704 }
705 if( riemanFitAngle.GetN()>3 ) riemanFitAngle.Update();
706 }
707
708 // do fit
709
46e89793 710 AliRieman riemanFit(2*kDelta);
711 AliRieman riemanFitW(2*kDelta);
712 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
713 // loop over irow +- kDelta rows (neighboured rows)
714 //
715 //
716 if (idelta == 0) continue; // don't use center cluster
717 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
718 AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
719 if (!currentCluster) continue;
720 if (currentCluster->GetType() < 0) continue;
721 if (currentCluster->GetDetector() != sector) continue;
722 nclFound++;
723 if (idelta < 0){
724 ncl0++;
9f0beaf7 725 }
46e89793 726 if (idelta > 0){
727 ncl1++;
10757ee9 728 }
7d14c1c1 729 //SG!!!
730 double dY = 0;
731 if( fClusterParam ){
732 Int_t padSize = 0; // short pads
733 if (currentCluster->GetDetector() >= 36) {
734 padSize = 1; // medium pads
735 if (currentCluster->GetRow() > 63) padSize = 2; // long pads
736 }
737 dY = fClusterParam->GetWaveCorrection( padSize,
738 currentCluster->GetZ(),
739 currentCluster->GetMax(),
740 currentCluster->GetPad(),
741 riemanFitAngle.GetDYat(currentCluster->GetX())
742 );
743 }
744 riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), csigmaY,csigmaZ);
745 riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), csigmaY*TMath::Sqrt(1+TMath::Abs(idelta)),csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta)));
46e89793 746 } // loop over neighbourhood for fitter filling
747 if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
748 if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
749 if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow
750 riemanFit.Update();
751 riemanFitW.Update();
752 Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5));
753 Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5));
754 Double_t paramYR[3], paramZR[3];
755 Double_t paramYRW[3], paramZRW[3];
756 //
757 paramYR[0] = riemanFit.GetYat(xref);
758 paramZR[0] = riemanFit.GetZat(xref);
759 paramYRW[0] = riemanFitW.GetYat(xref);
760 paramZRW[0] = riemanFitW.GetZat(xref);
761 //
762 paramYR[1] = riemanFit.GetDYat(xref);
763 paramZR[1] = riemanFit.GetDZat(xref);
764 paramYRW[1] = riemanFitW.GetDYat(xref);
765 paramZRW[1] = riemanFitW.GetDZat(xref);
766 //
767 Int_t reject=0;
768 if (chi2R>kChi2Cut) reject+=1;
769 if (chi2RW>kChi2Cut) reject+=2;
770 if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4;
771 if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8;
772 if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16;
773 if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32;
774 //
775 TTreeSRedirector *cstream = GetDebugStreamer();
776 // get fit parameters from pol2 fit:
777 Double_t tracky = paramYR[0];
778 Double_t trackz = paramZR[0];
779 Float_t deltay = cluster0->GetY()-tracky;
780 Float_t deltaz = cluster0->GetZ()-trackz;
781 Float_t angley = paramYR[1];
782 Float_t anglez = paramZR[1];
783
784 Int_t padSize = 0; // short pads
785 if (cluster0->GetDetector() >= 36) {
10757ee9 786 padSize = 1; // medium pads
787 if (cluster0->GetRow() > 63) padSize = 2; // long pads
46e89793 788 }
789 Int_t ipad = TMath::Nint(cluster0->GetPad());
790 if (cstream){
791 Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
792 (*cstream)<<"Resol0"<<
793 "run="<<fRun<< // run number
794 "event="<<fEvent<< // event number
795 "time="<<fTime<< // time stamp of event
796 "trigger="<<fTrigger<< // trigger
797 "mag="<<fMagF<< // magnetic field
798 "padSize="<<padSize<<
799 //
800 "reject="<<reject<<
801 "cl.="<<cluster0<< // cluster info
802 "xref="<<xref<< // cluster reference X
803 //rieman fit
804 "yr="<<paramYR[0]<< // track position y - rieman fit
805 "zr="<<paramZR[0]<< // track position z - rieman fit
806 "yrW="<<paramYRW[0]<< // track position y - rieman fit - weight
807 "zrW="<<paramZRW[0]<< // track position z - rieman fit - weight
808 "dyr="<<paramYR[1]<< // track position y - rieman fit
809 "dzr="<<paramZR[1]<< // track position z - rieman fit
810 "dyrW="<<paramYRW[1]<< // track position y - rieman fit - weight
811 "dzrW="<<paramZRW[1]<< // track position z - rieman fit - weight
812 //
813 "angley="<<angley<< // angle par fit
814 "anglez="<<anglez<< // angle par fit
815 "zdr="<<zdrift<< //
816 "dy="<<deltay<<
817 "dz="<<deltaz<<
818 "sy="<<csigmaY<<
819 "sz="<<csigmaZ<<
820 "chi2R="<<chi2R<<
821 "chi2RW="<<chi2RW<<
822 "\n";
823 }
824 if (reject) continue;
825
826
827 // =========================================
828 // wirte collected information to histograms
829 // =========================================
830
831 Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
832 fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
833
834
835 TH3F * his3 = 0;
836 his3 = (TH3F*)fRMSY->At(padSize);
7e22eccb 837 if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
46e89793 838 his3 = (TH3F*)fRMSZ->At(padSize);
7e22eccb 839 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
46e89793 840
841 his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
7e22eccb 842 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
46e89793 843 his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
7e22eccb 844 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
46e89793 845
846
847 his3 = (TH3F*)fResolY->At(padSize);
848 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
849 his3 = (TH3F*)fResolZ->At(padSize);
850 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
851 his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
852 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
853 his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
854 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
855 //=============================================================================================
856 //
857 // Fill THN histograms
858 //
859 Double_t xvar[9];
60e31b96 860 xvar[1]=padSize; // pad type
861 xvar[2]=cluster0->GetZ(); //
46e89793 862 xvar[3]=cluster0->GetMax();
863
864 xvar[0]=deltay;
7d14c1c1 865 double cog = cluster0->GetPad()-Int_t(cluster0->GetPad());// distance to the center of the pad
866 xvar[4] = cog;
867 xvar[5]=angley;
868
869 if( TMath::Abs(cog-0.5)<1.e-8 ) xvar[4]=-1; // fill one-pad clusters in the underflow bin
870 fHisDeltaY->Fill(xvar);
871
872 xvar[4]=cog;
46e89793 873 xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
874 fHisRMSY->Fill(xvar);
875
876 xvar[0]=deltaz;
60e31b96 877 xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin()); // distance to the center of the time bin
46e89793 878 xvar[5]=anglez;
879 fHisDeltaZ->Fill(xvar);
880 xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
881 fHisRMSZ->Fill(xvar);
882
883 } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
884} // FillResolutionHistoLocal(...)
885
10757ee9 886
887
888
889
890
10757ee9 891
892
893void AliTPCcalibTracks::SetStyle() const {
894 //
895 // set style, can be called by all draw functions
896 //
897 gROOT->SetStyle("Plain");
898 gStyle->SetFillColor(10);
899 gStyle->SetPadColor(10);
900 gStyle->SetCanvasColor(10);
901 gStyle->SetStatColor(10);
902 gStyle->SetPalette(1,0);
903 gStyle->SetNumberContours(60);
904}
905
906
10757ee9 907
f78da5ae 908void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
10757ee9 909 //
910 // all functions are called, that produce pictures
911 // the histograms are written to the directory 'pathName'
912 // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
913 // 'stat' is also the number of minEntries for MakeResPlotsQTree
914 //
915
ae28e92e 916 if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
46e89793 917 MakeResPlotsQTree(stat, pathName);
10757ee9 918}
919
920
10757ee9 921
922
f78da5ae 923void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
10757ee9 924 //
925 // Make tree with resolution parameters
926 // the result is written to 'resol.root' in directory 'pathname'
927 // file information are available in fileInfo
928 // available variables in the tree 'Resol':
929 // Entries: number of entries for this resolution point
930 // nbins: number of bins that were accumulated
931 // Dim: direction, Dim==0: y-direction, Dim==1: z-direction
932 // Pad: padSize; short, medium and long
933 // Length: pad length, 0.75, 1, 1.5
934 // QMean: mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
935 // Zc: center of middle bin in drift direction
936 // Zm: mean dirftlength for accumulated Delta-Histograms
937 // Zs: width of driftlength bin
938 // AngleC: center of middle bin in Angle-Direction
939 // AngleM: mean angle for accumulated Delta-Histograms
940 // AngleS: width of Angle-bin
941 // Resol: sigma for gaus fit through Delta-Histograms
942 // Sigma: error of sigma for gaus fit through Delta Histograms
943 // MeanR: mean of the Delta-Histogram
944 // SigmaR: rms of the Delta-Histogram
945 // RMSm: mean of the gaus fit through RMS-Histogram
946 // RMS: sigma of the gaus fit through RMS-Histogram
947 // RMSe0: error of mean of gaus fit in RMS-Histogram
948 // RMSe1: error of sigma of gaus fit in RMS-Histogram
949 //
950
ae28e92e 951 if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
952 if (GetDebugLevel() > -1) cout << " relax, the calculation will take a while..." << endl;
10757ee9 953
954 gSystem->MakeDirectory(pathName);
955 gSystem->ChangeDirectory(pathName);
956 TString kFileName = "resol.root";
957 TTreeSRedirector fTreeResol(kFileName.Data());
958
959 TH3F *resArray[2][3][11];
960 TH3F *rmsArray[2][3][11];
961
962 // load histograms from fArrayQDY and fArrayQDZ
963 // into resArray and rmsArray
964 // that is all we need here
965 for (Int_t idim = 0; idim < 2; idim++){
966 for (Int_t ipad = 0; ipad < 3; ipad++){
967 for (Int_t iq = 0; iq <= 10; iq++){
968 rmsArray[idim][ipad][iq]=0;
969 resArray[idim][ipad][iq]=0;
970 Int_t bin = GetBin(iq,ipad);
971 TH3F *hresl = 0;
972 if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
973 if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
974 if (!hresl) continue;
975 resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
976 resArray[idim][ipad][iq]->SetDirectory(0);
977 TH3F * hreslRMS = 0;
978 if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
979 if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
980 if (!hreslRMS) continue;
981 rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
982 rmsArray[idim][ipad][iq]->SetDirectory(0);
983 }
984 }
985 }
ae28e92e 986 if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
10757ee9 987
988 //--------------------------------------------------------------------------------------------
989
990 char name[200];
991 Double_t qMean;
992 Double_t zMean, angleMean, zCenter, angleCenter;
993 Double_t zSigma, angleSigma;
994 TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
995 TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
996 TF1 *fitFunction = new TF1("fitFunction", "gaus");
997 Float_t entriesQ = 0;
998 Int_t loopCounter = 1;
999
1000 for (Int_t idim = 0; idim < 2; idim++){
1001 // Loop y-z corrdinate
1002 for (Int_t ipad = 0; ipad < 3; ipad++){
1003 // loop pad type
1004 for (Int_t iq = -1; iq < 10; iq++){
1005 // LOOP Q
ae28e92e 1006 if (GetDebugLevel() > -1)
10757ee9 1007 cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
1008 << (Int_t)((loopCounter)/66.*100) << "% done), "
1009 << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << " \r" << std::flush;
1010 loopCounter++;
1011 TH3F *hres = 0;
1012 TH3F *hrms = 0;
1013 qMean = 0;
1014 entriesQ = 0;
1015
1016 // calculate qMean
1017 if (iq == -1){
1018 // integrated spectra
1019 for (Int_t iql = 0; iql < 10; iql++){
1020 Int_t bin = GetBin(iql,ipad);
1021 TH3F *hresl = resArray[idim][ipad][iql];
1022 TH3F *hrmsl = rmsArray[idim][ipad][iql];
1023 if (!hresl) continue;
1024 if (!hrmsl) continue;
1025 entriesQ += hresl->GetEntries();
1026 qMean += hresl->GetEntries() * GetQ(bin);
1027 if (!hres) {
1028 hres = (TH3F*)hresl->Clone();
1029 hrms = (TH3F*)hrmsl->Clone();
1030 }
1031 else{
1032 hres->Add(hresl);
1033 hrms->Add(hrmsl);
1034 }
1035 }
1036 qMean /= entriesQ;
1037 qMean *= -1.; // integral mean charge
1038 }
1039 else {
1040 // loop over neighboured Q-bins
1041 // accumulate entries from neighboured Q-bins
1042 for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
1043 if (iql < 0) continue;
1044 Int_t bin = GetBin(iql,ipad);
1045 TH3F * hresl = resArray[idim][ipad][iql];
1046 TH3F * hrmsl = rmsArray[idim][ipad][iql];
1047 if (!hresl) continue;
1048 if (!hrmsl) continue;
1049 entriesQ += hresl->GetEntries();
1050 qMean += hresl->GetEntries() * GetQ(bin);
1051 if (!hres) {
1052 hres = (TH3F*) hresl->Clone();
1053 hrms = (TH3F*) hrmsl->Clone();
1054 }
1055 else{
1056 hres->Add(hresl);
1057 hrms->Add(hrmsl);
1058 }
1059 }
1060 qMean/=entriesQ;
1061 }
f69a0aa7 1062 if (!hres) continue;
1063 if (!hrms) continue;
1064
10757ee9 1065 TAxis *xAxisDriftLength = hres->GetXaxis(); // driftlength / z - axis
1066 TAxis *yAxisAngle = hres->GetYaxis(); // angle axis
1067 TAxis *zAxisDelta = hres->GetZaxis(); // delta axis
1068 TAxis *zAxisRms = hrms->GetZaxis(); // rms axis
1069
1070 // loop over all angle bins
1071 for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1072 angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1073 // loop over all driftlength bins
1074 for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1075 zCenter = xAxisDriftLength->GetBinCenter(ibinxDL);
1076 zSigma = xAxisDriftLength->GetBinWidth(ibinxDL);
1077 angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
1078 zMean = zCenter; // changens, when more statistic is accumulated
1079 angleMean = angleCenter; // changens, when more statistic is accumulated
1080
1081 // create 2 1D-Histograms, projectionRes and projectionRms
1082 // these histograms are delta histograms for given direction, padSize, chargeBin,
1083 // angleBin and driftLengthBin
1084 // later on they will be fitted with a gausian, its sigma is the resoltuion...
4aa37f93 1085 snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
10757ee9 1086 // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1087 projectionRes->SetNameTitle(name, name);
4aa37f93 1088 snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
10757ee9 1089 // TH1D * projectionRms = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1090 projectionRms->SetNameTitle(name, name);
1091
1092 projectionRes->Reset();
1093 projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1094 projectionRms->Reset();
1095 projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1096 projectionRes->SetDirectory(0);
1097 projectionRms->SetDirectory(0);
1098
1099 Double_t entries = 0;
1100 Int_t nbins = 0; // counts, how many bins were accumulated
1101 zMean = 0;
1102 angleMean = 0;
1103 entries = 0;
1104
1105 // fill projectionRes and projectionRms for given dim, ipad and iq,
1106 // as well as for given angleBin and driftlengthBin
1107 // if this gives not enough statistic, include neighbourhood
1108 // (angle and driftlength) successifely
1109 for (Int_t dbin = 0; dbin <= 8; dbin++){ // delta-bins around centered angleBin and driftlengthBin
1110 for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) { // delta-bins in angle direction
1111 for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
1112 if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; // add each bin only one time !
1113 Int_t binx2 = ibinxDL + dbinx2; // position variable in x (driftlength) direction
1114 Int_t biny2 = ibinyAngle + dbiny2; // position variable in y (angle) direction
1115 if (binx2 < 1 || biny2 < 1) continue; // don't go out of the histogram!
1116 if (binx2 >= xAxisDriftLength->GetNbins()) continue; // don't go out of the histogram!
1117 if (biny2 >= yAxisAngle->GetNbins()) continue; // don't go out of the histogram!
1118 nbins++; // count the number of accumulated bins
1119 // Fill resolution histo
1120 for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
1121 // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3); // unused variable
1122 projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
1123 entries += hres->GetBinContent(binx2, biny2, ibin3);
1124 zMean += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
1125 angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
1126 } // ibin3 loop
1127 // fill RMS histo
1128 for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
1129 projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
1130 }
1131 } //dbinx2 loop
1132 if (entries > minEntries) break; // enough statistic accumulated
1133 } // dbiny2 loop
1134 if (entries > minEntries) break; // enough statistic accumulated
1135 } // dbin loop
1136 if ( entries< minEntries) continue; // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
1137 zMean /= entries;
1138 angleMean /= entries;
1139
1140 if (entries > minEntries) {
1141 // when enough statistic is accumulated
1142 // fit Delta histograms with a gausian
1143 // of the gausian is the resolution (resol), its fit error is sigma
1144 // store also mean and RMS of the histogram
1145 Float_t xmin = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1146 Float_t xmax = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1147
1148// projectionRes->Fit("gaus", "q0", "", xmin, xmax);
1149// Float_t resol = projectionRes->GetFunction("gaus")->GetParameter(2);
1150// Float_t sigma = projectionRes->GetFunction("gaus")->GetParError(2);
1151 fitFunction->SetMaximum(xmax);
1152 fitFunction->SetMinimum(xmin);
1153 projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
1154 Float_t resol = fitFunction->GetParameter(2);
1155 Float_t sigma = fitFunction->GetParError(2);
1156
1157 Float_t meanR = projectionRes->GetMean();
1158 Float_t sigmaR = projectionRes->GetRMS();
1159 // fit also RMS histograms with a gausian
1160 // store mean and sigma of the gausian in rmsMean and rmsSigma
1161 // store also the fit errors in errorRMS and errorSigma
1162 xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1163 xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1164
1165// projectionRms->Fit("gaus","q0","",xmin,xmax);
1166// Float_t rmsMean = projectionRms->GetFunction("gaus")->GetParameter(1);
1167// Float_t rmsSigma = projectionRms->GetFunction("gaus")->GetParameter(2);
1168// Float_t errorRMS = projectionRms->GetFunction("gaus")->GetParError(1);
1169// Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
1170 projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
1171 Float_t rmsMean = fitFunction->GetParameter(1);
1172 Float_t rmsSigma = fitFunction->GetParameter(2);
1173 Float_t errorRMS = fitFunction->GetParError(1);
1174 Float_t errorSigma = fitFunction->GetParError(2);
1175
1176 Float_t length = 0.75;
1177 if (ipad == 1) length = 1;
1178 if (ipad == 2) length = 1.5;
1179
1180 fTreeResol<<"Resol"<<
1181 "Entries="<<entries<< // number of entries for this resolution point
1182 "nbins="<<nbins<< // number of bins that were accumulated
1183 "Dim="<<idim<< // direction, Dim==0: y-direction, Dim==1: z-direction
1184 "Pad="<<ipad<< // padSize; short, medium and long
1185 "Length="<<length<< // pad length, 0.75, 1, 1.5
1186 "QMean="<<qMean<< // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1187 "Zc="<<zCenter<< // center of middle bin in drift direction
1188 "Zm="<<zMean<< // mean dirftlength for accumulated Delta-Histograms
1189 "Zs="<<zSigma<< // width of driftlength bin
1190 "AngleC="<<angleCenter<< // center of middle bin in Angle-Direction
1191 "AngleM="<<angleMean<< // mean angle for accumulated Delta-Histograms
1192 "AngleS="<<angleSigma<< // width of Angle-bin
1193 "Resol="<<resol<< // sigma for gaus fit through Delta-Histograms
1194 "Sigma="<<sigma<< // error of sigma for gaus fit through Delta Histograms
1195 "MeanR="<<meanR<< // mean of the Delta-Histogram
1196 "SigmaR="<<sigmaR<< // rms of the Delta-Histogram
1197 "RMSm="<<rmsMean<< // mean of the gaus fit through RMS-Histogram
1198 "RMSs="<<rmsSigma<< // sigma of the gaus fit through RMS-Histogram
1199 "RMSe0="<<errorRMS<< // error of mean of gaus fit in RMS-Histogram
1200 "RMSe1="<<errorSigma<< // error of sigma of gaus fit in RMS-Histogram
1201 "\n";
ae28e92e 1202 if (GetDebugLevel() > 5) {
10757ee9 1203 projectionRes->SetDirectory(fTreeResol.GetFile());
1204 projectionRes->Write(projectionRes->GetName());
1205 projectionRes->SetDirectory(0);
1206 projectionRms->SetDirectory(fTreeResol.GetFile());
1207 projectionRms->Write(projectionRms->GetName());
1208 projectionRes->SetDirectory(0);
1209 }
1210 } // if (projectionRes->GetSum() > minEntries)
1211 } // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
1212 } // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
1213
1214 } // iq-loop
1215 } // ipad-loop
1216 } // idim-loop
1217 delete projectionRes;
1218 delete projectionRms;
1219
1220// TFile resolFile(fTreeResol.GetFile());
1221 TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
1222 fileInfo.Write("fileInfo");
1223// resolFile.Close();
1224// fTreeResol.GetFile()->Close();
ae28e92e 1225 if (GetDebugLevel() > -1) cout << endl;
1226 if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
10757ee9 1227 gSystem->ChangeDirectory("..");
1228}
1229
1230
10757ee9 1231
10757ee9 1232
1233
1234Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
1235 //
1236 // function to merge several AliTPCcalibTracks objects after PROOF calculation
1237 // The object's histograms are merged via their merge functions
1238 // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
1239 //
1240
ae28e92e 1241 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl;
10757ee9 1242 if (!collectionList) return 0;
1243 if (collectionList->IsEmpty()) return -1;
1244
ae28e92e 1245 if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
1246 if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
10757ee9 1247 collectionList->Print();
1248
1249 // create a list for each data member
1250 TList* deltaYList = new TList;
1251 TList* deltaZList = new TList;
1252 TList* arrayAmpRowList = new TList;
1253 TList* rejectedTracksList = new TList;
10757ee9 1254 TList* clusterCutHistoList = new TList;
1255 TList* arrayAmpList = new TList;
1256 TList* arrayQDYList = new TList;
1257 TList* arrayQDZList = new TList;
1258 TList* arrayQRMSYList = new TList;
1259 TList* arrayQRMSZList = new TList;
10757ee9 1260 TList* resolYList = new TList;
1261 TList* resolZList = new TList;
1262 TList* rMSYList = new TList;
1263 TList* rMSZList = new TList;
1264
1265// TList* nRowsList = new TList;
1266// TList* nSectList = new TList;
1267// TList* fileNoList = new TList;
1268
1269 TIterator *listIterator = collectionList->MakeIterator();
1270 AliTPCcalibTracks *calibTracks = 0;
ae28e92e 1271 if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
10757ee9 1272 Int_t counter = 0;
a4c5fca9 1273 while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
10757ee9 1274 // loop over all entries in the collectionList and get dataMembers into lists
10757ee9 1275
10757ee9 1276 arrayQDYList->Add(calibTracks->GetfArrayQDY());
1277 arrayQDZList->Add(calibTracks->GetfArrayQDZ());
1278 arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
1279 arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
1280 resolYList->Add(calibTracks->GetfResolY());
1281 resolZList->Add(calibTracks->GetfResolZ());
1282 rMSYList->Add(calibTracks->GetfRMSY());
1283 rMSZList->Add(calibTracks->GetfRMSZ());
10757ee9 1284 rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
1285 clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
a4c5fca9 1286 //
1287 if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
1288 fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
1289 // fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
10757ee9 1290 counter++;
ae28e92e 1291 if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
af6a50bb 1292 AddHistos(calibTracks);
10757ee9 1293 }
1294
10757ee9 1295
1296 // merge data members
ae28e92e 1297 if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
10757ee9 1298 fClusterCutHisto->Merge(clusterCutHistoList);
1299 fRejectedTracksHisto->Merge(rejectedTracksList);
10757ee9 1300
1301 TObjArray* objarray = 0;
1302 TH1* hist = 0;
1303 TList* histList = 0;
1304 TIterator *objListIterator = 0;
1305
46e89793 1306
ae28e92e 1307 if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
10757ee9 1308 // merge fArrayQDY
1309 for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
1310 objListIterator = arrayQDYList->MakeIterator();
1311 histList = new TList;
1312 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1313 // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1314 hist = (TH3F*)(objarray->At(i));
1315 histList->Add(hist);
1316 }
1317 ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
1318 delete histList;
1319 delete objListIterator;
1320 }
1321
ae28e92e 1322 if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
10757ee9 1323 // merge fArrayQDZ
1324 for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1325 objListIterator = arrayQDZList->MakeIterator();
1326 histList = new TList;
1327 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1328 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1329 hist = (TH3F*)(objarray->At(i));
1330 histList->Add(hist);
1331 }
1332 ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
1333 delete histList;
1334 delete objListIterator;
1335 }
1336
ae28e92e 1337 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
10757ee9 1338 // merge fArrayQRMSY
1339 for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
1340 objListIterator = arrayQRMSYList->MakeIterator();
1341 histList = new TList;
1342 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1343 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
03533457 1344 // if (!objarray) continue; // removed for coverity -> JMT
10757ee9 1345 hist = (TH3F*)(objarray->At(i));
1346 histList->Add(hist);
1347 }
1348 ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
1349 delete histList;
1350 delete objListIterator;
1351 }
1352
ae28e92e 1353 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
10757ee9 1354 // merge fArrayQRMSZ
1355 for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1356 objListIterator = arrayQRMSZList->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
10757ee9 1360 hist = (TH3F*)(objarray->At(i));
1361 histList->Add(hist);
1362 }
1363 ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
1364 delete histList;
1365 delete objListIterator;
1366 }
1367
10757ee9 1368
10757ee9 1369
1370
1371
1372
ae28e92e 1373 if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
10757ee9 1374 // merge fResolY
1375 for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
1376 objListIterator = resolYList->MakeIterator();
1377 histList = new TList;
1378 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1379 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1380 hist = (TH3F*)(objarray->At(i));
1381 histList->Add(hist);
1382 }
1383 ((TH3F*)(fResolY->At(i)))->Merge(histList);
1384 delete histList;
1385 delete objListIterator;
1386 }
1387
1388 // merge fResolZ
1389 for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1390 objListIterator = resolZList->MakeIterator();
1391 histList = new TList;
1392 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1393 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1394 hist = (TH3F*)(objarray->At(i));
1395 histList->Add(hist);
1396 }
1397 ((TH3F*)(fResolZ->At(i)))->Merge(histList);
1398 delete histList;
1399 delete objListIterator;
1400 }
1401
1402 // merge fRMSY
1403 for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
1404 objListIterator = rMSYList->MakeIterator();
1405 histList = new TList;
1406 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1407 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1408 hist = (TH3F*)(objarray->At(i));
1409 histList->Add(hist);
1410 }
1411 ((TH3F*)(fRMSY->At(i)))->Merge(histList);
1412 delete histList;
1413 delete objListIterator;
1414 }
1415
1416 // merge fRMSZ
1417 for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1418 objListIterator = rMSZList->MakeIterator();
1419 histList = new TList;
1420 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1421 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1422 hist = (TH3F*)(objarray->At(i));
1423 histList->Add(hist);
1424 }
1425 ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
1426 delete histList;
1427 delete objListIterator;
1428 }
1429
1430 delete deltaYList;
1431 delete deltaZList;
1432 delete arrayAmpRowList;
1433 delete arrayAmpList;
1434 delete arrayQDYList;
1435 delete arrayQDZList;
1436 delete arrayQRMSYList;
1437 delete arrayQRMSZList;
1438 delete resolYList;
1439 delete resolZList;
1440 delete rMSYList;
1441 delete rMSZList;
10757ee9 1442 delete listIterator;
1443
ae28e92e 1444 if (GetDebugLevel() > 0) cout << "merging done!" << endl;
10757ee9 1445
1446 return 1;
1447}
1448
1449
10757ee9 1450
1451
af6a50bb 1452void AliTPCcalibTracks::MakeHistos(){
1453 //
1454 ////make THnSparse
1455 //
1456 //THnSparse *fHisDeltaY; // THnSparse - delta Y
1457 //THnSparse *fHisDeltaZ; // THnSparse - delta Z
1458 //THnSparse *fHisRMSY; // THnSparse - rms Y
1459 //THnSparse *fHisRMSZ; // THnSparse - rms Z
1460 //THnSparse *fHisQmax; // THnSparse - qmax
1461 //THnSparse *fHisQtot; // THnSparse - qtot
1462 // cluster performance bins
1463 // 0 - variable of interest
1464 // 1 - pad type - 0- short 1-medium 2-long pads
1465 // 2 - drift length - drift length -0-1
1466 // 3 - Qmax - Qmax - 2- 400
1467 // 4 - cog - COG position - 0-1
1468 // 5 - tan(phi) - local y angle
1469 // 6 - tan(theta) - local z angle
1470 // 7 - sector - sector number
1471 Double_t xminTrack[8], xmaxTrack[8];
1472 Int_t binsTrack[8];
1473 TString axisName[8];
1474
1475 //
46e89793 1476 binsTrack[0]=200;
af6a50bb 1477 axisName[0] ="var";
1478
1479 binsTrack[1] =3;
1480 xminTrack[1] =0; xmaxTrack[1]=3;
1481 axisName[1] ="pad type";
1482 //
46e89793 1483 binsTrack[2] =20;
1484 xminTrack[2] =-250; xmaxTrack[2]=250;
1485 axisName[2] ="z";
af6a50bb 1486 //
1487 binsTrack[3] =10;
1488 xminTrack[3] =1; xmaxTrack[3]=400;
1489 axisName[3] ="Qmax";
1490 //
46e89793 1491 binsTrack[4] =20;
af6a50bb 1492 xminTrack[4] =0; xmaxTrack[4]=1;
1493 axisName[4] ="cog";
1494 //
46e89793 1495 binsTrack[5] =15;
1496 xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
1497 axisName[5] ="tan(angle)";
af6a50bb 1498 //
af6a50bb 1499 //
46e89793 1500 xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1501 fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1502 xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1503 fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1504 xminTrack[0] =0.; xmaxTrack[0]=0.5;
46e89793 1505 fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1506 xminTrack[0] =0.; xmaxTrack[0]=0.5;
46e89793 1507 fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1508 xminTrack[0] =0.; xmaxTrack[0]=100;
46e89793 1509 fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1510
1511 xminTrack[0] =0.; xmaxTrack[0]=250;
46e89793 1512 fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1513
1514
1515 for (Int_t ivar=0;ivar<6;ivar++){
1516 fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1517 fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1518 fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1519 fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1520 fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1521 fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1522 fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
1523 fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1524 fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1525 fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1526 fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1527 fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1528 }
1529
1530
af6a50bb 1531 BinLogX(fHisDeltaY,3);
1532 BinLogX(fHisDeltaZ,3);
1533 BinLogX(fHisRMSY,3);
1534 BinLogX(fHisRMSZ,3);
1535 BinLogX(fHisQmax,3);
1536 BinLogX(fHisQtot,3);
1537
1538}
1539
1540void AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
1541 //
1542 // Add histograms
1543 //
1544 if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
1545 if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
1546 if (calib->fHisRMSY) fHisRMSY->Add(calib->fHisRMSY);
1547 if (calib->fHisRMSZ) fHisRMSZ->Add(calib->fHisRMSZ);
1548}
46e89793 1549
1550
1551
1552void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
1553 //
1554 // Dump summary info
1555 //
1556 // 0.OBJ: TAxis var
1557 // 1.OBJ: TAxis pad type
1558 // 2.OBJ: TAxis z
1559 // 3.OBJ: TAxis Qmax
1560 // 4.OBJ: TAxis cog
1561 // 5.OBJ: TAxis tan(angle)
1562 //
1563 if (ptype>3) return;
1564 Int_t idim[6]={0,1,2,3,4,5};
1565 TString hname[4]={"dy","dz","rmsy","rmsz"};
1566 //
1567 Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
1568 Int_t first5=hisInput->GetAxis(5)->GetFirst();
1569 Int_t last5 =hisInput->GetAxis(5)->GetLast();
1570 //
1571 for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){ // axis 5 - local angle
1572 Bool_t bin5All=kTRUE;
1573 if (ibin5>=first5){
1574 hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
1575 if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
1576 bin5All=kFALSE;
1577 }
1578 Double_t x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
1579 THnSparse * his5= hisInput->Projection(5,idim); //projected histogram according selection 5
1580 Int_t nbins4=his5->GetAxis(4)->GetNbins();
1581 Int_t first4=his5->GetAxis(4)->GetFirst();
1582 Int_t last4 =his5->GetAxis(4)->GetLast();
1583
1584 for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){ // axis 4 - cog
1585 Bool_t bin4All=kTRUE;
1586 if (ibin4>=first4){
1587 his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
1588 if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
1589 bin4All=kFALSE;
1590 }
1591 Double_t x4= his5->GetAxis(4)->GetBinCenter(ibin4);
1592 THnSparse * his4= his5->Projection(4,idim); //projected histogram according selection 4
1593 //
1594 Int_t nbins3=his4->GetAxis(3)->GetNbins();
1595 Int_t first3=his4->GetAxis(3)->GetFirst();
1596 Int_t last3 =his4->GetAxis(3)->GetLast();
1597 //
1598 for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - Qmax
1599 Bool_t bin3All=kTRUE;
1600 if (ibin3>=first3){
1601 his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
1602 bin3All=kFALSE;
1603 }
1604 Double_t x3= his4->GetAxis(3)->GetBinCenter(ibin3);
1605 THnSparse * his3= his4->Projection(3,idim); //projected histogram according selection 3
1606 //
1607 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1608 Int_t first2 = his3->GetAxis(2)->GetFirst();
1609 Int_t last2 = his3->GetAxis(2)->GetLast();
1610 //
1611 for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){ // axis 2 - z
1612 Bool_t bin2All=kTRUE;
1613 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1614 if (ibin2>=first2){
1615 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
1616 if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
1617 bin2All=kFALSE;
1618 }
1619 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
1620 //
1621 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
1622 Int_t first1 = his2->GetAxis(1)->GetFirst();
1623 Int_t last1 = his2->GetAxis(1)->GetLast();
1624 for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){ //axis 1 - pad type
1625 Bool_t bin1All=kTRUE;
1626 if (ibin1>=first1){
1627 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1628 bin1All=kFALSE;
1629 }
1630 Double_t x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
1631 TH1 * hisDelta = his2->Projection(0);
1632 Double_t entries = hisDelta->GetEntries();
1633 Double_t mean=0, rms=0;
1634 if (entries>10){
1635 mean = hisDelta->GetMean();
1636 rms = hisDelta->GetRMS();
1637 hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
1638 mean = hisDelta->GetMean();
1639 rms = hisDelta->GetRMS();
1640 }
1641 //
1642 (*pcstream)<<hname[ptype].Data()<<
1643 // flag - intgrated values for given bin
60e31b96 1644 "angleA="<<bin5All<<
46e89793 1645 "cogA="<<bin4All<<
1646 "qmaxA="<<bin3All<<
1647 "zA="<<bin2All<<
1648 "ipadA="<<bin1All<<
1649 // center of bin value
60e31b96 1650 "angle="<<x5<< // local angle
1651 "cog="<<x4<< // distance cluster to center
1652 "qmax="<<x3<< // qmax
1653 "z="<<x2<< // z of the cluster
1654 "ipad="<<x1<< // type of the region
46e89793 1655 // mean values
1656 //
1657 "entries="<<entries<<
1658 "mean="<<mean<<
1659 "rms="<<rms<<
1660 "ptype="<<ptype<< //
1661 "\n";
1662 delete hisDelta;
1663 printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);
1664 }//loop z
1665 delete his2;
1666 } // loop Q max
1667 delete his3;
1668 } // loop COG
1669 delete his4;
1670 }//loop local angle
1671 delete his5;
1672 }
1673}
7d14c1c1 1674
1675
1676int AliTPCcalibTracks::GetTHnStat(const THnBase *H, THnBase *&Mean, THnBase *&Sigma, THnBase *&Entr )
1677{
1678 // H is THnSparse:
1679 //
1680 // OBJ: TAxis var var
1681 // OBJ: TAxis axis 1
1682 // OBJ: TAxis axis 2
1683 // ...
1684
1685 // Output:
1686
1687 // Mean, Sigma and Entr is THnF of mean, RMS and entries of var:
1688 //
1689 // OBJ: TAxis axis 1
1690 // OBJ: TAxis axis 2
1691 // ..
1692
1693 Mean = 0;
1694 Sigma = 0;
1695 Entr = 0;
1696 if( !H ) return -1;
1697
1698 Bool_t ok = kTRUE;
1699
1700 int nDim = H->GetNdimensions();
1701 Long64_t nFilledBins = H->GetNbins();
1702 Long64_t nStatBins = 0;
1703
1704 Int_t *nBins = new Int_t [nDim];
1705 Double_t *xMin = new Double_t [nDim];
1706 Double_t *xMax = new Double_t [nDim];
1707 Int_t *idx = new Int_t [nDim];
1708
1709 TString nameMean = H->GetName();
1710 TString nameSigma = H->GetName();
1711 TString nameEntr = H->GetName();
1712 nameMean+="_Mean";
1713 nameSigma+="_Sigma";
1714 nameEntr+="_Entr";
1715
1716 ok = ok && ( nBins && xMin && xMax && idx );
1717
1718 if( ok ){
1719 for( int i=0; i<nDim; i++){
1720 xMin[i] = H->GetAxis(i)->GetXmin();
1721 xMax[i] = H->GetAxis(i)->GetXmax();
1722 nBins[i] = H->GetAxis(i)->GetNbins();
1723 }
1724
1725 Mean = new THnSparseF( nameMean.Data(), nameMean.Data(), nDim-1, nBins+1, xMin+1, xMax+1 );
1726 Sigma = new THnSparseF( nameSigma.Data(), nameSigma.Data(), nDim-1, nBins+1, xMin+1, xMax+1 );
1727 Entr = new THnSparseF( nameEntr.Data(), nameEntr.Data(), nDim-1, nBins+1, xMin+1, xMax+1 );
1728 }
1729
1730 ok = ok && ( Mean && Sigma && Entr );
1731
1732 if( ok ){
1733 for( int i=0; i<nDim-1; i++){
1734 const TAxis *ax = H->GetAxis(i+1);
1735 Mean->GetAxis(i)->SetName(ax->GetName());
1736 Mean->GetAxis(i)->SetTitle(ax->GetTitle());
1737 Sigma->GetAxis(i)->SetName(ax->GetName());
1738 Sigma->GetAxis(i)->SetTitle(ax->GetTitle());
1739 Entr->GetAxis(i)->SetName(ax->GetName());
1740 Entr->GetAxis(i)->SetTitle(ax->GetTitle());
1741 if( ax->GetXbins()->fN>0 ){
1742 Mean->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1743 Sigma->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1744 Entr->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1745 }
1746 }
1747
1748 // Allocate bins
1749
1750 for( Long64_t binS=0; binS<nFilledBins; binS++){
1751 H->GetBinContent(binS,idx);
1752 Mean->GetBin(idx+1,kTRUE);
1753 Sigma->GetBin(idx+1,kTRUE);
1754 Entr->GetBin(idx+1,kTRUE);
1755 }
1756
1757 nStatBins = Mean->GetNbins();
1758
1759 }
1760
1761
1762 Long64_t *hMap = new Long64_t[nFilledBins];
1763 Double_t *hVal = new Double_t[nFilledBins];
1764 Double_t *hEntr = new Double_t[nFilledBins];
1765 Double_t *meanD = new Double_t[nStatBins];
1766 Double_t *sigmaD = new Double_t[nStatBins];
1767
1768 ok = ok && ( hMap && hVal && hEntr && meanD && sigmaD );
1769
1770 // Map bins to Stat bins
1771
1772 if( ok ){
1773 for( Long64_t binS=0; binS<nFilledBins; binS++){
1774 double val = H->GetBinContent(binS,idx);
1775 Long64_t bin = Mean->GetBin(idx+1,kFALSE);
1776 ok = ok && ( bin>=0 && bin<nStatBins && bin==Sigma->GetBin(idx+1,kFALSE) && bin== Entr->GetBin(idx+1,kFALSE) );
1777 if( !ok ) break;
1778 if( val < 0 ){
1779 cout << "AliTPCcalibTracks: GetSparseStat : Unexpected content of the input THn"<<endl;
1780 bin = -1;
1781 }
1782 if( idx[0]<1 || idx[0]>nBins[0] ) bin = -1;
1783 hMap[binS] = bin;
1784 hVal[binS] = idx[0];
1785 hEntr[binS] = val;
1786 }
1787 }
1788
1789 // do 2 iteration with cut at 4 Sigma
1790
1791 for( int iter=0; ok && iter<2; iter++ ){
1792
1793 // clean up entries
1794
1795 for( Long64_t bin=0; bin < nStatBins; bin++){
1796 Entr->SetBinContent(bin, 0);
1797 meanD[bin]=0;
1798 sigmaD[bin]=0;
1799 }
1800
1801 // get Entries and Mean
1802
1803 for( Long64_t binS=0; binS<nFilledBins; binS++){
1804 Long64_t bin = hMap[binS];
1805 Double_t val = hVal[binS];
1806 Double_t entr = hEntr[binS];
1807 if( bin < 0 ) continue;
1808 if( iter!=0 ){ // cut
1809 double s2 = Sigma->GetBinContent(bin);
1810 double d = val - Mean->GetBinContent(bin);
1811 if( d*d>s2*16 ) continue;
1812 }
1813 meanD[bin]+= entr*val;
1814 Entr->AddBinContent(bin,entr);
1815 }
1816
1817 for( Long64_t bin=0; bin<nStatBins; bin++){
1818 double val = meanD[bin];
1819 double sum = Entr->GetBinContent(bin);
1820 if( sum>=1 ){
1821 Mean->SetBinContent(bin, val/sum );
1822 }
1823 else Mean->SetBinContent(bin, 0);
1824 Entr->SetBinContent(bin, 0);
1825 }
1826
1827 // get RMS
1828
1829 for( Long64_t binS=0; binS<nFilledBins; binS++){
1830 Long64_t bin = hMap[binS];
1831 Double_t val = hVal[binS];
1832 Double_t entr = hEntr[binS];
1833 if( bin < 0 ) continue;
1834 double d2 = val - Mean->GetBinContent(bin);
1835 d2 *= d2;
1836 if( iter!=0 ){ // cut
1837 double s2 = Sigma->GetBinContent(bin);
1838 if( d2>s2*16 ) continue;
1839 }
1840 sigmaD[bin] += entr*d2;
1841 Entr->AddBinContent(bin,entr);
1842 }
1843
1844 for( Long64_t bin=0; bin<nStatBins; bin++ ){
1845 double val = sigmaD[bin];
1846 double sum = Entr->GetBinContent(bin);
1847 if( sum>=1 && val>=0 ){
1848 Sigma->SetBinContent(bin, val/sum );
1849 }
1850 else Sigma->SetBinContent(bin, 0);
1851 }
1852 }
1853
1854 // scale the Mean and the Sigma
1855
1856 if( ok ){
1857 double v0 = H->GetAxis(0)->GetBinCenter(1);
1858 double vscale = H->GetAxis(0)->GetBinWidth(1);
1859
1860 for( Long64_t bin=0; bin<nStatBins; bin++){
1861 double m = Mean->GetBinContent(bin);
1862 double s2 = Sigma->GetBinContent(bin);
1863 double sum = Entr->GetBinContent(bin);
1864 if( sum>0 && s2>=0 ){
1865 Mean->SetBinContent(bin, v0 + (m-1)*vscale );
1866 Sigma->SetBinContent(bin, sqrt(s2)*vscale );
1867 }
1868 else{
1869 Mean->SetBinContent(bin, 0);
1870 Sigma->SetBinContent(bin, 0);
1871 Entr->SetBinContent(bin, 0);
1872 }
1873 }
1874 }
1875
d4fa6a95 1876 delete[] nBins;
1877 delete[] xMin;
1878 delete[] xMax;
1879 delete[] idx;
1880 delete[] hMap;
1881 delete[] hVal;
1882 delete[] hEntr;
1883 delete[] meanD;
1884 delete[] sigmaD;
7d14c1c1 1885
1886 if( !ok ){
1887 cout << "AliTPCcalibTracks: GetSparseStat : No memory or internal Error "<<endl;
1888 delete Mean;
1889 delete Sigma;
1890 delete Entr;
1891 Mean =0;
1892 Sigma = 0;
1893 Entr = 0;
1894 return -1;
1895 }
1896
1897 return 0;
1898}
1899
1900
1901
1902int AliTPCcalibTracks::CreateWaveCorrection(const THnBase *DeltaY, THnBase *&MeanY, THnBase *&SigmaY, THnBase *&EntrY,
1903 Bool_t MirrorZ, Bool_t MirrorPad, Bool_t MirrorAngle, Int_t MinStat )
1904{
1905 // DeltaY is THnSparse:
1906 //
1907 // OBJ: TAxis 0 var var
1908 // OBJ: TAxis 1 pad type pad type
1909 // OBJ: TAxis 2 z z
1910 // OBJ: TAxis 3 Qmax Qmax
1911 // OBJ: TAxis 4 cog cog
1912 // OBJ: TAxis 5 tan(angle) tan(angle)
1913 // ...
1914
1915 MeanY = 0;
1916 SigmaY = 0;
1917 EntrY = 0;
1918
1919 int nDim = DeltaY->GetNdimensions();
d4fa6a95 1920 if( nDim<6 ){
1921 cout << "AliTPCcalibTracks: CreateWaveCorrection: Unknown input"<<endl;
1922 return -1;
7d14c1c1 1923 }
1924
d4fa6a95 1925 Int_t *nBins = new Int_t[nDim];
1926 Int_t *nBinsNew = new Int_t[nDim];
1927 Double_t *xMin = new Double_t[nDim];
1928 Double_t *xMax = new Double_t[nDim];
1929 Int_t *idx = new Int_t[nDim];
1930 THnSparseD *mergedDeltaY = 0;
7d14c1c1 1931
d4fa6a95 1932 int ret = 0;
7d14c1c1 1933
d4fa6a95 1934 if( !nBins || !nBinsNew || !xMin || !xMax || !idx ){
1935 ret = -1;
1936 cout << "AliTPCcalibTracks: CreateWaveCorrection: Out of memory"<<endl;
7d14c1c1 1937 }
d4fa6a95 1938
1939 if( ret==0 ){
1940
1941 for( int i=0; i<nDim; i++){
1942 xMin[i] = DeltaY->GetAxis(i)->GetXmin();
1943 xMax[i] = DeltaY->GetAxis(i)->GetXmax();
1944 nBins[i] = DeltaY->GetAxis(i)->GetNbins();
1945 nBinsNew[i] = nBins[i];
1946 }
7d14c1c1 1947
d4fa6a95 1948 // Merge cog axis
1949 if( MirrorPad ){
1950 Int_t centralBin = DeltaY->GetAxis(4)->FindFixBin(0.5);
1951 xMin[4] = DeltaY->GetAxis(4)->GetBinLowEdge(centralBin);
1952 nBinsNew[4] = nBinsNew[4]-centralBin+1;
7d14c1c1 1953 }
7d14c1c1 1954
d4fa6a95 1955 // Merge Z axis
1956 if( MirrorZ ){
1957 Int_t centralBin = DeltaY->GetAxis(2)->FindFixBin(0.0);
1958 xMin[2] = DeltaY->GetAxis(2)->GetBinLowEdge(centralBin)-0.0;
1959 nBinsNew[2] = nBinsNew[2]-centralBin+1;
1960 }
7d14c1c1 1961
d4fa6a95 1962 // Merge Angle axis
1963 if( MirrorAngle ){
1964 Int_t centralBin = DeltaY->GetAxis(5)->FindFixBin(0.0);
1965 xMin[5] = DeltaY->GetAxis(5)->GetBinLowEdge(centralBin)-0.0;
1966 nBinsNew[5] = nBinsNew[5]-centralBin+1;
7d14c1c1 1967 }
1968
d4fa6a95 1969 // Merge Sparse
1970
1971 mergedDeltaY = new THnSparseD("mergedDeltaY", "mergedDeltaY",nDim, nBinsNew, xMin, xMax );
1972
1973 if( !mergedDeltaY ){
1974 cout << "AliTPCcalibTracks: CreateWaveCorrection: Can not copy a Sparse"<<endl;
1975 ret = -1;
7d14c1c1 1976 }
d4fa6a95 1977 }
1978
1979 if( ret == 0 ){
1980
1981 for( int i=0; i<nDim; i++){
1982 const TAxis *ax = DeltaY->GetAxis(i);
1983 mergedDeltaY->GetAxis(i)->SetName(ax->GetName());
1984 mergedDeltaY->GetAxis(i)->SetTitle(ax->GetTitle());
1985 if( ax->GetXbins()->fN>0 ){
1986 mergedDeltaY->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1987 }
7d14c1c1 1988 }
d4fa6a95 1989
1990 for( Long64_t binS=0; binS<DeltaY->GetNbins(); binS++){
1991 Double_t stat = DeltaY->GetBinContent(binS,idx);
1992 if( stat < 1 ) continue;
1993 bool swap0=0;
1994
1995 if( MirrorPad && idx[4]>0){ // underflow reserved for contains one-pad clusters, don't mirror
1996 Double_t v = DeltaY->GetAxis(4)->GetBinCenter(idx[4]);
1997 if( v < 0.5 ) swap0 = !swap0;
1998 idx[4] = mergedDeltaY->GetAxis(4)->FindFixBin( 0.5 + TMath::Abs(0.5 - v) );
1999 }
2000
2001 if( MirrorZ ){
2002 Double_t v = DeltaY->GetAxis(2)->GetBinCenter(idx[2]);
2003 if( v < 0.0 ) swap0 = !swap0;
2004 if( idx[2]<=0 ) idx[2] = nBinsNew[2]+1;
2005 else idx[2] = mergedDeltaY->GetAxis(2)->FindFixBin( TMath::Abs(v) );
2006 }
2007
2008 if( MirrorAngle ){
2009 Double_t v = DeltaY->GetAxis(5)->GetBinCenter(idx[5]);
2010 if( idx[5]<=0 ) idx[5] = nBinsNew[5]+1;
2011 else idx[5] = mergedDeltaY->GetAxis(5)->FindFixBin( TMath::Abs(v) );
2012 }
7d14c1c1 2013
d4fa6a95 2014 if( swap0 ){
2015 if( idx[0]<=0 ) idx[0] = nBinsNew[0]+1;
2016 else if( idx[0] >= nBins[0]+1 ) idx[0] = 0;
2017 else {
2018 Double_t v = DeltaY->GetAxis(0)->GetBinCenter(idx[0]);
2019 idx[0] = mergedDeltaY->GetAxis(0)->FindFixBin(-v);
2020 }
7d14c1c1 2021 }
d4fa6a95 2022
2023 Long64_t bin = mergedDeltaY->GetBin(idx,kTRUE);
2024 if( bin<0 ){
2025 cout << "AliTPCcalibTracks: CreateWaveCorrection : wrong bining"<<endl;
2026 continue;
2027 }
2028 mergedDeltaY->AddBinContent(bin,stat);
7d14c1c1 2029 }
2030
d4fa6a95 2031 ret = GetTHnStat(mergedDeltaY, MeanY, SigmaY, EntrY );
7d14c1c1 2032 }
2033
d4fa6a95 2034 if( ret==0 ){
2035
2036 MeanY->SetName("TPCWaveCorrectionMap");
2037 MeanY->SetTitle("TPCWaveCorrectionMap");
2038 SigmaY->SetName("TPCResolutionMap");
2039 SigmaY->SetTitle("TPCResolutionMap");
2040 EntrY->SetName("TPCWaveCorrectionEntr");
2041 EntrY->SetTitle("TPCWaveCorrectionEntr");
7d14c1c1 2042
d4fa6a95 2043 for( Long64_t bin=0; bin<MeanY->GetNbins(); bin++){
2044 Double_t stat = EntrY->GetBinContent(bin,idx);
2045
2046 // Normalize : Set no correction for one pad clusters
2047
2048 if( idx[3]<=0 ) MeanY->SetBinContent(bin,0);
2049
2050 // Suppress bins with low statistic
2051
2052 if( stat<MinStat ){
2053 EntrY->SetBinContent(bin,0);
2054 MeanY->SetBinContent(bin,0);
2055 SigmaY->SetBinContent(bin,-1);
2056 }
2057 }
2058
2059 }
7d14c1c1 2060
d4fa6a95 2061 delete[] nBins;
2062 delete[] nBinsNew;
2063 delete[] xMin;
2064 delete[] xMax;
2065 delete[] idx;
2066 delete mergedDeltaY;
7d14c1c1 2067
d4fa6a95 2068 if( ret!=0 ){
2069 delete MeanY;
2070 delete SigmaY;
2071 delete EntrY;
2072 MeanY = 0;
2073 SigmaY = 0;
2074 EntrY = 0;
7d14c1c1 2075 }
7d14c1c1 2076
d4fa6a95 2077 return ret;
7d14c1c1 2078}
2079
2080int AliTPCcalibTracks::UpdateClusterParam( AliTPCClusterParam *cParam, Bool_t MirrorZ, Bool_t MirrorPad, Bool_t MirrorAngle, Int_t MinStat )
2081{
2082 if( !cParam || !fHisDeltaY) return -1;
2083
2084 THnBase *meanY = 0;
2085 THnBase *sigmaY = 0;
2086 THnBase *entrY = 0;
2087 int ret = CreateWaveCorrection(fHisDeltaY, meanY, sigmaY, entrY, MirrorZ, MirrorAngle, MirrorPad, MinStat );
2088 if( ret<0 ) return ret;
2089 cParam->SetWaveCorrectionMap(meanY);
2090 cParam->SetResolutionYMap(sigmaY);
2091 delete meanY;
2092 delete sigmaY;
2093 delete entrY;
2094 return 0;
2095}