]> 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
Clean-up of include files
[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
1223 if (!objarray) continue;
1224 hist = (TH3F*)(objarray->At(i));
1225 histList->Add(hist);
1226 }
1227 ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
1228 delete histList;
1229 delete objListIterator;
1230 }
1231
ae28e92e 1232 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
10757ee9 1233 // merge fArrayQRMSY
1234 for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
1235 objListIterator = arrayQRMSYList->MakeIterator();
1236 histList = new TList;
1237 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1238 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1239 if (!objarray) continue;
1240 hist = (TH3F*)(objarray->At(i));
1241 histList->Add(hist);
1242 }
1243 ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
1244 delete histList;
1245 delete objListIterator;
1246 }
1247
ae28e92e 1248 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
10757ee9 1249 // merge fArrayQRMSZ
1250 for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1251 objListIterator = arrayQRMSZList->MakeIterator();
1252 histList = new TList;
1253 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1254 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1255 if (!objarray) continue;
1256 hist = (TH3F*)(objarray->At(i));
1257 histList->Add(hist);
1258 }
1259 ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
1260 delete histList;
1261 delete objListIterator;
1262 }
1263
10757ee9 1264
10757ee9 1265
1266
1267
1268
ae28e92e 1269 if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
10757ee9 1270 // merge fResolY
1271 for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
1272 objListIterator = resolYList->MakeIterator();
1273 histList = new TList;
1274 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1275 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1276 if (!objarray) continue;
1277 hist = (TH3F*)(objarray->At(i));
1278 histList->Add(hist);
1279 }
1280 ((TH3F*)(fResolY->At(i)))->Merge(histList);
1281 delete histList;
1282 delete objListIterator;
1283 }
1284
1285 // merge fResolZ
1286 for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1287 objListIterator = resolZList->MakeIterator();
1288 histList = new TList;
1289 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1290 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1291 if (!objarray) continue;
1292 hist = (TH3F*)(objarray->At(i));
1293 histList->Add(hist);
1294 }
1295 ((TH3F*)(fResolZ->At(i)))->Merge(histList);
1296 delete histList;
1297 delete objListIterator;
1298 }
1299
1300 // merge fRMSY
1301 for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
1302 objListIterator = rMSYList->MakeIterator();
1303 histList = new TList;
1304 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1305 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1306 if (!objarray) continue;
1307 hist = (TH3F*)(objarray->At(i));
1308 histList->Add(hist);
1309 }
1310 ((TH3F*)(fRMSY->At(i)))->Merge(histList);
1311 delete histList;
1312 delete objListIterator;
1313 }
1314
1315 // merge fRMSZ
1316 for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1317 objListIterator = rMSZList->MakeIterator();
1318 histList = new TList;
1319 while (( objarray = (TObjArray*)objListIterator->Next() )) {
1320 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1321 if (!objarray) continue;
1322 hist = (TH3F*)(objarray->At(i));
1323 histList->Add(hist);
1324 }
1325 ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
1326 delete histList;
1327 delete objListIterator;
1328 }
1329
1330 delete deltaYList;
1331 delete deltaZList;
1332 delete arrayAmpRowList;
1333 delete arrayAmpList;
1334 delete arrayQDYList;
1335 delete arrayQDZList;
1336 delete arrayQRMSYList;
1337 delete arrayQRMSZList;
1338 delete resolYList;
1339 delete resolZList;
1340 delete rMSYList;
1341 delete rMSZList;
10757ee9 1342 delete listIterator;
1343
ae28e92e 1344 if (GetDebugLevel() > 0) cout << "merging done!" << endl;
10757ee9 1345
1346 return 1;
1347}
1348
1349
10757ee9 1350
1351
af6a50bb 1352void AliTPCcalibTracks::MakeHistos(){
1353 //
1354 ////make THnSparse
1355 //
1356 //THnSparse *fHisDeltaY; // THnSparse - delta Y
1357 //THnSparse *fHisDeltaZ; // THnSparse - delta Z
1358 //THnSparse *fHisRMSY; // THnSparse - rms Y
1359 //THnSparse *fHisRMSZ; // THnSparse - rms Z
1360 //THnSparse *fHisQmax; // THnSparse - qmax
1361 //THnSparse *fHisQtot; // THnSparse - qtot
1362 // cluster performance bins
1363 // 0 - variable of interest
1364 // 1 - pad type - 0- short 1-medium 2-long pads
1365 // 2 - drift length - drift length -0-1
1366 // 3 - Qmax - Qmax - 2- 400
1367 // 4 - cog - COG position - 0-1
1368 // 5 - tan(phi) - local y angle
1369 // 6 - tan(theta) - local z angle
1370 // 7 - sector - sector number
1371 Double_t xminTrack[8], xmaxTrack[8];
1372 Int_t binsTrack[8];
1373 TString axisName[8];
1374
1375 //
46e89793 1376 binsTrack[0]=200;
af6a50bb 1377 axisName[0] ="var";
1378
1379 binsTrack[1] =3;
1380 xminTrack[1] =0; xmaxTrack[1]=3;
1381 axisName[1] ="pad type";
1382 //
46e89793 1383 binsTrack[2] =20;
1384 xminTrack[2] =-250; xmaxTrack[2]=250;
1385 axisName[2] ="z";
af6a50bb 1386 //
1387 binsTrack[3] =10;
1388 xminTrack[3] =1; xmaxTrack[3]=400;
1389 axisName[3] ="Qmax";
1390 //
46e89793 1391 binsTrack[4] =20;
af6a50bb 1392 xminTrack[4] =0; xmaxTrack[4]=1;
1393 axisName[4] ="cog";
1394 //
46e89793 1395 binsTrack[5] =15;
1396 xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
1397 axisName[5] ="tan(angle)";
af6a50bb 1398 //
af6a50bb 1399 //
46e89793 1400 xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1401 fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1402 xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1403 fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1404 xminTrack[0] =0.; xmaxTrack[0]=0.5;
46e89793 1405 fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1406 xminTrack[0] =0.; xmaxTrack[0]=0.5;
46e89793 1407 fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1408 xminTrack[0] =0.; xmaxTrack[0]=100;
46e89793 1409 fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
af6a50bb 1410
1411 xminTrack[0] =0.; xmaxTrack[0]=250;
46e89793 1412 fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1413
1414
1415 for (Int_t ivar=0;ivar<6;ivar++){
1416 fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1417 fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1418 fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1419 fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1420 fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1421 fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1422 fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
1423 fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1424 fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1425 fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1426 fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1427 fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1428 }
1429
1430
af6a50bb 1431 BinLogX(fHisDeltaY,3);
1432 BinLogX(fHisDeltaZ,3);
1433 BinLogX(fHisRMSY,3);
1434 BinLogX(fHisRMSZ,3);
1435 BinLogX(fHisQmax,3);
1436 BinLogX(fHisQtot,3);
1437
1438}
1439
1440void AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
1441 //
1442 // Add histograms
1443 //
1444 if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
1445 if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
1446 if (calib->fHisRMSY) fHisRMSY->Add(calib->fHisRMSY);
1447 if (calib->fHisRMSZ) fHisRMSZ->Add(calib->fHisRMSZ);
1448}
46e89793 1449
1450
1451
1452void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
1453 //
1454 // Dump summary info
1455 //
1456 // 0.OBJ: TAxis var
1457 // 1.OBJ: TAxis pad type
1458 // 2.OBJ: TAxis z
1459 // 3.OBJ: TAxis Qmax
1460 // 4.OBJ: TAxis cog
1461 // 5.OBJ: TAxis tan(angle)
1462 //
1463 if (ptype>3) return;
1464 Int_t idim[6]={0,1,2,3,4,5};
1465 TString hname[4]={"dy","dz","rmsy","rmsz"};
1466 //
1467 Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
1468 Int_t first5=hisInput->GetAxis(5)->GetFirst();
1469 Int_t last5 =hisInput->GetAxis(5)->GetLast();
1470 //
1471 for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){ // axis 5 - local angle
1472 Bool_t bin5All=kTRUE;
1473 if (ibin5>=first5){
1474 hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
1475 if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
1476 bin5All=kFALSE;
1477 }
1478 Double_t x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
1479 THnSparse * his5= hisInput->Projection(5,idim); //projected histogram according selection 5
1480 Int_t nbins4=his5->GetAxis(4)->GetNbins();
1481 Int_t first4=his5->GetAxis(4)->GetFirst();
1482 Int_t last4 =his5->GetAxis(4)->GetLast();
1483
1484 for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){ // axis 4 - cog
1485 Bool_t bin4All=kTRUE;
1486 if (ibin4>=first4){
1487 his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
1488 if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
1489 bin4All=kFALSE;
1490 }
1491 Double_t x4= his5->GetAxis(4)->GetBinCenter(ibin4);
1492 THnSparse * his4= his5->Projection(4,idim); //projected histogram according selection 4
1493 //
1494 Int_t nbins3=his4->GetAxis(3)->GetNbins();
1495 Int_t first3=his4->GetAxis(3)->GetFirst();
1496 Int_t last3 =his4->GetAxis(3)->GetLast();
1497 //
1498 for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - Qmax
1499 Bool_t bin3All=kTRUE;
1500 if (ibin3>=first3){
1501 his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
1502 bin3All=kFALSE;
1503 }
1504 Double_t x3= his4->GetAxis(3)->GetBinCenter(ibin3);
1505 THnSparse * his3= his4->Projection(3,idim); //projected histogram according selection 3
1506 //
1507 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
1508 Int_t first2 = his3->GetAxis(2)->GetFirst();
1509 Int_t last2 = his3->GetAxis(2)->GetLast();
1510 //
1511 for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){ // axis 2 - z
1512 Bool_t bin2All=kTRUE;
1513 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1514 if (ibin2>=first2){
1515 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
1516 if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
1517 bin2All=kFALSE;
1518 }
1519 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
1520 //
1521 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
1522 Int_t first1 = his2->GetAxis(1)->GetFirst();
1523 Int_t last1 = his2->GetAxis(1)->GetLast();
1524 for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){ //axis 1 - pad type
1525 Bool_t bin1All=kTRUE;
1526 if (ibin1>=first1){
1527 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1528 bin1All=kFALSE;
1529 }
1530 Double_t x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
1531 TH1 * hisDelta = his2->Projection(0);
1532 Double_t entries = hisDelta->GetEntries();
1533 Double_t mean=0, rms=0;
1534 if (entries>10){
1535 mean = hisDelta->GetMean();
1536 rms = hisDelta->GetRMS();
1537 hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
1538 mean = hisDelta->GetMean();
1539 rms = hisDelta->GetRMS();
1540 }
1541 //
1542 (*pcstream)<<hname[ptype].Data()<<
1543 // flag - intgrated values for given bin
1544 "angleA="<<bin5All<<
1545 "cogA="<<bin4All<<
1546 "qmaxA="<<bin3All<<
1547 "zA="<<bin2All<<
1548 "ipadA="<<bin1All<<
1549 // center of bin value
1550 "angle="<<x5<<
1551 "cog="<<x4<<
1552 "qmax="<<x3<<
1553 "z="<<x2<<
1554 "ipad="<<x1<<
1555 // mean values
1556 //
1557 "entries="<<entries<<
1558 "mean="<<mean<<
1559 "rms="<<rms<<
1560 "ptype="<<ptype<< //
1561 "\n";
1562 delete hisDelta;
1563 printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);
1564 }//loop z
1565 delete his2;
1566 } // loop Q max
1567 delete his3;
1568 } // loop COG
1569 delete his4;
1570 }//loop local angle
1571 delete his5;
1572 }
1573}