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