Fixed memory leaks (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTracksGain.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////////////////////////////////////////////////////////////////////////////
0d3279d4 17//
18//
19// Gain calibration using tracks
20//
21// The main goal:
22// 1.) Inner TPC gain alignement - (parabolic) parameterization (inside of one sector)
23// 2.) Angular and z-position correction (parabolic) parameterization
24// 3.) Sector gain alignment
25//
26// Following histograms are accumulated
27// a.) Simple 1D histograms per chamber
28// b.) Profile histograms per chamber - local x dependence
29// c.) 2D Profile histograms - local x - fi dependence
30//
31// To get the gain map - the simple solution - use the histograms - is not enough
32// The resulting mean amplitude map depends strongly on the track topology
33// These dependence can be reduced, taking into account angular effect, and diffusion effect
34// Using proper fit modeles
35//
36//
37//
10757ee9 38//
39// === Calibration class for gain calibration using tracks ===
40//
41// 1) Genereal idea
42// ================
43// A 6-parametric parabolic function
44//
45// G(x, y) = p0 + p1*x + p2*y + p3*x^2 + p4*y^2 + p5 * x*y
46//
47// is fitted to the maximum charge values or total charge values of
48// all the clusters contained in the tracks that are added to this
49// object. This fit is performed for each read out chamber, in fact even
50// for each type of pad sizes (thus for one segment, which consists of
51// an IROC and an OROC, there are three fitters used, corresponding to
52// the three pad sizes). The coordinate origin is at the center of the
53// particular pad size region on each ROC.
54//
55// Because of the Landau nature of the charge deposition we use
56// different "types" of fitters instead of one to minimize the effect
57// of the long Landau tail. The difference between the fitters is only
58// the charge value, that is put into them, i.e. the charge is subject
59// to a transformation. At this point we use three different fit types:
60//
61// a) simple: the charge is put in as it is
62// b) sqrt: the square root of the charge is put into the fitter
63// c) log: fgkM * Log(1+q/fgkM) is put into the fitter, with
64// q being the untransformed charge and fgkM=25
65//
66// The results of the fits may be visualized and further used by
67// creating an AliTPCCalROC or AliTPCCalPad. You may specify to undo
68// the transformation and/or to normalize to the pad size.
69//
70// Not every track you add to this object is actually used for
71// calibration. There are some cuts and conditions to exclude bad
72// tracks, e.g. a pt cut to cut out tracks with too much charge
73// deposition or a cut on edge clusters which are not fully
74// registered and don't give a usable signal.
75//
76// 2) Interface / usage
77// ====================
78// For each track to be added you need to call Process().
79// This method expects an AliTPCseed, which contains the necessary
80// cluster information. At the moment of writing this information
81// is stored in an AliESDfriend corresponding to an AliESD.
82// You may also call AddTrack() if you don't want the cuts and
83// other quality conditions to kick in (thus forcing the object to
84// accept the track) or AddCluster() for adding single clusters.
85// Call one of the Evaluate functions to evaluate the fitter(s) and
86// to retrieve the fit parameters, erros and so on. You can also
87// do this later on by using the different Getters.
88//
89// The visualization methods CreateFitCalPad() and CreateFitCalROC()
90// are straight forward to use.
91//
92// Note: If you plan to write this object to a ROOT file, make sure
93// you evaluate all the fitters *before* writing, because due
94// to a bug in the fitter component writing fitters doesn't
95// work properly (yet). Be aware that you cannot re-evaluate
96// the fitters after loading this object from file.
97// (This will be gone for a new ROOT version > v5-17-05)
98//
99////////////////////////////////////////////////////////////////////////////
100
101#include <TPDGCode.h>
102#include <TStyle.h>
103#include "TSystem.h"
104#include "TMatrixD.h"
105#include "TTreeStream.h"
106#include "TF1.h"
107#include "AliTPCParamSR.h"
108#include "AliTPCClusterParam.h"
109#include "AliTrackPointArray.h"
110#include "TCint.h"
111#include "AliTPCcalibTracksGain.h"
112#include <TH1.h>
113#include <TH3F.h>
114#include <TLinearFitter.h>
115#include <TTreeStream.h>
116#include <TFile.h>
117#include <TCollection.h>
118#include <TIterator.h>
119
120//
121// AliRoot includes
122//
123#include "AliMagF.h"
124#include "AliMathBase.h"
125//
126#include "AliTPCROC.h"
127#include "AliTPCParamSR.h"
128#include "AliTPCCalROC.h"
129#include "AliTPCCalPad.h"
8076baa0 130#include "AliTPCClusterParam.h"
10757ee9 131//
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 "AliTPCcalibTracksCuts.h"
140#include "AliTPCFitPad.h"
141
142// REMOVE ALL OF THIS
143#include <TTree.h>
144#include "AliESDEvent.h"
145
146/*
147
148TFile f("TPCCalibTracksGain.root")
149
150gSystem->Load("libPWG1.so")
151AliTreeDraw comp
152comp.SetTree(dEdx)
153Double_t chi2;
154TVectorD vec(3)
155TMatrixD mat(3,3)
156TString * str = comp.FitPlane("Cl.fQ/dedxQ.fElements[0]","Cl.fY++Cl.fX","Cl.fDetector<36",chi2,vec,mat)
157
158*/
159
160ClassImp(AliTPCcalibTracksGain)
161
162const Bool_t AliTPCcalibTracksGain::fgkUseTotalCharge = kTRUE;
163const Double_t AliTPCcalibTracksGain::fgkM = 25.;
164const char* AliTPCcalibTracksGain::fgkDebugStreamFileName = "TPCCalibTracksGain.root";
165AliTPCParamSR* AliTPCcalibTracksGain::fgTPCparam = new AliTPCParamSR();
166
167AliTPCcalibTracksGain::AliTPCcalibTracksGain() :
b8601924 168 AliTPCcalibBase(),
0d3279d4 169 fCuts(0), // cuts that are used for sieving the tracks used for calibration
170 //
171 // Simple Array of histograms
172 //
173 fArrayQM(0), // Qmax normalized
174 fArrayQT(0), // Qtot normalized
175 fProfileArrayQM(0), // Qmax normalized versus local X
176 fProfileArrayQT(0), // Qtot normalized versus local X
177 fProfileArrayQM2D(0), // Qmax normalized versus local X and phi
178 fProfileArrayQT2D(0), // Qtot normalized versus local X and phi
179 //
180 // Fitters
181 //
182 fSimpleFitter(0), // simple fitter for short pads
183 fSqrtFitter(0), // sqrt fitter for medium pads
184 fLogFitter(0), // log fitter for long pads
185 fFitter0M(0), // fitting of the atenuation, angular correction, and mean chamber gain
186 fFitter1M(0), // fitting of the atenuation, angular correction, and mean chamber gain
187 fFitter2M(0), // fitting of the atenuation, angular correction, and mean chamber gain
188 fFitter0T(0), // fitting of the atenuation, angular correction, and mean chamber gain
189 fFitter1T(0), // fitting of the atenuation, angular correction, and mean chamber gain
190 fFitter2T(0), // fitting of the atenuation, angular correction, and mean chamber gain
191 fSingleSectorFitter(0), // just for debugging
192 //
193 // Counters
194 //
195 fTotalTracks(0), // just for debugging
196 fAcceptedTracks(0), // just for debugging
197 fDebugCalPadRaw(0), // just for debugging
198 fDebugCalPadCorr(0), // just for debugging
199 fPrevIter(0) // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!)
200
10757ee9 201{
202 //
203 // Default constructor.
204 //
205}
206
207AliTPCcalibTracksGain::AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj) :
b8601924 208 AliTPCcalibBase(obj),
0d3279d4 209 fCuts(obj.fCuts), // cuts that are used for sieving the tracks used for calibration
f1c2a4a3 210 fArrayQM(0), // Qmax normalized
211 fArrayQT(0), // Qtot normalized
0d3279d4 212 //
213 // Simple Histograms
214 //
215 fProfileArrayQM(obj.fProfileArrayQM), // Qmax normalized versus local X
216 fProfileArrayQT(obj.fProfileArrayQT), // Qtot normalized versus local X
217 fProfileArrayQM2D(obj.fProfileArrayQM2D), // Qmax normalized versus local X and phi
218 fProfileArrayQT2D(obj.fProfileArrayQT2D), // Qtot normalized versus local X and phi
219 //
220 // Fitters
221 //
222 fSimpleFitter(obj.fSimpleFitter), // simple fitter for short pads
223 fSqrtFitter(obj.fSqrtFitter), // sqrt fitter for medium pads
224 fLogFitter(obj.fLogFitter), // log fitter for long pads
225 fFitter0M(obj.fFitter0M),
226 fFitter1M(obj.fFitter1M),
227 fFitter2M(obj.fFitter2M),
228 fFitter0T(obj.fFitter0T),
229 fFitter1T(obj.fFitter1T),
230 fFitter2T(obj.fFitter2T),
231 fSingleSectorFitter(obj.fSingleSectorFitter), // just for debugging
232 //
233 // Counters
234 //
235 fTotalTracks(obj.fTotalTracks), // just for debugging
236 fAcceptedTracks(obj.fAcceptedTracks), // just for debugging
237 fDebugCalPadRaw(obj.fDebugCalPadRaw), // just for debugging
238 fDebugCalPadCorr(obj.fDebugCalPadCorr), // just for debugging
239 fPrevIter(obj.fPrevIter) // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!)
240
10757ee9 241{
242 //
243 // Copy constructor.
244 //
10757ee9 245}
246
247AliTPCcalibTracksGain& AliTPCcalibTracksGain::operator=(const AliTPCcalibTracksGain& rhs) {
248 //
249 // Assignment operator.
250 //
251
252 if (this != &rhs) {
253 TNamed::operator=(rhs);
254 fDebugCalPadRaw = new AliTPCCalPad(*(rhs.fDebugCalPadRaw));
255 fDebugCalPadCorr = new AliTPCCalPad(*(rhs.fDebugCalPadCorr));
256 fSimpleFitter = new AliTPCFitPad(*(rhs.fSimpleFitter));
257 fSqrtFitter = new AliTPCFitPad(*(rhs.fSqrtFitter));
258 fLogFitter = new AliTPCFitPad(*(rhs.fLogFitter));
259 fSingleSectorFitter = new AliTPCFitPad(*(rhs.fSingleSectorFitter));
260 fPrevIter = new AliTPCcalibTracksGain(*(rhs.fPrevIter));
10757ee9 261 fCuts = new AliTPCcalibTracksCuts(*(rhs.fCuts));
262 }
263 return *this;
264}
265
f1c2a4a3 266AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts, TNamed* /*debugStreamPrefix*/, AliTPCcalibTracksGain* prevIter) :
b8601924 267 AliTPCcalibBase(),
0d3279d4 268 fCuts(0), // cuts that are used for sieving the tracks used for calibration
f1c2a4a3 269 fArrayQM(0), // Qmax normalized
270 fArrayQT(0), // Qtot normalized
0d3279d4 271 //
272 // Simple Histograms
273 //
274 fProfileArrayQM(0), // Qmax normalized versus local X
275 fProfileArrayQT(0), // Qtot normalized versus local X
276 fProfileArrayQM2D(0), // Qmax normalized versus local X and phi
277 fProfileArrayQT2D(0), // Qtot normalized versus local X and phi
278 //
279 // Fitters
280 //
281 fSimpleFitter(0), // simple fitter for short pads
282 fSqrtFitter(0), // sqrt fitter for medium pads
283 fLogFitter(0), // log fitter for long pads
284 fFitter0M(0), // fitting of the atenuation, angular correction, and mean chamber gain
285 fFitter1M(0), // fitting of the atenuation, angular correction, and mean chamber gain
286 fFitter2M(0), // fitting of the atenuation, angular correction, and mean chamber gain
287 fFitter0T(0), // fitting of the atenuation, angular correction, and mean chamber gain
288 fFitter1T(0), // fitting of the atenuation, angular correction, and mean chamber gain
289 fFitter2T(0), // fitting of the atenuation, angular correction, and mean chamber gain
290 fSingleSectorFitter(0), // just for debugging
291 //
292 // Counters
293 //
294 fTotalTracks(0), // just for debugging
295 fAcceptedTracks(0), // just for debugging
296 fDebugCalPadRaw(0), // just for debugging
297 fDebugCalPadCorr(0), // just for debugging
298 fPrevIter(0) // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!)
299
10757ee9 300{
301 //
302 // Constructor.
0d3279d4 303 //
10757ee9 304 G__SetCatchException(0);
b8601924 305 this->SetNameTitle(name, title);
10757ee9 306 fCuts = cuts;
10757ee9 307 fPrevIter = prevIter;
0d3279d4 308 //
309 // Fitter initialization
310 //
10757ee9 311 fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
312 fSqrtFitter = new AliTPCFitPad(8, "hyp7", "");
313 fLogFitter = new AliTPCFitPad(8, "hyp7", "");
314 fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
8076baa0 315 //
316 fFitter0M = new TLinearFitter(45,"hyp44");
317 fFitter1M = new TLinearFitter(45,"hyp44");
318 fFitter2M = new TLinearFitter(45,"hyp44");
319 fFitter0T = new TLinearFitter(45,"hyp44");
320 fFitter1T = new TLinearFitter(45,"hyp44");
321 fFitter2T = new TLinearFitter(45,"hyp44");
0d3279d4 322 //
cbc19295 323 //
324 fFitter0M->StoreData(kFALSE);
325 fFitter1M->StoreData(kFALSE);
326 fFitter2M->StoreData(kFALSE);
327 fFitter0T->StoreData(kFALSE);
328 fFitter1T->StoreData(kFALSE);
329 fFitter2T->StoreData(kFALSE);
330 //
0d3279d4 331 //
332 // Add profile histograms -JUST for visualization - Not used for real calibration
333 //
334 //
335 fArrayQM=new TObjArray(73); // Qmax normalized
336 fArrayQT=new TObjArray(73); // Qtot normalized
337 fProfileArrayQM = new TObjArray(37); // Qmax normalized versus local X
338 fProfileArrayQT = new TObjArray(37); // Qtot normalized versus local X
339 fProfileArrayQM2D = new TObjArray(37); // Qmax normalized versus local X and phi
340 fProfileArrayQT2D = new TObjArray(37); // Qtot normalized versus local X and phi
341 char hname[1000];
342 for (Int_t i=0; i<73; i++){
343 sprintf(hname,"QM_%d",i);
344 fArrayQM->AddAt(new TH1F(hname,hname,200,0,1000),i);
345 sprintf(hname,"QT_%d",i);
346 fArrayQT->AddAt(new TH1F(hname,hname,200,0,1000),i);
347 }
10757ee9 348
0d3279d4 349 for (Int_t i=0; i<37;i++){
350 sprintf(hname,"QMvsx_%d",i);
351 fProfileArrayQM->AddAt(new TProfile(hname,hname,50,89,250),i);
352 sprintf(hname,"QTvsx_%d",i);
353 fProfileArrayQT->AddAt(new TProfile(hname,hname,50,89,250),i);
354 sprintf(hname,"QM2D_%d",i);
355 fProfileArrayQM2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
356 sprintf(hname,"QT2D_%d",i);
357 fProfileArrayQT2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
358 }
359 //
360 // just for debugging -counters
361 //
10757ee9 362 fTotalTracks = 0;
363 fAcceptedTracks = 0;
364 fDebugCalPadRaw = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
0d3279d4 365 fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
10757ee9 366 // this will be gone for the a new ROOT version > v5-17-05
367 for (UInt_t i = 0; i < 36; i++) {
368 fNShortClusters[i] = 0;
369 fNMediumClusters[i] = 0;
370 fNLongClusters[i] = 0;
371 }
372 }
373
374AliTPCcalibTracksGain::~AliTPCcalibTracksGain() {
375 //
376 // Destructor.
377 //
378
379 Info("Destructor","");
380 if (fSimpleFitter) delete fSimpleFitter;
381 if (fSqrtFitter) delete fSqrtFitter;
382 if (fLogFitter) delete fLogFitter;
383 if (fSingleSectorFitter) delete fSingleSectorFitter;
384
10757ee9 385 if (fDebugCalPadRaw) delete fDebugCalPadRaw;
386 if (fDebugCalPadCorr) delete fDebugCalPadCorr;
387}
388
389void AliTPCcalibTracksGain::Terminate(){
390 //
391 // Evaluate fitters and close the debug stream.
392 // Also move or copy the debug stream, if a debugStreamPrefix is provided.
393 //
394
395 Evaluate();
ae28e92e 396 AliTPCcalibBase::Terminate();
10757ee9 397}
398
10757ee9 399
10757ee9 400
401void AliTPCcalibTracksGain::Process(AliTPCseed* seed) {
402 //
403 // Main method to be called when a new seed is supposed to be processed
404 // and be used for gain calibration. Its quality is checked before it
405 // is added.
406 //
407
408 fTotalTracks++;
b8601924 409 if (!fCuts->AcceptTrack(seed)) return;
10757ee9 410 fAcceptedTracks++;
411 AddTrack(seed);
412}
413
414Long64_t AliTPCcalibTracksGain::Merge(TCollection *list) {
415 //
416 // Merge() merges the results of all AliTPCcalibTracksGain objects contained in
417 // list, thus allowing a distributed computation of several files, e.g. on PROOF.
418 // The merged results are merged with the data members of the AliTPCcalibTracksGain
419 // object used for calling the Merge method.
420 // The return value is 0 /*the total number of tracks used for calibration*/ if the merge
421 // is successful, otherwise it is -1.
422 //
423
424 if (!list || list->IsEmpty()) return -1;
425
426 if (!fSimpleFitter) fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
427 if (!fSqrtFitter) fSqrtFitter = new AliTPCFitPad(8, "hyp7", "");
428 if (!fLogFitter) fLogFitter = new AliTPCFitPad(8, "hyp7", "");
429 if (!fSingleSectorFitter) fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
430
10757ee9 431
432 // just for debugging
0d3279d4 433 if (!fDebugCalPadRaw) fDebugCalPadRaw = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
434 if (!fDebugCalPadCorr) fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
10757ee9 435
436 TIterator* iter = list->MakeIterator();
437 AliTPCcalibTracksGain* cal = 0;
438
439 while ((cal = (AliTPCcalibTracksGain*)iter->Next())) {
440 if (!cal->InheritsFrom(AliTPCcalibTracksGain::Class())) {
441 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
442 return -1;
443 }
0d3279d4 444
10757ee9 445 Add(cal);
446 }
447 return 0;
448}
449
450void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
451 //
452 // Adds another AliTPCcalibTracksGain object to this object.
453 //
454
455 fSimpleFitter->Add(cal->fSimpleFitter);
456 fSqrtFitter->Add(cal->fSqrtFitter);
457 fLogFitter->Add(cal->fLogFitter);
458 fSingleSectorFitter->Add(cal->fSingleSectorFitter);
0d3279d4 459 //
460 //
461 //
462 fFitter0M->Add(cal->fFitter0M);
463 fFitter1M->Add(cal->fFitter1M);
464 fFitter2M->Add(cal->fFitter2M);
465 fFitter0T->Add(cal->fFitter0T);
466 fFitter1T->Add(cal->fFitter1T);
467 fFitter2T->Add(cal->fFitter2T);
468 //
469 //
470 // histograms
471 //
472 for (Int_t i=0; i<73; i++){
473 TH1F *his,*hism;
474 his = (TH1F*)fArrayQM->At(i);
475 hism = (TH1F*)cal->fArrayQM->At(i);
476 if (his && hism) his->Add(hism);
477 his = (TH1F*)fArrayQT->At(i);
478 hism = (TH1F*)cal->fArrayQT->At(i);
479 if (his && hism) his->Add(hism);
480 }
481 //
482 //
483 for (Int_t i=0; i<37; i++){
484 TProfile *his,*hism;
485 his = (TProfile*)fProfileArrayQM->At(i);
486 hism = (TProfile*)cal->fProfileArrayQM->At(i);
487 if (his && hism) his->Add(hism);
488 his = (TProfile*)fProfileArrayQT->At(i);
489 hism = (TProfile*)cal->fProfileArrayQT->At(i);
490 if (his && hism) his->Add(hism);
491 }
492 //
493 //
494 for (Int_t i=0; i<37; i++){
495 TProfile2D *his,*hism;
496 his = (TProfile2D*)fProfileArrayQM2D->At(i);
497 hism = (TProfile2D*)cal->fProfileArrayQM2D->At(i);
498 if (his && hism) his->Add(hism);
499 his = (TProfile2D*)fProfileArrayQT2D->At(i);
500 hism = (TProfile2D*)cal->fProfileArrayQT2D->At(i);
501 if (his && hism) his->Add(hism);
502 }
503 //
10757ee9 504 // this will be gone for the a new ROOT version > v5-17-05
505 for (UInt_t iSegment = 0; iSegment < 36; iSegment++) {
506 fNShortClusters[iSegment] += cal->fNShortClusters[iSegment];
507 fNMediumClusters[iSegment] += cal->fNMediumClusters[iSegment];
508 fNLongClusters[iSegment] += cal->fNLongClusters[iSegment];
509 }
510
511 // just for debugging, remove me
512 fTotalTracks += cal->fTotalTracks;
513 fAcceptedTracks += cal->fAcceptedTracks;
514 fDebugCalPadRaw->Add(cal->fDebugCalPadRaw);
515 fDebugCalPadCorr->Add(cal->fDebugCalPadCorr);
516
10757ee9 517}
518
519void AliTPCcalibTracksGain::AddTrack(AliTPCseed* seed) {
520 //
521 // The clusters making up the track (seed) are added to various fit functions.
522 // See AddCluster(...) for more detail.
523 //
524
0d3279d4 525 DumpTrack(seed);
526 //
527 // simple histograming part
528 for (Int_t i=0; i<159; i++){
529 AliTPCclusterMI* cluster = seed->GetClusterPointer(i);
530 if (cluster) AddCluster(cluster);
531 }
10757ee9 532}
533
0d3279d4 534void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster){
535 //
536 // Adding cluster information to the simple histograms
537 // No correction, fittings are applied
538 //
539 Float_t kThreshold=5;
540 if (cluster->GetX()<=0) return;
541 if (cluster->GetQ()<=kThreshold) return;
542 //
543
544
545 Int_t sector = cluster->GetDetector();
546 TH1F * his;
547 his = GetQT(sector);
548 if (his) his->Fill(cluster->GetQ());
549 his = GetQT(-1);
550 if (his) his->Fill(cluster->GetQ());
551 his = GetQM(sector);
552 if (his) his->Fill(cluster->GetMax());
553 his = GetQM(-1);
554 if (his) his->Fill(cluster->GetMax());
555 //
556 sector = sector%36;
557 TProfile *prof;
558 prof = GetProfileQT(sector);
559 if (prof) prof->Fill(cluster->GetX(),cluster->GetQ());
560 prof = GetProfileQT(-1);
561 if (prof) prof->Fill(cluster->GetX(),cluster->GetQ());
562 prof = GetProfileQM(sector);
563 if (prof) prof->Fill(cluster->GetX(),cluster->GetMax());
564 prof = GetProfileQM(-1);
565 if (prof) prof->Fill(cluster->GetX(),cluster->GetMax());
566 //
567 Float_t phi = cluster->GetY()/cluster->GetX();
568 TProfile2D *prof2;
569 prof2 = GetProfileQT2D(sector);
570 if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetQ());
571 prof2 = GetProfileQT2D(-1);
572 if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetQ());
573 prof2 = GetProfileQM2D(sector);
574 if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetMax());
575 prof2 = GetProfileQM2D(-1);
576 if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetMax());
577
578 //
579}
580
581
582
f1c2a4a3 583void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t /*momenta*/, Float_t/* mdedx*/, Int_t padType,
584 Float_t xcenter, TVectorD& dedxQ, TVectorD& /*dedxM*/, Float_t /*fraction*/, Float_t fraction2, Float_t dedge,
585 TVectorD& /*parY*/, TVectorD& /*parZ*/, TVectorD& meanPos) {
10757ee9 586 //
587 // Adds cluster to the appropriate fitter for later analysis.
588 // The charge used for the fit is the maximum charge for this specific cluster or the
589 // accumulated charge per cluster, depending on the value of fgkUseTotalCharge.
590 // Depending on the pad size where the cluster is registered, the value will be put in
591 // the appropriate fitter. Furthermore, for each pad size three different types of fitters
592 // are used. The fit functions are the same for all fitters (parabolic functions), but the value
593 // added to each fitter is different. The simple fitter gets the charge plugged in as is, the sqrt fitter
594 // gets the square root of the charge, and the log fitter gets fgkM*(1+q/fgkM), where q is the original charge
595 // and fgkM==25.
596 //
597
598 if (!cluster) {
599 Error("AddCluster", "Cluster not valid.");
600 return;
601 }
602
603 if (dedge < 3.) return;
604 if (fraction2 > 0.7) return;
605
606 //Int_t padType = GetPadType(cluster->GetX());
607 Double_t xx[7];
608 //Double_t centerPad[2] = {0};
609 //AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
610 //xx[0] = cluster->GetX() - centerPad[0];
611 //xx[1] = cluster->GetY() - centerPad[1];
612 xx[0] = cluster->GetX() - xcenter;
613 xx[1] = cluster->GetY();
614 xx[2] = xx[0] * xx[0];
615 xx[3] = xx[1] * xx[1];
616 xx[4] = xx[0] * xx[1];
617 xx[5] = TMath::Abs(cluster->GetZ()) - TMath::Abs(meanPos[4]);
618 xx[6] = xx[5] * xx[5];
0d3279d4 619 //
620 // Update profile histograms
621 //
10757ee9 622
0d3279d4 623 //
624 // Update fitters
625 //
10757ee9 626 Int_t segment = cluster->GetDetector() % 36;
627 Double_t q = fgkUseTotalCharge ? ((Double_t)(cluster->GetQ())) : ((Double_t)(cluster->GetMax())); // note: no normalization to pad size!
628
629 // just for debugging
630 Int_t row = 0;
631 Int_t pad = 0;
632 GetRowPad(cluster->GetX(), cluster->GetY(), row, pad);
633 fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
634
635 // correct charge by normalising to mean charge per track
636 q /= dedxQ[2];
637
638 // just for debugging
639 fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
640
641 Double_t sqrtQ = TMath::Sqrt(q);
642 Double_t logQ = fgkM * TMath::Log(1 + q / fgkM);
643 fSimpleFitter->GetFitter(segment, padType)->AddPoint(xx, q);
644 fSqrtFitter->GetFitter(segment, padType)->AddPoint(xx, sqrtQ);
645 fLogFitter->GetFitter(segment, padType)->AddPoint(xx, logQ);
646 fSingleSectorFitter->GetFitter(0, padType)->AddPoint(xx, q);
647
648 // this will be gone for the a new ROOT version > v5-17-05
649 if (padType == kShortPads)
650 fNShortClusters[segment]++;
0d3279d4 651 if (padType == kMediumPads)
10757ee9 652 fNMediumClusters[segment]++;
0d3279d4 653 if (padType == kLongPads)
10757ee9 654 fNLongClusters[segment]++;
655}
656
657void AliTPCcalibTracksGain::Evaluate(Bool_t robust, Double_t frac) {
658 //
659 // Evaluates all fitters contained in this object.
660 // If the robust option is set to kTRUE a robust fit is performed with frac as
661 // the minimal fraction of good points (see TLinearFitter::EvalRobust for details).
662 // Beware: Robust fitting is much slower!
663 //
664
665 fSimpleFitter->Evaluate(robust, frac);
666 fSqrtFitter->Evaluate(robust, frac);
667 fLogFitter->Evaluate(robust, frac);
668 fSingleSectorFitter->Evaluate(robust, frac);
0d3279d4 669 fFitter0M->Eval();
670 fFitter1M->Eval();
671 fFitter2M->Eval();
672 fFitter0T->Eval();
673 fFitter1T->Eval();
674 fFitter2T->Eval();
10757ee9 675}
676
677AliTPCCalPad* AliTPCcalibTracksGain::CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
678 //
679 // Creates the calibration object AliTPCcalPad using fitted parameterization
680 //
681 TObjArray tpc(72);
682 for (UInt_t iSector = 0; iSector < 72; iSector++)
683 tpc.Add(CreateFitCalROC(iSector, fitType, undoTransformation, normalizeToPadSize));
684 return new AliTPCCalPad(&tpc);
685}
686
687AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
688 //
689 // Create the AliTPCCalROC with the values per pad
690 // sector - sector of interest
691 // fitType -
692 //
693
694 TVectorD par(8);
695 if (sector < 36) {
696 GetParameters(sector % 36, 0, fitType, par);
697 return CreateFitCalROC(sector, 0, par, fitType, undoTransformation, normalizeToPadSize);
698 }
699 else {
700 GetParameters(sector % 36, 1, fitType, par);
701 AliTPCCalROC* roc1 = CreateFitCalROC(sector, 1, par, fitType, undoTransformation, normalizeToPadSize);
702 GetParameters(sector % 36, 2, fitType, par);
703 AliTPCCalROC* roc2 = CreateFitCalROC(sector, 2, par, fitType, undoTransformation, normalizeToPadSize);
704 AliTPCCalROC* roc3 = CreateCombinedCalROC(roc1, roc2);
705 delete roc1;
706 delete roc2;
707 return roc3;
708 }
709}
710
711AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t padType, TVectorD &fitParam, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
712 //
713 // This function is essentially a copy of AliTPCCalROC::CreateGlobalFitCalROC(...), with the
714 // modifications, that the center of the region of same pad size is used as the origin
715 // of the fit function instead of the center of the ROC.
716 // The possibility of a linear fit is removed as well because it is not needed.
717 // Only values for pads with the given pad size are calculated, the rest is 0.
718 // Set undoTransformation for undoing the transformation that was applied to the
719 // charge values before they were put into the fitter (thus allowing comparison to the original
720 // charge values). For fitType use 0 for the simple fitter, 1 for the sqrt fitter, 2 for the log fitter.
721 // If normalizeToPadSize is true, the values are normalized to the pad size.
722 // Please be aware, that you even need to specify the fitType if you want to normalize to the pad size without
723 // undoing the transformation (because normalizing involves undoing the trafo first, then normalizing, then
724 // applying the trafo again).
725 // Please note: The normalization to the pad size is a simple linear scaling with the pad length, which
726 // actually doesn't describe reality!
727 //
728
729 Float_t dlx, dly;
730 Double_t centerPad[2] = {0};
731 Float_t localXY[3] = {0};
732 AliTPCROC* tpcROC = AliTPCROC::Instance();
733 if ((padType == 0 && sector >= tpcROC->GetNInnerSector()) || (padType > 0 && sector < tpcROC->GetNInnerSector()) || sector >= tpcROC->GetNSector())
734 return 0;
735 AliTPCCalROC* lROCfitted = new AliTPCCalROC(sector);
736 //tpcROC->GetPositionLocal(sector, lROCfitted->GetNrows()/2, lROCfitted->GetNPads(lROCfitted->GetNrows()/2)/2, centerPad); // use this instead of the switch statement if you want to calculate the center of the ROC and not the center of the regions with the same pad size
737 UInt_t startRow = 0;
738 UInt_t endRow = 0;
739 switch (padType) {
740 case kShortPads:
741 startRow = 0;
742 endRow = lROCfitted->GetNrows();
743 break;
744 case kMediumPads:
745 startRow = 0;
746 endRow = 64;
747 break;
748 case kLongPads:
749 startRow = 64;
750 endRow = lROCfitted->GetNrows();
751 break;
752 }
753
754 AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
755 Double_t value = 0;
756 for (UInt_t irow = startRow; irow < endRow; irow++) {
757 for (UInt_t ipad = 0; ipad < lROCfitted->GetNPads(irow); ipad++) {
758 tpcROC->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number
759 dlx = localXY[0] - centerPad[0];
760 dly = localXY[1] - centerPad[1];
761 value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
762
763 // Let q' = value be the transformed value without any pad size corrections,
764 // let T be the transformation and let l be the pad size
765 // 1) don't undo transformation, don't normalize: return q'
766 // 2) undo transformation, don't normalize: return T^{-1} q'
767 // 3) undo transformation, normalize: return (T^{-1} q') / l
768 // 4) don't undo transformation, normalize: return T((T^{-1} q') / l)
769 if (!undoTransformation && !normalizeToPadSize) {/* value remains unchanged */} // (1)
770 else { // (2), (3), (4)
771 //calculate T^{-1}
772 switch (fitType) {
773 case 0: /* value remains unchanged */ break;
774 case 1: value = value * value; break;
775 case 2: value = (TMath::Exp(value / fgkM) - 1) * fgkM; break;
776 default: Error("CreateFitCalROC", "Wrong fit type."); break;
777 }
778 if (normalizeToPadSize) value /= GetPadLength(localXY[0]); // (3)
779 }
780 if (!undoTransformation && normalizeToPadSize) { // (4)
781 // calculate T
782 switch (fitType) {
783 case 0: /* value remains unchanged */ break;
784 case 1: value = TMath::Sqrt(value); break;
785 case 2: value = fgkM * TMath::Log(1 + value / fgkM); break;
786 default: Error("CreateFitCalROC", "Wrong fit type."); break;
787 }
788 }
789 lROCfitted->SetValue(irow, ipad, value);
790 }
791 }
792 return lROCfitted;
793}
794
795AliTPCCalROC* AliTPCcalibTracksGain::CreateCombinedCalROC(const AliTPCCalROC* roc1, const AliTPCCalROC* roc2) {
796 //
797 // Combines the medium pad size values of roc1 with the long pad size values of roc2 into a new
798 // AliTPCCalROC. Returns a null pointer if any one of the ROCs is an IROC; issues a warning message
799 // if the sectors of roc1 and roc2 don't match, but still continue and use the sector of roc1 as the
800 // sector of the new ROC.
801 //
802
803 if (!roc1 || !roc2) return 0;
804 if ((Int_t)(roc1->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
805 if ((Int_t)(roc2->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
806 if (roc1->GetSector() != roc2->GetSector()) Warning("CreateCombinedCalROC", "Sector number mismatch.");
807 AliTPCCalROC* roc = new AliTPCCalROC(roc1->GetSector());
808
809 for (UInt_t iRow = 0; iRow < 64; iRow++) {
810 for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
811 roc->SetValue(iRow, iPad, roc1->GetValue(iRow, iPad));
812 }
813 for (UInt_t iRow = 64; iRow < roc->GetNrows(); iRow++) {
814 for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
815 roc->SetValue(iRow, iPad, roc2->GetValue(iRow, iPad));
816 }
817 return roc;
818}
819
820void AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) {
821 //
822 // Puts the fit parameters for the specified segment (IROC & OROC), padType and fitType
823 // into the fitParam TVectorD (which should contain 8 elements).
824 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
825 // Note: The fitter has to be evaluated first!
826 //
827
828 GetFitter(segment, padType, fitType)->GetParameters(fitParam);
829}
830
831void AliTPCcalibTracksGain::GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError) {
832 //
833 // Puts the fit parameter errors for the specified segment (IROC & OROC), padType and fitType
834 // into the fitError TVectorD (which should contain 8 elements).
835 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
836 // Note: The fitter has to be evaluated first!
837 //
838
839 GetFitter(segment, padType, fitType)->GetErrors(fitError);
840 fitError *= TMath::Sqrt(GetRedChi2(segment, padType, fitType));
841}
842
843Double_t AliTPCcalibTracksGain::GetRedChi2(UInt_t segment, UInt_t padType, UInt_t fitType) {
844 //
845 // Returns the reduced chi^2 value for the specified segment, padType and fitType.
846 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
847 // Note: The fitter has to be evaluated first!
848 //
849
850 // this will be gone for the a new ROOT version > v5-17-05
851 Int_t lNClusters = 0;
852 switch (padType) {
853 case kShortPads:
854 lNClusters = fNShortClusters[segment];
855 break;
856 case kMediumPads:
857 lNClusters = fNMediumClusters[segment];
858 break;
859 case kLongPads:
860 lNClusters = fNLongClusters[segment];
861 break;
862 }
863 return GetFitter(segment, padType, fitType)->GetChisquare()/(lNClusters - 8);
864}
865
866void AliTPCcalibTracksGain::GetCovarianceMatrix(UInt_t segment, UInt_t padType, UInt_t fitType, TMatrixD& covMatrix) {
867 //
868 // Returns the covariance matrix for the specified segment, padType, fitType.
869 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
870 //
871
872 GetFitter(segment, padType, fitType)->GetCovarianceMatrix(covMatrix);
873}
874
875TLinearFitter* AliTPCcalibTracksGain::GetFitter(UInt_t segment, UInt_t padType, UInt_t fitType) {
876 //
877 // Returns the TLinearFitter object for the specified segment, padType, fitType.
878 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
879 //
880
881 switch (fitType) {
882 case kSimpleFitter:
883 return fSimpleFitter->GetFitter(segment, padType);
884 case kSqrtFitter:
885 return fSqrtFitter->GetFitter(segment, padType);
886 case kLogFitter:
887 return fLogFitter->GetFitter(segment, padType);
888 case 3:
889 return fSingleSectorFitter->GetFitter(0, padType);
890 }
891 return 0;
892}
893
894Double_t AliTPCcalibTracksGain::GetPadLength(Double_t lx) {
895 //
896 // The function returns 0.75 for an IROC, 1. for an OROC at medium pad size position,
897 // 1.5 for an OROC at long pad size position, -1 if out of bounds.
898 //
899
900 Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
901 Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
902 Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
903 Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
904 Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
905 Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
906
907 // if IROC
908 if (lx >= irocLow && lx <= irocUp) return 0.75;
909 // if OROC medium pads
910 if (lx >= orocLow1 && lx <= orocUp1) return 1.;
911 // if OROC long pads
912 if (lx >= orocLow2 && lx <= orocUp2) return 1.5;
913 // if out of bounds
914 return -1;
915}
916
917Int_t AliTPCcalibTracksGain::GetPadType(Double_t lx) {
918 //
919 // The function returns 0 for an IROC, 1 for an OROC at medium pad size position,
920 // 2 for an OROC at long pad size position, -1 if out of bounds.
921 //
922
923 if (GetPadLength(lx) == 0.75) return 0;
924 else if (GetPadLength(lx) == 1.) return 1;
925 else if (GetPadLength(lx) == 1.5) return 2;
926 return -1;
927}
928
929// ONLY FOR DEBUGGING PURPOSES - REMOVE ME WHEN NOT NEEDED ANYMORE
930Bool_t AliTPCcalibTracksGain::GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad) {
931 //
932 // Calculate the row and pad number when the local coordinates are given.
933 // Returns kFALSE if the position is out of range, otherwise return kTRUE.
934 // WARNING: This function is preliminary and probably isn't very accurate!!
935 //
936
937 Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
938 //Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
939 Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
940 //Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
941 Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
942 //Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
943
944 if (GetPadType(lx) == 0) {
945 row = (Int_t)((lx - irocLow) / fgTPCparam->GetInnerPadPitchLength());
946 pad = (Int_t)((ly + fgTPCparam->GetYInner(row)) / fgTPCparam->GetInnerPadPitchWidth());
947 } else if (GetPadType(lx) == 1) {
948 row = (Int_t)((lx - orocLow1) / fgTPCparam->GetOuter1PadPitchLength());
949 pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
950 } else if (GetPadType(lx) == 2) {
951 row = fgTPCparam->GetNRowUp1() + (Int_t)((lx - orocLow2) / fgTPCparam->GetOuter2PadPitchLength());
952 pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
953 }
954 else return kFALSE;
955 return kTRUE;
956}
957
958void AliTPCcalibTracksGain::DumpTrack(AliTPCseed* track) {
959 //
960 // Dump track information to the debug stream
961 //
962
10757ee9 963 Int_t rows[200];
0d3279d4 964 Int_t sector[3];
965 Int_t npoints[3];
966 TVectorD dedxM[3];
967 TVectorD dedxQ[3];
968 TVectorD parY[3];
969 TVectorD parZ[3];
970 TVectorD meanPos[3];
971 //
972 Int_t count=0;
10757ee9 973 for (Int_t ipad = 0; ipad < 3; ipad++) {
0d3279d4 974 dedxM[ipad].ResizeTo(5);
975 dedxQ[ipad].ResizeTo(5);
976 parY[ipad].ResizeTo(3);
977 parZ[ipad].ResizeTo(3);
978 meanPos[ipad].ResizeTo(6);
979 Bool_t isOK = GetDedx(track, ipad, rows, sector[ipad], npoints[ipad], dedxM[ipad], dedxQ[ipad], parY[ipad], parZ[ipad], meanPos[ipad]);
980 if (isOK)
981 AddTracklet(sector[ipad],ipad, dedxQ[ipad], dedxM[ipad], parY[ipad], parZ[ipad], meanPos[ipad] );
982 if (isOK) count++;
10757ee9 983 }
ae28e92e 984
985 TTreeSRedirector * cstream = GetDebugStreamer();
986 if (cstream){
987 (*cstream) << "Track" <<
0d3279d4 988 "Track.=" << track << // track information
0d3279d4 989 "\n";
ae28e92e 990 //
991 //
992 //
993 if ( GetStreamLevel()>1 && count>1){
994 (*cstream) << "TrackG" <<
995 "Track.=" << track << // track information
996 "ncount="<<count<<
997 // info for pad type 0
998 "sector0="<<sector[0]<<
999 "npoints0="<<npoints[0]<<
1000 "dedxM0.="<<&dedxM[0]<<
1001 "dedxQ0.="<<&dedxQ[0]<<
1002 "parY0.="<<&parY[0]<<
1003 "parZ0.="<<&parZ[0]<<
1004 "meanPos0.="<<&meanPos[0]<<
1005 //
1006 // info for pad type 1
1007 "sector1="<<sector[1]<<
1008 "npoints1="<<npoints[1]<<
1009 "dedxM1.="<<&dedxM[1]<<
1010 "dedxQ1.="<<&dedxQ[1]<<
1011 "parY1.="<<&parY[1]<<
1012 "parZ1.="<<&parZ[1]<<
1013 "meanPos1.="<<&meanPos[1]<<
1014 //
1015 // info for pad type 2
1016 "sector2="<<sector[2]<<
1017 "npoints2="<<npoints[2]<<
1018 "dedxM2.="<<&dedxM[2]<<
1019 "dedxQ2.="<<&dedxQ[2]<<
1020 "parY2.="<<&parY[2]<<
1021 "parZ2.="<<&parZ[2]<<
1022 "meanPos2.="<<&meanPos[2]<<
1023 //
1024 "\n";
1025 }
0d3279d4 1026 }
0d3279d4 1027}
1028
f1c2a4a3 1029Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /*rows*/,
0d3279d4 1030 Int_t &sector, Int_t& npoints,
1031 TVectorD &dedxM, TVectorD &dedxQ,
1032 TVectorD &parY, TVectorD &parZ, TVectorD&meanPos)
1033{
1034 //
1035 // GetDedx for given sector for given track
1036 // padType - type of pads
1037 //
1038
1039 static TLinearFitter fitY(2, "pol1");
1040 static TLinearFitter fitZ(2, "pol1");
1041 //
1042 //
10757ee9 1043 Int_t firstRow = 0, lastRow = 0;
1044 Int_t minRow = 100;
1045 Float_t xcenter = 0;
1046 const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
1047 const Float_t kedgey = 4.;
1048 if (padType == 0) {
1049 firstRow = 0;
1050 lastRow = fgTPCparam->GetNRowLow();
1051 xcenter = 108.47;
1052 }
1053 if (padType == 1) {
1054 firstRow = fgTPCparam->GetNRowLow();
1055 lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1056 xcenter = 166.60;
1057 }
1058 if (padType == 2) {
1059 firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1060 lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
1061 xcenter = 222.6;
1062 }
1063 minRow = (lastRow - firstRow) / 2;
1064 //
1065 //
1066 Int_t nclusters = 0;
1067 Int_t nclustersNE = 0; // number of not edge clusters
1068 Int_t lastSector = -1;
1069 Float_t amplitudeQ[100];
1070 Float_t amplitudeM[100];
1071 Int_t rowIn[100];
1072 Int_t index[100];
1073 //
0d3279d4 1074 //
10757ee9 1075 fitY.ClearPoints();
1076 fitZ.ClearPoints();
10757ee9 1077
1078 for (Int_t iCluster = firstRow; iCluster < lastRow; iCluster++) {
1079 AliTPCclusterMI* cluster = track->GetClusterPointer(iCluster);
1080 if (cluster) {
1081 Int_t detector = cluster->GetDetector() ;
1082 if (lastSector == -1) lastSector = detector;
1083 if (lastSector != detector) continue;
1084 amplitudeQ[nclusters] = cluster->GetQ();
1085 amplitudeM[nclusters] = cluster->GetMax();
1086 rowIn[nclusters] = iCluster;
1087 nclusters++;
1088 Double_t dx = cluster->GetX() - xcenter;
1089 Double_t y = cluster->GetY();
1090 Double_t z = cluster->GetZ();
1091 fitY.AddPoint(&dx, y);
1092 fitZ.AddPoint(&dx, z);
1093 meanPos[0] += dx;
1094 meanPos[1] += dx;
1095 meanPos[2] += y;
1096 meanPos[3] += y*y;
1097 meanPos[4] += z;
1098 meanPos[5] += z*z;
1099 if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) nclustersNE++;
1100 }
1101 }
1102
1103 if (nclusters < minRow / 2) return kFALSE;
1104 if (nclustersNE < minRow / 2) return kFALSE;
1105 for (Int_t i = 0; i < 6; i++) meanPos[i] /= Double_t(nclusters);
1106 fitY.Eval();
1107 fitZ.Eval();
1108 fitY.GetParameters(parY);
1109 fitZ.GetParameters(parZ);
1110 //
1111 // calculate truncated mean
1112 //
1113 TMath::Sort(nclusters, amplitudeQ, index, kFALSE);
0d3279d4 1114 //
1115 //
1116 //
10757ee9 1117 Float_t ndedx[5];
1118 for (Int_t i = 0; i < 5; i++) {
1119 dedxQ[i] = 0;
1120 dedxM[i] = 0;
1121 ndedx[i] = 0;
1122 }
1123 //
1124 // dedx calculation
1125 //
1126 Int_t inonEdge = 0;
1127 for (Int_t i = 0; i < nclusters; i++) {
1128 Int_t rowSorted = rowIn[index[i]];
1129 AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1130
1131 if (TMath::Abs(cluster->GetY()) > cluster->GetX()*ktany - kedgey) continue; //don't take edge clusters
1132 inonEdge++;
1133 if (inonEdge < nclustersNE * 0.5) {
1134 ndedx[0]++;
1135 dedxQ[0] += amplitudeQ[index[i]];
1136 dedxM[0] += amplitudeM[index[i]];
1137 }
1138 if (inonEdge < nclustersNE * 0.6) {
1139 ndedx[1]++;
1140 dedxQ[1] += amplitudeQ[index[i]];
1141 dedxM[1] += amplitudeM[index[i]];
1142 }
1143 if (inonEdge < nclustersNE * 0.7) {
1144 ndedx[2]++;
1145 dedxQ[2] += amplitudeQ[index[i]];
1146 dedxM[2] += amplitudeM[index[i]];
1147 }
1148 if (inonEdge < nclustersNE * 0.8) {
1149 ndedx[3]++;
1150 dedxQ[3] += amplitudeQ[index[i]];
1151 dedxM[3] += amplitudeM[index[i]];
1152 }
1153 if (inonEdge < nclustersNE * 0.9) {
1154 ndedx[4]++;
1155 dedxQ[4] += amplitudeQ[index[i]];
1156 dedxM[4] += amplitudeM[index[i]];
1157 }
1158 }
1159 for (Int_t i = 0; i < 5; i++) {
1160 dedxQ[i] /= ndedx[i];
1161 dedxM[i] /= ndedx[i];
1162 }
ae28e92e 1163 TTreeSRedirector * cstream = GetDebugStreamer();
10757ee9 1164 inonEdge = 0;
0d3279d4 1165 Float_t momenta = track->GetP();
1166 Float_t mdedx = track->GetdEdx();
10757ee9 1167 for (Int_t i = 0; i < nclusters; i++) {
1168 Int_t rowSorted = rowIn[index[i]];
1169 AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1170 if (!cluster) {
1171 printf("Problem\n");
1172 continue;
1173 }
1174 if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) inonEdge++;
1175 Float_t dedge = cluster->GetX()*ktany - TMath::Abs(cluster->GetY());
1176 Float_t fraction = Float_t(i) / Float_t(nclusters);
1177 Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
10757ee9 1178
1179 AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
1180
ae28e92e 1181 if (cstream) (*cstream) << "dEdx" <<
10757ee9 1182 "Cl.=" << cluster << // cluster of interest
1183 "P=" << momenta << // track momenta
1184 "dedx=" << mdedx << // mean dedx - corrected for angle
1185 "IPad=" << padType << // pad type 0..2
1186 "xc=" << xcenter << // x center of chamber
1187 "dedxQ.=" << &dedxQ << // dedxQ - total charge
1188 "dedxM.=" << &dedxM << // dedxM - maximal charge
1189 "fraction=" << fraction << // fraction - order in statistic (0,1)
1190 "fraction2=" << fraction2 << // fraction - order in statistic (0,1)
1191 "dedge=" << dedge << // distance to the edge
1192 "parY.=" << &parY << // line fit
1193 "parZ.=" << &parZ << // line fit
1194 "meanPos.=" << &meanPos << // mean position (dx, dx^2, y,y^2, z, z^2)
1195 "\n";
1196 }
0d3279d4 1197
ae28e92e 1198 if (cstream) (*cstream) << "dEdxT" <<
0d3279d4 1199 "P=" << momenta << // track momenta
1200 "npoints="<<inonEdge<< // number of points
1201 "sector="<<lastSector<< // sector number
1202 "dedx=" << mdedx << // mean dedx - corrected for angle
1203 "IPad=" << padType << // pad type 0..2
1204 "xc=" << xcenter << // x center of chamber
1205 "dedxQ.=" << &dedxQ << // dedxQ - total charge
1206 "dedxM.=" << &dedxM << // dedxM - maximal charge
1207 "parY.=" << &parY << // line fit
1208 "parZ.=" << &parZ << // line fit
1209 "meanPos.=" << &meanPos << // mean position (dx, dx^2, y,y^2, z, z^2)
1210 "\n";
1211
1212 sector = lastSector;
1213 npoints = inonEdge;
10757ee9 1214 return kTRUE;
1215}
0d3279d4 1216
1217void AliTPCcalibTracksGain::AddTracklet(UInt_t sector, UInt_t padType,TVectorD &dedxQ, TVectorD &dedxM,TVectorD& parY, TVectorD& parZ, TVectorD& meanPos){
1218 //
1219 // Add measured point - dedx to the fitter
1220 //
1221 //
8076baa0 1222 //chain->SetAlias("dr","(250-abs(meanPos.fElements[4]))/250");
0d3279d4 1223 //chain->SetAlias("tz","(0+abs(parZ.fElements[1]))");
1224 //chain->SetAlias("ty","(0+abs(parY.fElements[1]))");
1225 //chain->SetAlias("corrg","sqrt((1+ty^2)*(1+tz^2))");
8076baa0 1226 //expession fast - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000);
0d3279d4 1227
1228 Double_t xxx[100];
1229 //
1230 // z and angular part
1231 //
8076baa0 1232
1233 xxx[0] = (250.-TMath::Abs(meanPos[4]))/250.;
0d3279d4 1234 xxx[1] = TMath::Abs(parY[1]);
1235 xxx[2] = TMath::Abs(parZ[1]);
1236 xxx[3] = xxx[0]*xxx[1];
1237 xxx[4] = xxx[0]*xxx[2];
1238 xxx[5] = xxx[1]*xxx[2];
1239 xxx[6] = xxx[0]*xxx[0];
8076baa0 1240 xxx[7] = xxx[1]*xxx[1];
1241 xxx[8] = xxx[2]*xxx[2];
0d3279d4 1242 //
1243 // chamber part
1244 //
1245 Int_t tsector = sector%36;
1246 for (Int_t i=0;i<35;i++){
8076baa0 1247 xxx[9+i]=(i==tsector)?1:0;
0d3279d4 1248 }
1249 TLinearFitter *fitterM = fFitter0M;
1250 if (padType==1) fitterM=fFitter1M;
1251 if (padType==2) fitterM=fFitter2M;
1252 fitterM->AddPoint(xxx,dedxM[1]);
1253 //
1254 TLinearFitter *fitterT = fFitter0T;
1255 if (padType==1) fitterT = fFitter1T;
1256 if (padType==2) fitterT = fFitter2T;
1257 fitterT->AddPoint(xxx,dedxQ[1]);
1258}
8076baa0 1259
1260
1261TGraph *AliTPCcalibTracksGain::CreateAmpGraph(Int_t ipad, Bool_t qmax){
1262 //
1263 // create the amplitude graph
1264 // The normalized amplitudes are extrapolated to the 0 angle (y,z) and 0 drift length
1265 //
1266
1267 TVectorD vec;
1268 if (qmax){
1269 if (ipad==0) fFitter0M->GetParameters(vec);
1270 if (ipad==1) fFitter1M->GetParameters(vec);
1271 if (ipad==2) fFitter2M->GetParameters(vec);
1272 }else{
1273 if (ipad==0) fFitter0T->GetParameters(vec);
1274 if (ipad==1) fFitter1T->GetParameters(vec);
1275 if (ipad==2) fFitter2T->GetParameters(vec);
1276 }
1277
1278 Float_t amp[36];
1279 Float_t sec[36];
1280 for (Int_t i=0;i<35;i++){
1281 sec[i]=i;
1282 amp[i]=vec[10+i]+vec[0];
1283 }
1284 amp[35]=vec[0];
1285 Float_t mean = TMath::Mean(36,amp);
1286 for (Int_t i=0;i<36;i++){
1287 sec[i]=i;
1288 amp[i]=(amp[i]-mean)/mean;
1289 }
1290 TGraph *gr = new TGraph(36,sec,amp);
1291 return gr;
1292}
1293
1294
1295void AliTPCcalibTracksGain::UpdateClusterParam(AliTPCClusterParam* clparam){
1296 //
1297 // SetQ normalization parameters
1298 //
1299 // void SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm);
1300
1301 TVectorD vec;
1302 //
1303 fFitter0T->GetParameters(vec);
1304 clparam->SetQnorm(0,0,&vec);
1305 fFitter1T->GetParameters(vec);
1306 clparam->SetQnorm(1,0,&vec);
1307 fFitter2T->GetParameters(vec);
1308 clparam->SetQnorm(2,0,&vec);
1309 //
1310 fFitter0M->GetParameters(vec);
1311 clparam->SetQnorm(0,1,&vec);
1312 fFitter1M->GetParameters(vec);
1313 clparam->SetQnorm(1,1,&vec);
1314 fFitter2M->GetParameters(vec);
1315 clparam->SetQnorm(2,1,&vec);
1316 //
1317
1318}
b8601924 1319
1320
1321void AliTPCcalibTracksGain::Analyze(){
1322
1323 Evaluate();
1324
1325}
1326
1327