]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix for coverity issues : 17923
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTracks.cxx
CommitLineData
10757ee9 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17///////////////////////////////////////////////////////////////////////////////
18// //
19// Class to analyse tracks for calibration //
20// to be used as a component in AliTPCSelectorTracks //
21// In the constructor you have to specify name and title //
22// to get the Object out of a file. //
23// The parameter 'clusterParam', a AliTPCClusterParam object //
24// (needed for TPC cluster error and shape parameterization) //
25// Normally you get this object out of the file 'TPCClusterParam.root' //
26// In the parameter 'cuts' the cuts are specified, that decide //
27// weather a track will be accepted for calibration or not. //
28// //
29//
30// The data flow:
31//
32/*
33 Raw Data -> Local Reconstruction -> Tracking -> Calibration -> RefData (component itself)
34 Offline/HLT Offline/HLT OCDB entries (AliTPCClusterParam)
35*/
36
9f0beaf7 37
10757ee9 38 //
39///////////////////////////////////////////////////////////////////////////////
40
41//
42// ROOT includes
43//
44#include <iostream>
45#include <fstream>
46using namespace std;
47
48#include <TROOT.h>
49#include <TChain.h>
50#include <TFile.h>
51#include <TH3F.h>
5b00528f 52#include <TProfile.h>
53
10757ee9 54//
55//#include <TPDGCode.h>
56#include <TStyle.h>
57#include "TLinearFitter.h"
58//#include "TMatrixD.h"
59#include "TTreeStream.h"
60#include "TF1.h"
61#include <TCanvas.h>
62#include <TGraph2DErrors.h>
63#include "TPostScript.h"
64#include "TCint.h"
65
66#include <TH2D.h>
67#include <TF2.h>
68#include <TSystem.h>
69#include <TCollection.h>
70#include <iostream>
71#include <TLinearFitter.h>
2e5bcb67 72#include <TString.h>
10757ee9 73
74//
75// AliROOT includes
76//
77#include "AliMagF.h"
78#include "AliTracker.h"
79#include "AliESD.h"
80#include "AliESDtrack.h"
81#include "AliESDfriend.h"
82#include "AliESDfriendTrack.h"
83#include "AliTPCseed.h"
84#include "AliTPCclusterMI.h"
85#include "AliTPCROC.h"
86
87#include "AliTPCParamSR.h"
88#include "AliTrackPointArray.h"
89#include "AliTPCcalibTracks.h"
90#include "AliTPCClusterParam.h"
91#include "AliTPCcalibTracksCuts.h"
10757ee9 92#include "AliTPCCalPad.h"
93#include "AliTPCCalROC.h"
94#include "TText.h"
95#include "TPaveText.h"
96#include "TSystem.h"
2e5bcb67 97#include "TStatToolkit.h"
98#include "TCut.h"
af6a50bb 99#include "THnSparse.h"
46e89793 100#include "AliRieman.h"
10757ee9 101
10757ee9 102
103
104ClassImp(AliTPCcalibTracks)
105
106
107AliTPCcalibTracks::AliTPCcalibTracks():
c32da879 108 AliTPCcalibBase(),
2c632057 109 fClusterParam(0),
2c632057 110 fROC(0),
af6a50bb 111 fHisDeltaY(0), // THnSparse - delta Y
112 fHisDeltaZ(0), // THnSparse - delta Z
113 fHisRMSY(0), // THnSparse - rms Y
114 fHisRMSZ(0), // THnSparse - rms Z
115 fHisQmax(0), // THnSparse - qmax
116 fHisQtot(0), // THnSparse - qtot
2c632057 117 fArrayQDY(0),
118 fArrayQDZ(0),
119 fArrayQRMSY(0),
120 fArrayQRMSZ(0),
2c632057 121 fResolY(0),
122 fResolZ(0),
123 fRMSY(0),
124 fRMSZ(0),
125 fCuts(0),
2c632057 126 fRejectedTracksHisto(0),
2c632057 127 fClusterCutHisto(0),
128 fCalPadClusterPerPad(0),
9f0beaf7 129 fCalPadClusterPerPadRaw(0)
10757ee9 130{
131 //
132 // AliTPCcalibTracks default constructor
133 //
ae28e92e 134 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;
10757ee9 135}
136
137
b42cf9ed 138
139AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
ae28e92e 140 AliTPCcalibBase(calibTracks),
2c632057 141 fClusterParam(0),
2c632057 142 fROC(0),
af6a50bb 143 fHisDeltaY(0), // THnSparse - delta Y
144 fHisDeltaZ(0), // THnSparse - delta Z
145 fHisRMSY(0), // THnSparse - rms Y
146 fHisRMSZ(0), // THnSparse - rms Z
147 fHisQmax(0), // THnSparse - qmax
148 fHisQtot(0), // THnSparse - qtot
2c632057 149 fArrayQDY(0),
150 fArrayQDZ(0),
151 fArrayQRMSY(0),
152 fArrayQRMSZ(0),
2c632057 153 fResolY(0),
154 fResolZ(0),
155 fRMSY(0),
156 fRMSZ(0),
157 fCuts(0),
2c632057 158 fRejectedTracksHisto(0),
2c632057 159 fClusterCutHisto(0),
160 fCalPadClusterPerPad(0),
9f0beaf7 161 fCalPadClusterPerPadRaw(0)
2c632057 162{
10757ee9 163 //
164 // AliTPCcalibTracks copy constructor
165 //
ae28e92e 166 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
10757ee9 167
168 Bool_t dirStatus = TH1::AddDirectoryStatus();
169 TH1::AddDirectory(kFALSE);
170
171 Int_t length = -1;
10757ee9 172
b42cf9ed 173 (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
10757ee9 174 fArrayQDY= new TObjArray(length);
175 fArrayQDZ= new TObjArray(length);
176 fArrayQRMSY= new TObjArray(length);
177 fArrayQRMSZ= new TObjArray(length);
178 for (Int_t i = 0; i < length; i++) {
b42cf9ed 179 fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
180 fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
181 fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
182 fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
10757ee9 183 }
184
b42cf9ed 185 (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
10757ee9 186 fResolY = new TObjArray(length);
187 fResolZ = new TObjArray(length);
188 fRMSY = new TObjArray(length);
189 fRMSZ = new TObjArray(length);
190 for (Int_t i = 0; i < length; i++) {
b42cf9ed 191 fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
192 fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
193 fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
194 fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
10757ee9 195 }
196
10757ee9 197
b42cf9ed 198 fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
199 fRejectedTracksHisto = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
b42cf9ed 200 fCalPadClusterPerPad = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
201 fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
202
203 fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(),
204 calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
b42cf9ed 205 SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
10757ee9 206 TH1::AddDirectory(dirStatus); // set status back to original status
207// cout << "+++++ end of copy constructor +++++" << endl; // TO BE REMOVED
208}
209
210
b42cf9ed 211AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
212 //
213 // assgnment operator
214 //
215 if (this != &calibTracks) {
216 new (this) AliTPCcalibTracks(calibTracks);
217 }
218 return *this;
219
220}
221
222
10757ee9 223AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts, Int_t logLevel) :
c32da879 224 AliTPCcalibBase(),
2c632057 225 fClusterParam(0),
2c632057 226 fROC(0),
af6a50bb 227 fHisDeltaY(0), // THnSparse - delta Y
228 fHisDeltaZ(0), // THnSparse - delta Z
229 fHisRMSY(0), // THnSparse - rms Y
230 fHisRMSZ(0), // THnSparse - rms Z
231 fHisQmax(0), // THnSparse - qmax
232 fHisQtot(0), // THnSparse - qtot
2c632057 233 fArrayQDY(0),
234 fArrayQDZ(0),
235 fArrayQRMSY(0),
236 fArrayQRMSZ(0),
2c632057 237 fResolY(0),
238 fResolZ(0),
239 fRMSY(0),
240 fRMSZ(0),
241 fCuts(0),
2c632057 242 fRejectedTracksHisto(0),
2c632057 243 fClusterCutHisto(0),
244 fCalPadClusterPerPad(0),
9f0beaf7 245 fCalPadClusterPerPadRaw(0)
2c632057 246 {
10757ee9 247 //
248 // AliTPCcalibTracks constructor
249 // specify 'name' and 'title' of your object
250 // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
251 // In the parameter 'cuts' the cuts are specified, that decide
252 // weather a track will be accepted for calibration or not.
ae28e92e 253 //
254 // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
10757ee9 255 //
256 // All histograms are instatiated in this constructor.
257 //
c32da879 258 this->SetName(name);
259 this->SetTitle(title);
260
ae28e92e 261 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
10757ee9 262 G__SetCatchException(0);
263
264 fClusterParam = clusterParam;
265 if (fClusterParam){
266 fClusterParam->SetInstance(fClusterParam);
267 }
268 else {
269 Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
270 }
271 fCuts = cuts;
ae28e92e 272 SetDebugLevel(logLevel);
af6a50bb 273 MakeHistos();
10757ee9 274
275 TH1::AddDirectory(kFALSE);
276
8b9b3187 277 fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
10757ee9 278 fCalPadClusterPerPad = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
279 fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
280 fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
281
10757ee9 282
283 TH1::AddDirectory(kFALSE);
284
10757ee9 285
286 fResolY = new TObjArray(3);
287 fResolZ = new TObjArray(3);
288 fRMSY = new TObjArray(3);
289 fRMSZ = new TObjArray(3);
290 TH3F * his3D;
291 //
292 his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
293 fResolY->AddAt(his3D,0);
294 his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
295 fResolY->AddAt(his3D,1);
296 his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
297 fResolY->AddAt(his3D,2);
298 //
299 his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
300 fResolZ->AddAt(his3D,0);
301 his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
302 fResolZ->AddAt(his3D,1);
303 his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
304 fResolZ->AddAt(his3D,2);
305 //
306 his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
307 fRMSY->AddAt(his3D,0);
308 his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
309 fRMSY->AddAt(his3D,1);
310 his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
311 fRMSY->AddAt(his3D,2);
312 //
313 his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
314 fRMSZ->AddAt(his3D,0);
315 his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
316 fRMSZ->AddAt(his3D,1);
317 his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
318 fRMSZ->AddAt(his3D,2);
319 //
320
321 TH1::AddDirectory(kFALSE);
322
323 fArrayQDY = new TObjArray(300);
324 fArrayQDZ = new TObjArray(300);
325 fArrayQRMSY = new TObjArray(300);
326 fArrayQRMSZ = new TObjArray(300);
327 for (Int_t iq = 0; iq <= 10; iq++){
328 for (Int_t ipad = 0; ipad < 3; ipad++){
329 Int_t bin = GetBin(iq, ipad);
330 Float_t qmean = GetQ(bin);
d23650ed 331 char hname[200];
f69a0aa7 332 snprintf(hname,100,"ResolY Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 333 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
10757ee9 334 fArrayQDY->AddAt(his3D, bin);
f69a0aa7 335 snprintf(hname,100,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 336 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
10757ee9 337 fArrayQDZ->AddAt(his3D, bin);
f69a0aa7 338 snprintf(hname,100,"RMSY Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 339 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
10757ee9 340 fArrayQRMSY->AddAt(his3D, bin);
f69a0aa7 341 snprintf(hname,100,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
a4c5fca9 342 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
10757ee9 343 fArrayQRMSZ->AddAt(his3D, bin);
344 }
345 }
346
10757ee9 347
10757ee9 348
ae28e92e 349 if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl;
10757ee9 350 cout << "end of main constructor" << endl; // TO BE REMOVED
351}
352
353
354AliTPCcalibTracks::~AliTPCcalibTracks() {
355 //
356 // AliTPCcalibTracks destructor
357 //
358
ae28e92e 359 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
10757ee9 360 Int_t length = 0;
10757ee9 361
10757ee9 362
f69a0aa7 363 if (fResolY) {
364 length = fResolY->GetEntriesFast();
365 for (Int_t i = 0; i < length; i++){
366 delete fResolY->At(i);
367 delete fResolZ->At(i);
368 delete fRMSY->At(i);
369 delete fRMSZ->At(i);
370 }
371 delete fResolY;
372 delete fResolZ;
373 delete fRMSY;
374 delete fRMSZ;
10757ee9 375 }
10757ee9 376
f69a0aa7 377 if (fArrayQDY) {
378 length = fArrayQDY->GetEntriesFast();
379 for (Int_t i = 0; i < length; i++){
380 delete fArrayQDY->At(i);
381 delete fArrayQDZ->At(i);
382 delete fArrayQRMSY->At(i);
383 delete fArrayQRMSZ->At(i);
384 }
10757ee9 385 }
386
10757ee9 387
9f0beaf7 388
10757ee9 389 delete fArrayQDY;
390 delete fArrayQDZ;
391 delete fArrayQRMSY;
392 delete fArrayQRMSZ;
10757ee9 393
10757ee9 394 delete fRejectedTracksHisto;
395 delete fClusterCutHisto;
10757ee9 396 if (fCalPadClusterPerPad) delete fCalPadClusterPerPad;
397 if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
af6a50bb 398 delete fHisDeltaY; // THnSparse - delta Y
399 delete fHisDeltaZ; // THnSparse - delta Z
400 delete fHisRMSY; // THnSparse - rms Y
401 delete fHisRMSZ; // THnSparse - rms Z
402 delete fHisQmax; // THnSparse - qmax
403 delete fHisQtot; // THnSparse - qtot
404
10757ee9 405}
406
407
10757ee9 408
c32da879 409void AliTPCcalibTracks::Process(AliTPCseed *track){
10757ee9 410 //
411 // To be called in the selector
412 // first AcceptTrack is evaluated, then calls all the following analyse functions:
413 // FillResolutionHistoLocal(track)
af6a50bb 414
10757ee9 415 //
ae28e92e 416 if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
10757ee9 417 Int_t accpetStatus = AcceptTrack(track);
418 if (accpetStatus == 0) {
419 FillResolutionHistoLocal(track);
10757ee9 420 }
421 else fRejectedTracksHisto->Fill(accpetStatus);
422}
423
424
425
426Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
427 //
428 // calculate bins for given q and pad type
429 // used in TObjArray
430 //
431 Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );
432 res *= 3;
433 res += pad;
434 return res;
435}
436
437
438Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
439 //
440 // calculate bins for given iq and pad type
441 // used in TObjArray
442 //
443 return iq * 3 + pad;;
444}
445
446
447Float_t AliTPCcalibTracks::GetQ(Int_t bin){
448 //
449 // returns to bin belonging charge
450 // (bin / 3 + 3)^2
451 //
452 Int_t bin0 = bin / 3;
453 bin0 += 3;
454 return bin0 * bin0;
455}
456
457
458Float_t AliTPCcalibTracks::GetPad(Int_t bin){
459 //
460 // returns to bin belonging pad
461 // bin % 3
462 //
463 return bin % 3;
464}
465
466
467
468Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
469 //
470 // Function, that decides wheather a given track is accepted for
471 // the analysis or not.
472 // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
473 // Returns 0 if a track is accepted or an integer different from 0
474 // to indicate the failed cut
475 //
476 const Int_t kMinClusters = fCuts->GetMinClusters();
477 const Float_t kMinRatio = fCuts->GetMinRatio();
478 const Float_t kMax1pt = fCuts->GetMax1pt();
479 const Float_t kEdgeYXCutNoise = fCuts->GetEdgeYXCutNoise();
480 const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
481
482 //
483 // edge induced noise tracks - NEXT RELEASE will be removed during tracking
484 if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
485 if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
486 if (track->GetNumberOfClusters() < kMinClusters) return 2;
487 Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
488 if (ratio < kMinRatio) return 3;
489// Float_t mpt = track->Get1Pt(); // Get1Pt() doesn't exist any more
490 Float_t mpt = track->GetSigned1Pt();
491 if (TMath::Abs(mpt) > kMax1pt) return 4;
492 //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
493 //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
494 //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
495
430a3403 496 if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");
10757ee9 497 return 0;
498}
499
500
501void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
502 //
503 // fill resolution histograms - localy - tracklet in the neighborhood
504 // write debug information to 'TPCSelectorDebug.root'
505 //
506 // _ the main function, called during track analysis _
507 //
508 // loop over all padrows along the track
509 // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
510 //
511 // loop again over all padrows along the track
512 // fit tracklet (clusters in given padrow +- kDelta padrows)
513 // with polynom of 2nd order and two polynoms of 1st order
514 // take both polynoms of 1st order, calculate difference of their parameters
515 // add covariance matrixes and calculate chi2 of this difference
516 // if this chi2 is bigger than a given threshold, assume that the current cluster is
517 // a kink an goto next padrow
518 // if not:
46e89793 519 // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
10757ee9 520 //
521 // write debug information to 'TPCSelectorDebug.root'
522 // only for every kDeltaWriteDebugStream'th padrow to reduce data volume
523 // and to avoid redundant data
524 //
525
9f0beaf7 526 static TLinearFitter fFitterParY(3,"pol2"); //
527 static TLinearFitter fFitterParZ(3,"pol2"); //
46e89793 528 static TLinearFitter fFitterParYWeight(3,"pol2"); //
529 static TLinearFitter fFitterParZWeight(3,"pol2"); //
530 fFitterParY.StoreData(kTRUE);
531 fFitterParZ.StoreData(kTRUE);
532 fFitterParYWeight.StoreData(kTRUE);
533 fFitterParZWeight.StoreData(kTRUE);
ae28e92e 534 if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
46e89793 535 const Int_t kDelta = 10; // delta rows to fit
536 const Float_t kMinRatio = 0.75; // minimal ratio
537 const Float_t kChi2Cut = 10.; // cut chi2 - left right
538 const Float_t kSigmaCut = 3.; //sigma cut
539 const Float_t kdEdxCut=300;
540 const Float_t kPtCut=0.300;
541
542 if (track->GetTPCsignal()>kdEdxCut) return;
543 if (TMath::Abs(AliTracker::GetBz())>0.1 &&TMath::Abs(track->Pt())<kPtCut) return;
544
545 // estimate mean error
546 Int_t nTrackletsAll = 0; // number of tracklets for given track
547 Float_t csigmaY = 0; // mean sigma for tracklet refit in Y direction
548 Float_t csigmaZ = 0; // mean sigma for tracklet refit in Z direction
549 Int_t nClusters = 0; // working variable, number of clusters per tracklet
550 Int_t sectorG = -1; // working variable, sector of tracklet, has to stay constant for one tracklet
551 Double_t refX=0;
552 // ---------------------------------------------------------------------
553 for (Int_t irow = 0; irow < 159; irow++){
554 // loop over all rows along the track
555 // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
556 // calculate mean chi^2 for this track-fit in Y and Z direction
557 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
558 if (!cluster0) continue; // no cluster found
559 Int_t sector = cluster0->GetDetector();
560
561 Int_t ipad = TMath::Nint(cluster0->GetPad());
562 Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
563 fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
564
565 if (sector != sectorG){
566 // track leaves sector before it crossed enough rows to fit / initialization
567 nClusters = 0;
9f0beaf7 568 fFitterParY.ClearPoints();
569 fFitterParZ.ClearPoints();
46e89793 570 sectorG = sector;
571 }
572 else {
573 nClusters++;
574 if (refX<1) refX=cluster0->GetX()+kDelta*0.5;
575 Double_t x = cluster0->GetX()-refX;
576 fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
577 fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
578 //
579 if ( nClusters >= kDelta + 3 ){
580 // if more than 13 (kDelta+3) clusters were added to the fitters
581 // fit the tracklet, increase trackletCounter
582 fFitterParY.Eval();
583 fFitterParZ.Eval();
584 nTrackletsAll++;
585 csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
586 csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
587 nClusters = -1;
588 fFitterParY.ClearPoints();
589 fFitterParZ.ClearPoints();
590 refX=0;
10757ee9 591 }
46e89793 592 }
593 } // for (Int_t irow = 0; irow < 159; irow++)
594 // mean chi^2 for all tracklet fits in Y and in Z direction:
595 csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
596 csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
597 // ---------------------------------------------------------------------
598 //
599 //
10757ee9 600
46e89793 601 for (Int_t irow = kDelta; irow < 159-kDelta; irow++){
602 // loop again over all rows along the track
603 // do analysis
604 //
605 Int_t nclFound = 0; // number of clusters in the neighborhood
606 Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster
607 Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster
608 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
609 if (!cluster0) continue;
610 Int_t sector = cluster0->GetDetector();
611 Float_t xref = cluster0->GetX();
612
613 // Make Fit
614 fFitterParY.ClearPoints();
615 fFitterParZ.ClearPoints();
616 fFitterParYWeight.ClearPoints();
617 fFitterParZWeight.ClearPoints();
618 // fit tracklet (clusters in given padrow +- kDelta padrows)
619 // with polynom of 2nd order and two polynoms of 1st order
620 // take both polynoms of 1st order, calculate difference of their parameters
621 // add covariance matrixes and calculate chi2 of this difference
622 // if this chi2 is bigger than a given threshold, assume that the current cluster is
623 // a kink an goto next padrow
624 AliRieman riemanFit(2*kDelta);
625 AliRieman riemanFitW(2*kDelta);
626 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
627 // loop over irow +- kDelta rows (neighboured rows)
628 //
629 //
630 if (idelta == 0) continue; // don't use center cluster
631 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
632 AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
633 if (!currentCluster) continue;
634 if (currentCluster->GetType() < 0) continue;
635 if (currentCluster->GetDetector() != sector) continue;
636 nclFound++;
637 if (idelta < 0){
638 ncl0++;
9f0beaf7 639 }
46e89793 640 if (idelta > 0){
641 ncl1++;
10757ee9 642 }
46e89793 643 riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ);
644 riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY*TMath::Sqrt(1+TMath::Abs(idelta)),csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta)));
645 } // loop over neighbourhood for fitter filling
646 if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
647 if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
648 if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow
649 riemanFit.Update();
650 riemanFitW.Update();
651 Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5));
652 Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5));
653 Double_t paramYR[3], paramZR[3];
654 Double_t paramYRW[3], paramZRW[3];
655 //
656 paramYR[0] = riemanFit.GetYat(xref);
657 paramZR[0] = riemanFit.GetZat(xref);
658 paramYRW[0] = riemanFitW.GetYat(xref);
659 paramZRW[0] = riemanFitW.GetZat(xref);
660 //
661 paramYR[1] = riemanFit.GetDYat(xref);
662 paramZR[1] = riemanFit.GetDZat(xref);
663 paramYRW[1] = riemanFitW.GetDYat(xref);
664 paramZRW[1] = riemanFitW.GetDZat(xref);
665 //
666 Int_t reject=0;
667 if (chi2R>kChi2Cut) reject+=1;
668 if (chi2RW>kChi2Cut) reject+=2;
669 if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4;
670 if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8;
671 if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16;
672 if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32;
673 //
674 TTreeSRedirector *cstream = GetDebugStreamer();
675 // get fit parameters from pol2 fit:
676 Double_t tracky = paramYR[0];
677 Double_t trackz = paramZR[0];
678 Float_t deltay = cluster0->GetY()-tracky;
679 Float_t deltaz = cluster0->GetZ()-trackz;
680 Float_t angley = paramYR[1];
681 Float_t anglez = paramZR[1];
682
683 Int_t padSize = 0; // short pads
684 if (cluster0->GetDetector() >= 36) {
10757ee9 685 padSize = 1; // medium pads
686 if (cluster0->GetRow() > 63) padSize = 2; // long pads
46e89793 687 }
688 Int_t ipad = TMath::Nint(cluster0->GetPad());
689 if (cstream){
690 Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
691 (*cstream)<<"Resol0"<<
692 "run="<<fRun<< // run number
693 "event="<<fEvent<< // event number
694 "time="<<fTime<< // time stamp of event
695 "trigger="<<fTrigger<< // trigger
696 "mag="<<fMagF<< // magnetic field
697 "padSize="<<padSize<<
698 //
699 "reject="<<reject<<
700 "cl.="<<cluster0<< // cluster info
701 "xref="<<xref<< // cluster reference X
702 //rieman fit
703 "yr="<<paramYR[0]<< // track position y - rieman fit
704 "zr="<<paramZR[0]<< // track position z - rieman fit
705 "yrW="<<paramYRW[0]<< // track position y - rieman fit - weight
706 "zrW="<<paramZRW[0]<< // track position z - rieman fit - weight
707 "dyr="<<paramYR[1]<< // track position y - rieman fit
708 "dzr="<<paramZR[1]<< // track position z - rieman fit
709 "dyrW="<<paramYRW[1]<< // track position y - rieman fit - weight
710 "dzrW="<<paramZRW[1]<< // track position z - rieman fit - weight
711 //
712 "angley="<<angley<< // angle par fit
713 "anglez="<<anglez<< // angle par fit
714 "zdr="<<zdrift<< //
715 "dy="<<deltay<<
716 "dz="<<deltaz<<
717 "sy="<<csigmaY<<
718 "sz="<<csigmaZ<<
719 "chi2R="<<chi2R<<
720 "chi2RW="<<chi2RW<<
721 "\n";
722 }
723 if (reject) continue;
724
725
726 // =========================================
727 // wirte collected information to histograms
728 // =========================================
729
730 Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
731 fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
732
733
734 TH3F * his3 = 0;
735 his3 = (TH3F*)fRMSY->At(padSize);
736 if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
737 his3 = (TH3F*)fRMSZ->At(padSize);
738 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
739
740 his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
741 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
742 his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
743 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
744
745
746 his3 = (TH3F*)fResolY->At(padSize);
747 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
748 his3 = (TH3F*)fResolZ->At(padSize);
749 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
750 his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
751 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
752 his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
753 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
754 //=============================================================================================
755 //
756 // Fill THN histograms
757 //
758 Double_t xvar[9];
759 xvar[1]=padSize;
760 xvar[2]=cluster0->GetZ();
761 xvar[3]=cluster0->GetMax();
762
763 xvar[0]=deltay;
764 xvar[4]=cluster0->GetPad()-Int_t(cluster0->GetPad());
765 xvar[5]=angley;
766 fHisDeltaY->Fill(xvar);
767 xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
768 fHisRMSY->Fill(xvar);
769
770 xvar[0]=deltaz;
771 xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin());
772 xvar[5]=anglez;
773 fHisDeltaZ->Fill(xvar);
774 xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
775 fHisRMSZ->Fill(xvar);
776
777 } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
778} // FillResolutionHistoLocal(...)
779
10757ee9 780
781
782
783
784
10757ee9 785
786
787void AliTPCcalibTracks::SetStyle() const {
788 //
789 // set style, can be called by all draw functions
790 //
791 gROOT->SetStyle("Plain");
792 gStyle->SetFillColor(10);
793 gStyle->SetPadColor(10);
794 gStyle->SetCanvasColor(10);
795 gStyle->SetStatColor(10);
796 gStyle->SetPalette(1,0);
797 gStyle->SetNumberContours(60);
798}
799
800
10757ee9 801
f78da5ae 802void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
10757ee9 803 //
804 // all functions are called, that produce pictures
805 // the histograms are written to the directory 'pathName'
806 // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
807 // 'stat' is also the number of minEntries for MakeResPlotsQTree
808 //
809
ae28e92e 810 if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
46e89793 811 MakeResPlotsQTree(stat, pathName);
10757ee9 812}
813
814
10757ee9 815
816
f78da5ae 817void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
10757ee9 818 //
819 // Make tree with resolution parameters
820 // the result is written to 'resol.root' in directory 'pathname'
821 // file information are available in fileInfo
822 // available variables in the tree 'Resol':
823 // Entries: number of entries for this resolution point
824 // nbins: number of bins that were accumulated
825 // Dim: direction, Dim==0: y-direction, Dim==1: z-direction
826 // Pad: padSize; short, medium and long
827 // Length: pad length, 0.75, 1, 1.5
828 // QMean: mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
829 // Zc: center of middle bin in drift direction
830 // Zm: mean dirftlength for accumulated Delta-Histograms
831 // Zs: width of driftlength bin
832 // AngleC: center of middle bin in Angle-Direction
833 // AngleM: mean angle for accumulated Delta-Histograms
834 // AngleS: width of Angle-bin
835 // Resol: sigma for gaus fit through Delta-Histograms
836 // Sigma: error of sigma for gaus fit through Delta Histograms
837 // MeanR: mean of the Delta-Histogram
838 // SigmaR: rms of the Delta-Histogram
839 // RMSm: mean of the gaus fit through RMS-Histogram
840 // RMS: sigma of the gaus fit through RMS-Histogram
841 // RMSe0: error of mean of gaus fit in RMS-Histogram
842 // RMSe1: error of sigma of gaus fit in RMS-Histogram
843 //
844
ae28e92e 845 if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
846 if (GetDebugLevel() > -1) cout << " relax, the calculation will take a while..." << endl;
10757ee9 847
848 gSystem->MakeDirectory(pathName);
849 gSystem->ChangeDirectory(pathName);
850 TString kFileName = "resol.root";
851 TTreeSRedirector fTreeResol(kFileName.Data());
852
853 TH3F *resArray[2][3][11];
854 TH3F *rmsArray[2][3][11];
855
856 // load histograms from fArrayQDY and fArrayQDZ
857 // into resArray and rmsArray
858 // that is all we need here
859 for (Int_t idim = 0; idim < 2; idim++){
860 for (Int_t ipad = 0; ipad < 3; ipad++){
861 for (Int_t iq = 0; iq <= 10; iq++){
862 rmsArray[idim][ipad][iq]=0;
863 resArray[idim][ipad][iq]=0;
864 Int_t bin = GetBin(iq,ipad);
865 TH3F *hresl = 0;
866 if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
867 if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
868 if (!hresl) continue;
869 resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
870 resArray[idim][ipad][iq]->SetDirectory(0);
871 TH3F * hreslRMS = 0;
872 if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
873 if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
874 if (!hreslRMS) continue;
875 rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
876 rmsArray[idim][ipad][iq]->SetDirectory(0);
877 }
878 }
879 }
ae28e92e 880 if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
10757ee9 881
882 //--------------------------------------------------------------------------------------------
883
884 char name[200];
885 Double_t qMean;
886 Double_t zMean, angleMean, zCenter, angleCenter;
887 Double_t zSigma, angleSigma;
888 TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
889 TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
890 TF1 *fitFunction = new TF1("fitFunction", "gaus");
891 Float_t entriesQ = 0;
892 Int_t loopCounter = 1;
893
894 for (Int_t idim = 0; idim < 2; idim++){
895 // Loop y-z corrdinate
896 for (Int_t ipad = 0; ipad < 3; ipad++){
897 // loop pad type
898 for (Int_t iq = -1; iq < 10; iq++){
899 // LOOP Q
ae28e92e 900 if (GetDebugLevel() > -1)
10757ee9 901 cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
902 << (Int_t)((loopCounter)/66.*100) << "% done), "
903 << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << " \r" << std::flush;
904 loopCounter++;
905 TH3F *hres = 0;
906 TH3F *hrms = 0;
907 qMean = 0;
908 entriesQ = 0;
909
910 // calculate qMean
911 if (iq == -1){
912 // integrated spectra
913 for (Int_t iql = 0; iql < 10; iql++){
914 Int_t bin = GetBin(iql,ipad);
915 TH3F *hresl = resArray[idim][ipad][iql];
916 TH3F *hrmsl = rmsArray[idim][ipad][iql];
917 if (!hresl) continue;
918 if (!hrmsl) continue;
919 entriesQ += hresl->GetEntries();
920 qMean += hresl->GetEntries() * GetQ(bin);
921 if (!hres) {
922 hres = (TH3F*)hresl->Clone();
923 hrms = (TH3F*)hrmsl->Clone();
924 }
925 else{
926 hres->Add(hresl);
927 hrms->Add(hrmsl);
928 }
929 }
930 qMean /= entriesQ;
931 qMean *= -1.; // integral mean charge
932 }
933 else {
934 // loop over neighboured Q-bins
935 // accumulate entries from neighboured Q-bins
936 for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
937 if (iql < 0) continue;
938 Int_t bin = GetBin(iql,ipad);
939 TH3F * hresl = resArray[idim][ipad][iql];
940 TH3F * hrmsl = rmsArray[idim][ipad][iql];
941 if (!hresl) continue;
942 if (!hrmsl) continue;
943 entriesQ += hresl->GetEntries();
944 qMean += hresl->GetEntries() * GetQ(bin);
945 if (!hres) {
946 hres = (TH3F*) hresl->Clone();
947 hrms = (TH3F*) hrmsl->Clone();
948 }
949 else{
950 hres->Add(hresl);
951 hrms->Add(hrmsl);
952 }
953 }
954 qMean/=entriesQ;
955 }
f69a0aa7 956 if (!hres) continue;
957 if (!hrms) continue;
958
10757ee9 959 TAxis *xAxisDriftLength = hres->GetXaxis(); // driftlength / z - axis
960 TAxis *yAxisAngle = hres->GetYaxis(); // angle axis
961 TAxis *zAxisDelta = hres->GetZaxis(); // delta axis
962 TAxis *zAxisRms = hrms->GetZaxis(); // rms axis
963
964 // loop over all angle bins
965 for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
966 angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
967 // loop over all driftlength bins
968 for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
969 zCenter = xAxisDriftLength->GetBinCenter(ibinxDL);
970 zSigma = xAxisDriftLength->GetBinWidth(ibinxDL);
971 angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
972 zMean = zCenter; // changens, when more statistic is accumulated
973 angleMean = angleCenter; // changens, when more statistic is accumulated
974
975 // create 2 1D-Histograms, projectionRes and projectionRms
976 // these histograms are delta histograms for given direction, padSize, chargeBin,
977 // angleBin and driftLengthBin
978 // later on they will be fitted with a gausian, its sigma is the resoltuion...
4aa37f93 979 snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
10757ee9 980 // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
981 projectionRes->SetNameTitle(name, name);
4aa37f93 982 snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
10757ee9 983 // TH1D * projectionRms = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
984 projectionRms->SetNameTitle(name, name);
985
986 projectionRes->Reset();
987 projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
988 projectionRms->Reset();
989 projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
990 projectionRes->SetDirectory(0);
991 projectionRms->SetDirectory(0);
992
993 Double_t entries = 0;
994 Int_t nbins = 0; // counts, how many bins were accumulated
995 zMean = 0;
996 angleMean = 0;
997 entries = 0;
998
999 // fill projectionRes and projectionRms for given dim, ipad and iq,
1000 // as well as for given angleBin and driftlengthBin
1001 // if this gives not enough statistic, include neighbourhood
1002 // (angle and driftlength) successifely
1003 for (Int_t dbin = 0; dbin <= 8; dbin++){ // delta-bins around centered angleBin and driftlengthBin
1004 for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) { // delta-bins in angle direction
1005 for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
1006 if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; // add each bin only one time !
1007 Int_t binx2 = ibinxDL + dbinx2; // position variable in x (driftlength) direction
1008 Int_t biny2 = ibinyAngle + dbiny2; // position variable in y (angle) direction
1009 if (binx2 < 1 || biny2 < 1) continue; // don't go out of the histogram!
1010 if (binx2 >= xAxisDriftLength->GetNbins()) continue; // don't go out of the histogram!
1011 if (biny2 >= yAxisAngle->GetNbins()) continue; // don't go out of the histogram!
1012 nbins++; // count the number of accumulated bins
1013 // Fill resolution histo
1014 for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
1015 // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3); // unused variable
1016 projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
1017 entries += hres->GetBinContent(binx2, biny2, ibin3);
1018 zMean += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
1019 angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
1020 } // ibin3 loop
1021 // fill RMS histo
1022 for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
1023 projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
1024 }
1025 } //dbinx2 loop
1026 if (entries > minEntries) break; // enough statistic accumulated
1027 } // dbiny2 loop
1028 if (entries > minEntries) break; // enough statistic accumulated
1029 } // dbin loop
1030 if ( entries< minEntries) continue; // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
1031 zMean /= entries;
1032 angleMean /= entries;
1033
1034 if (entries > minEntries) {
1035 // when enough statistic is accumulated
1036 // fit Delta histograms with a gausian
1037 // of the gausian is the resolution (resol), its fit error is sigma
1038 // store also mean and RMS of the histogram
1039 Float_t xmin = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1040 Float_t xmax = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1041
1042// projectionRes->Fit("gaus", "q0", "", xmin, xmax);
1043// Float_t resol = projectionRes->GetFunction("gaus")->GetParameter(2);
1044// Float_t sigma = projectionRes->GetFunction("gaus")->GetParError(2);
1045 fitFunction->SetMaximum(xmax);
1046 fitFunction->SetMinimum(xmin);
1047 projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
1048 Float_t resol = fitFunction->GetParameter(2);
1049 Float_t sigma = fitFunction->GetParError(2);
1050
1051 Float_t meanR = projectionRes->GetMean();
1052 Float_t sigmaR = projectionRes->GetRMS();
1053 // fit also RMS histograms with a gausian
1054 // store mean and sigma of the gausian in rmsMean and rmsSigma
1055 // store also the fit errors in errorRMS and errorSigma
1056 xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1057 xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1058
1059// projectionRms->Fit("gaus","q0","",xmin,xmax);
1060// Float_t rmsMean = projectionRms->GetFunction("gaus")->GetParameter(1);
1061// Float_t rmsSigma = projectionRms->GetFunction("gaus")->GetParameter(2);
1062// Float_t errorRMS = projectionRms->GetFunction("gaus")->GetParError(1);
1063// Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
1064 projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
1065 Float_t rmsMean = fitFunction->GetParameter(1);
1066 Float_t rmsSigma = fitFunction->GetParameter(2);
1067 Float_t errorRMS = fitFunction->GetParError(1);
1068 Float_t errorSigma = fitFunction->GetParError(2);
1069
1070 Float_t length = 0.75;
1071 if (ipad == 1) length = 1;
1072 if (ipad == 2) length = 1.5;
1073
1074 fTreeResol<<"Resol"<<
1075 "Entries="<<entries<< // number of entries for this resolution point
1076 "nbins="<<nbins<< // number of bins that were accumulated
1077 "Dim="<<idim<< // direction, Dim==0: y-direction, Dim==1: z-direction
1078 "Pad="<<ipad<< // padSize; short, medium and long
1079 "Length="<<length<< // pad length, 0.75, 1, 1.5
1080 "QMean="<<qMean<< // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1081 "Zc="<<zCenter<< // center of middle bin in drift direction
1082 "Zm="<<zMean<< // mean dirftlength for accumulated Delta-Histograms
1083 "Zs="<<zSigma<< // width of driftlength bin
1084 "AngleC="<<angleCenter<< // center of middle bin in Angle-Direction
1085 "AngleM="<<angleMean<< // mean angle for accumulated Delta-Histograms
1086 "AngleS="<<angleSigma<< // width of Angle-bin
1087 "Resol="<<resol<< // sigma for gaus fit through Delta-Histograms
1088 "Sigma="<<sigma<< // error of sigma for gaus fit through Delta Histograms
1089 "MeanR="<<meanR<< // mean of the Delta-Histogram
1090 "SigmaR="<<sigmaR<< // rms of the Delta-Histogram
1091 "RMSm="<<rmsMean<< // mean of the gaus fit through RMS-Histogram
1092 "RMSs="<<rmsSigma<< // sigma of the gaus fit through RMS-Histogram
1093 "RMSe0="<<errorRMS<< // error of mean of gaus fit in RMS-Histogram
1094 "RMSe1="<<errorSigma<< // error of sigma of gaus fit in RMS-Histogram
1095 "\n";
ae28e92e 1096 if (GetDebugLevel() > 5) {
10757ee9 1097 projectionRes->SetDirectory(fTreeResol.GetFile());
1098 projectionRes->Write(projectionRes->GetName());
1099 projectionRes->SetDirectory(0);
1100 projectionRms->SetDirectory(fTreeResol.GetFile());
1101 projectionRms->Write(projectionRms->GetName());
1102 projectionRes->SetDirectory(0);
1103 }
1104 } // if (projectionRes->GetSum() > minEntries)
1105 } // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
1106 } // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
1107
1108 } // iq-loop
1109 } // ipad-loop
1110 } // idim-loop
1111 delete projectionRes;
1112 delete projectionRms;
1113
1114// TFile resolFile(fTreeResol.GetFile());
1115 TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
1116 fileInfo.Write("fileInfo");
1117// resolFile.Close();
1118// fTreeResol.GetFile()->Close();
ae28e92e 1119 if (GetDebugLevel() > -1) cout << endl;
1120 if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
10757ee9 1121 gSystem->ChangeDirectory("..");
1122}
1123
1124
10757ee9 1125
10757ee9 1126
1127
1128Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
1129 //
1130 // function to merge several AliTPCcalibTracks objects after PROOF calculation
1131 // The object's histograms are merged via their merge functions
1132 // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
1133 //
1134
ae28e92e 1135 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl;
10757ee9 1136 if (!collectionList) return 0;
1137 if (collectionList->IsEmpty()) return -1;
1138
ae28e92e 1139 if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
1140 if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
10757ee9 1141 collectionList->Print();
1142
1143 // create a list for each data member
1144 TList* deltaYList = new TList;
1145 TList* deltaZList = new TList;
1146 TList* arrayAmpRowList = new TList;
1147 TList* rejectedTracksList = new TList;
10757ee9 1148 TList* clusterCutHistoList = new TList;
1149 TList* arrayAmpList = new TList;
1150 TList* arrayQDYList = new TList;
1151 TList* arrayQDZList = new TList;
1152 TList* arrayQRMSYList = new TList;
1153 TList* arrayQRMSZList = new TList;
10757ee9 1154 TList* resolYList = new TList;
1155 TList* resolZList = new TList;
1156 TList* rMSYList = new TList;
1157 TList* rMSZList = new TList;
1158
1159// TList* nRowsList = new TList;
1160// TList* nSectList = new TList;
1161// TList* fileNoList = new TList;
1162
1163 TIterator *listIterator = collectionList->MakeIterator();
1164 AliTPCcalibTracks *calibTracks = 0;
ae28e92e 1165 if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
10757ee9 1166 Int_t counter = 0;
a4c5fca9 1167 while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
10757ee9 1168 // loop over all entries in the collectionList and get dataMembers into lists
10757ee9 1169
10757ee9 1170 arrayQDYList->Add(calibTracks->GetfArrayQDY());
1171 arrayQDZList->Add(calibTracks->GetfArrayQDZ());
1172 arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
1173 arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
1174 resolYList->Add(calibTracks->GetfResolY());
1175 resolZList->Add(calibTracks->GetfResolZ());
1176 rMSYList->Add(calibTracks->GetfRMSY());
1177 rMSZList->Add(calibTracks->GetfRMSZ());
10757ee9 1178 rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
1179 clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
a4c5fca9 1180 //
1181 if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
1182 fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
1183 // fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
10757ee9 1184 counter++;
ae28e92e 1185 if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
af6a50bb 1186 AddHistos(calibTracks);
10757ee9 1187 }
1188
10757ee9 1189
1190 // merge data members
ae28e92e 1191 if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
10757ee9 1192 fClusterCutHisto->Merge(clusterCutHistoList);
1193 fRejectedTracksHisto->Merge(rejectedTracksList);
10757ee9 1194
1195 TObjArray* objarray = 0;
1196 TH1* hist = 0;
1197 TList* histList = 0;
1198 TIterator *objListIterator = 0;
1199
46e89793 1200
ae28e92e 1201 if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
10757ee9 1202 // merge fArrayQDY
1203 for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
1204 objListIterator = arrayQDYList->MakeIterator();
1205 histList = new TList;
1206 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1207 // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1208 hist = (TH3F*)(objarray->At(i));
1209 histList->Add(hist);
1210 }
1211 ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
1212 delete histList;
1213 delete objListIterator;
1214 }
1215
ae28e92e 1216 if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
10757ee9 1217 // merge fArrayQDZ
1218 for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1219 objListIterator = arrayQDZList->MakeIterator();
1220 histList = new TList;
1221 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1222 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
10757ee9 1223 hist = (TH3F*)(objarray->At(i));
1224 histList->Add(hist);
1225 }
1226 ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
1227 delete histList;
1228 delete objListIterator;
1229 }
1230
ae28e92e 1231 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
10757ee9 1232 // merge fArrayQRMSY
1233 for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
1234 objListIterator = arrayQRMSYList->MakeIterator();
1235 histList = new TList;
1236 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1237 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
03533457 1238 // if (!objarray) continue; // removed for coverity -> JMT
10757ee9 1239 hist = (TH3F*)(objarray->At(i));
1240 histList->Add(hist);
1241 }
1242 ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
1243 delete histList;
1244 delete objListIterator;
1245 }
1246
ae28e92e 1247 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
10757ee9 1248 // merge fArrayQRMSZ
1249 for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1250 objListIterator = arrayQRMSZList->MakeIterator();
1251 histList = new TList;
1252 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1253 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1254 if (!objarray) continue;
1255 hist = (TH3F*)(objarray->At(i));
1256 histList->Add(hist);
1257 }
1258 ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
1259 delete histList;
1260 delete objListIterator;
1261 }
1262
10757ee9 1263
10757ee9 1264
1265
1266
1267
ae28e92e 1268 if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
10757ee9 1269 // merge fResolY
1270 for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
1271 objListIterator = resolYList->MakeIterator();
1272 histList = new TList;
1273 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1274 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1275 if (!objarray) continue;
1276 hist = (TH3F*)(objarray->At(i));
1277 histList->Add(hist);
1278 }
1279 ((TH3F*)(fResolY->At(i)))->Merge(histList);
1280 delete histList;
1281 delete objListIterator;
1282 }
1283
1284 // merge fResolZ
1285 for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1286 objListIterator = resolZList->MakeIterator();
1287 histList = new TList;
1288 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1289 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1290 if (!objarray) continue;
1291 hist = (TH3F*)(objarray->At(i));
1292 histList->Add(hist);
1293 }
1294 ((TH3F*)(fResolZ->At(i)))->Merge(histList);
1295 delete histList;
1296 delete objListIterator;
1297 }
1298
1299 // merge fRMSY
1300 for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
1301 objListIterator = rMSYList->MakeIterator();
1302 histList = new TList;
1303 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1304 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1305 if (!objarray) continue;
1306 hist = (TH3F*)(objarray->At(i));
1307 histList->Add(hist);
1308 }
1309 ((TH3F*)(fRMSY->At(i)))->Merge(histList);
1310 delete histList;
1311 delete objListIterator;
1312 }
1313
1314 // merge fRMSZ
1315 for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1316 objListIterator = rMSZList->MakeIterator();
1317 histList = new TList;
1318 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1319 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1320 if (!objarray) continue;
1321 hist = (TH3F*)(objarray->At(i));
1322 histList->Add(hist);
1323 }
1324 ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
1325 delete histList;
1326 delete objListIterator;
1327 }
1328
1329 delete deltaYList;
1330 delete deltaZList;
1331 delete arrayAmpRowList;
1332 delete arrayAmpList;
1333 delete arrayQDYList;
1334 delete arrayQDZList;
1335 delete arrayQRMSYList;
1336 delete arrayQRMSZList;
1337 delete resolYList;
1338 delete resolZList;
1339 delete rMSYList;
1340 delete rMSZList;
10757ee9 1341 delete listIterator;
1342
ae28e92e 1343 if (GetDebugLevel() > 0) cout << "merging done!" << endl;
10757ee9 1344
1345 return 1;
1346}
1347
1348
10757ee9 1349
1350
af6a50bb 1351void AliTPCcalibTracks::MakeHistos(){
1352 //
1353 ////make THnSparse
1354 //
1355 //THnSparse *fHisDeltaY; // THnSparse - delta Y
1356 //THnSparse *fHisDeltaZ; // THnSparse - delta Z
1357 //THnSparse *fHisRMSY; // THnSparse - rms Y
1358 //THnSparse *fHisRMSZ; // THnSparse - rms Z
1359 //THnSparse *fHisQmax; // THnSparse - qmax
1360 //THnSparse *fHisQtot; // THnSparse - qtot
1361 // cluster performance bins
1362 // 0 - variable of interest
1363 // 1 - pad type - 0- short 1-medium 2-long pads
1364 // 2 - drift length - drift length -0-1
1365 // 3 - Qmax - Qmax - 2- 400
1366 // 4 - cog - COG position - 0-1
1367 // 5 - tan(phi) - local y angle
1368 // 6 - tan(theta) - local z angle
1369 // 7 - sector - sector number
1370 Double_t xminTrack[8], xmaxTrack[8];
1371 Int_t binsTrack[8];
1372 TString axisName[8];
1373
1374 //
46e89793 1375 binsTrack[0]=200;
af6a50bb 1376 axisName[0] ="var";
1377
1378 binsTrack[1] =3;
1379 xminTrack[1] =0; xmaxTrack[1]=3;
1380 axisName[1] ="pad type";
1381 //
46e89793 1382 binsTrack[2] =20;
1383 xminTrack[2] =-250; xmaxTrack[2]=250;
1384 axisName[2] ="z";
af6a50bb 1385 //
1386 binsTrack[3] =10;
1387 xminTrack[3] =1; xmaxTrack[3]=400;
1388 axisName[3] ="Qmax";
1389 //
46e89793 1390 binsTrack[4] =20;
af6a50bb 1391 xminTrack[4] =0; xmaxTrack[4]=1;
1392 axisName[4] ="cog";
1393 //
46e89793 1394 binsTrack[5] =15;
1395 xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
1396 axisName[5] ="tan(angle)";
af6a50bb 1397 //
af6a50bb 1398 //
46e89793 1399 xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1400 fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1401 xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1402 fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1403 xminTrack[0] =0.; xmaxTrack[0]=0.5;
46e89793 1404 fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1405 xminTrack[0] =0.; xmaxTrack[0]=0.5;
46e89793 1406 fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1407 xminTrack[0] =0.; xmaxTrack[0]=100;
46e89793 1408 fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1409
1410 xminTrack[0] =0.; xmaxTrack[0]=250;
46e89793 1411 fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1412
1413
1414 for (Int_t ivar=0;ivar<6;ivar++){
1415 fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1416 fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1417 fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1418 fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1419 fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1420 fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1421 fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
1422 fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1423 fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1424 fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1425 fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1426 fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1427 }
1428
1429
af6a50bb 1430 BinLogX(fHisDeltaY,3);
1431 BinLogX(fHisDeltaZ,3);
1432 BinLogX(fHisRMSY,3);
1433 BinLogX(fHisRMSZ,3);
1434 BinLogX(fHisQmax,3);
1435 BinLogX(fHisQtot,3);
1436
1437}
1438
1439void AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
1440 //
1441 // Add histograms
1442 //
1443 if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
1444 if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
1445 if (calib->fHisRMSY) fHisRMSY->Add(calib->fHisRMSY);
1446 if (calib->fHisRMSZ) fHisRMSZ->Add(calib->fHisRMSZ);
1447}
46e89793 1448
1449
1450
1451void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
1452 //
1453 // Dump summary info
1454 //
1455 // 0.OBJ: TAxis var
1456 // 1.OBJ: TAxis pad type
1457 // 2.OBJ: TAxis z
1458 // 3.OBJ: TAxis Qmax
1459 // 4.OBJ: TAxis cog
1460 // 5.OBJ: TAxis tan(angle)
1461 //
1462 if (ptype>3) return;
1463 Int_t idim[6]={0,1,2,3,4,5};
1464 TString hname[4]={"dy","dz","rmsy","rmsz"};
1465 //
1466 Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
1467 Int_t first5=hisInput->GetAxis(5)->GetFirst();
1468 Int_t last5 =hisInput->GetAxis(5)->GetLast();
1469 //
1470 for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){ // axis 5 - local angle
1471 Bool_t bin5All=kTRUE;
1472 if (ibin5>=first5){
1473 hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
1474 if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
1475 bin5All=kFALSE;
1476 }
1477 Double_t x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
1478 THnSparse * his5= hisInput->Projection(5,idim); //projected histogram according selection 5
1479 Int_t nbins4=his5->GetAxis(4)->GetNbins();
1480 Int_t first4=his5->GetAxis(4)->GetFirst();
1481 Int_t last4 =his5->GetAxis(4)->GetLast();
1482
1483 for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){ // axis 4 - cog
1484 Bool_t bin4All=kTRUE;
1485 if (ibin4>=first4){
1486 his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
1487 if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
1488 bin4All=kFALSE;
1489 }
1490 Double_t x4= his5->GetAxis(4)->GetBinCenter(ibin4);
1491 THnSparse * his4= his5->Projection(4,idim); //projected histogram according selection 4
1492 //
1493 Int_t nbins3=his4->GetAxis(3)->GetNbins();
1494 Int_t first3=his4->GetAxis(3)->GetFirst();
1495 Int_t last3 =his4->GetAxis(3)->GetLast();
1496 //
1497 for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - Qmax
1498 Bool_t bin3All=kTRUE;
1499 if (ibin3>=first3){
1500 his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
1501 bin3All=kFALSE;
1502 }
1503 Double_t x3= his4->GetAxis(3)->GetBinCenter(ibin3);
1504 THnSparse * his3= his4->Projection(3,idim); //projected histogram according selection 3
1505 //
1506 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1507 Int_t first2 = his3->GetAxis(2)->GetFirst();
1508 Int_t last2 = his3->GetAxis(2)->GetLast();
1509 //
1510 for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){ // axis 2 - z
1511 Bool_t bin2All=kTRUE;
1512 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1513 if (ibin2>=first2){
1514 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
1515 if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
1516 bin2All=kFALSE;
1517 }
1518 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
1519 //
1520 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
1521 Int_t first1 = his2->GetAxis(1)->GetFirst();
1522 Int_t last1 = his2->GetAxis(1)->GetLast();
1523 for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){ //axis 1 - pad type
1524 Bool_t bin1All=kTRUE;
1525 if (ibin1>=first1){
1526 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1527 bin1All=kFALSE;
1528 }
1529 Double_t x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
1530 TH1 * hisDelta = his2->Projection(0);
1531 Double_t entries = hisDelta->GetEntries();
1532 Double_t mean=0, rms=0;
1533 if (entries>10){
1534 mean = hisDelta->GetMean();
1535 rms = hisDelta->GetRMS();
1536 hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
1537 mean = hisDelta->GetMean();
1538 rms = hisDelta->GetRMS();
1539 }
1540 //
1541 (*pcstream)<<hname[ptype].Data()<<
1542 // flag - intgrated values for given bin
1543 "angleA="<<bin5All<<
1544 "cogA="<<bin4All<<
1545 "qmaxA="<<bin3All<<
1546 "zA="<<bin2All<<
1547 "ipadA="<<bin1All<<
1548 // center of bin value
1549 "angle="<<x5<<
1550 "cog="<<x4<<
1551 "qmax="<<x3<<
1552 "z="<<x2<<
1553 "ipad="<<x1<<
1554 // mean values
1555 //
1556 "entries="<<entries<<
1557 "mean="<<mean<<
1558 "rms="<<rms<<
1559 "ptype="<<ptype<< //
1560 "\n";
1561 delete hisDelta;
1562 printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);
1563 }//loop z
1564 delete his2;
1565 } // loop Q max
1566 delete his3;
1567 } // loop COG
1568 delete his4;
1569 }//loop local angle
1570 delete his5;
1571 }
1572}