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