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