Adding stream level and debug streamer to the base class
[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 //
323 //
324 // Add profile histograms -JUST for visualization - Not used for real calibration
325 //
326 //
327 fArrayQM=new TObjArray(73); // Qmax normalized
328 fArrayQT=new TObjArray(73); // Qtot normalized
329 fProfileArrayQM = new TObjArray(37); // Qmax normalized versus local X
330 fProfileArrayQT = new TObjArray(37); // Qtot normalized versus local X
331 fProfileArrayQM2D = new TObjArray(37); // Qmax normalized versus local X and phi
332 fProfileArrayQT2D = new TObjArray(37); // Qtot normalized versus local X and phi
333 char hname[1000];
334 for (Int_t i=0; i<73; i++){
335 sprintf(hname,"QM_%d",i);
336 fArrayQM->AddAt(new TH1F(hname,hname,200,0,1000),i);
337 sprintf(hname,"QT_%d",i);
338 fArrayQT->AddAt(new TH1F(hname,hname,200,0,1000),i);
339 }
10757ee9 340
0d3279d4 341 for (Int_t i=0; i<37;i++){
342 sprintf(hname,"QMvsx_%d",i);
343 fProfileArrayQM->AddAt(new TProfile(hname,hname,50,89,250),i);
344 sprintf(hname,"QTvsx_%d",i);
345 fProfileArrayQT->AddAt(new TProfile(hname,hname,50,89,250),i);
346 sprintf(hname,"QM2D_%d",i);
347 fProfileArrayQM2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
348 sprintf(hname,"QT2D_%d",i);
349 fProfileArrayQT2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
350 }
351 //
352 // just for debugging -counters
353 //
10757ee9 354 fTotalTracks = 0;
355 fAcceptedTracks = 0;
356 fDebugCalPadRaw = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
0d3279d4 357 fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
10757ee9 358 // this will be gone for the a new ROOT version > v5-17-05
359 for (UInt_t i = 0; i < 36; i++) {
360 fNShortClusters[i] = 0;
361 fNMediumClusters[i] = 0;
362 fNLongClusters[i] = 0;
363 }
364 }
365
366AliTPCcalibTracksGain::~AliTPCcalibTracksGain() {
367 //
368 // Destructor.
369 //
370
371 Info("Destructor","");
372 if (fSimpleFitter) delete fSimpleFitter;
373 if (fSqrtFitter) delete fSqrtFitter;
374 if (fLogFitter) delete fLogFitter;
375 if (fSingleSectorFitter) delete fSingleSectorFitter;
376
10757ee9 377 if (fDebugCalPadRaw) delete fDebugCalPadRaw;
378 if (fDebugCalPadCorr) delete fDebugCalPadCorr;
379}
380
381void AliTPCcalibTracksGain::Terminate(){
382 //
383 // Evaluate fitters and close the debug stream.
384 // Also move or copy the debug stream, if a debugStreamPrefix is provided.
385 //
386
387 Evaluate();
ae28e92e 388 AliTPCcalibBase::Terminate();
10757ee9 389}
390
10757ee9 391
10757ee9 392
393void AliTPCcalibTracksGain::Process(AliTPCseed* seed) {
394 //
395 // Main method to be called when a new seed is supposed to be processed
396 // and be used for gain calibration. Its quality is checked before it
397 // is added.
398 //
399
400 fTotalTracks++;
b8601924 401 if (!fCuts->AcceptTrack(seed)) return;
10757ee9 402 fAcceptedTracks++;
403 AddTrack(seed);
404}
405
406Long64_t AliTPCcalibTracksGain::Merge(TCollection *list) {
407 //
408 // Merge() merges the results of all AliTPCcalibTracksGain objects contained in
409 // list, thus allowing a distributed computation of several files, e.g. on PROOF.
410 // The merged results are merged with the data members of the AliTPCcalibTracksGain
411 // object used for calling the Merge method.
412 // The return value is 0 /*the total number of tracks used for calibration*/ if the merge
413 // is successful, otherwise it is -1.
414 //
415
416 if (!list || list->IsEmpty()) return -1;
417
418 if (!fSimpleFitter) fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
419 if (!fSqrtFitter) fSqrtFitter = new AliTPCFitPad(8, "hyp7", "");
420 if (!fLogFitter) fLogFitter = new AliTPCFitPad(8, "hyp7", "");
421 if (!fSingleSectorFitter) fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
422
10757ee9 423
424 // just for debugging
0d3279d4 425 if (!fDebugCalPadRaw) fDebugCalPadRaw = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
426 if (!fDebugCalPadCorr) fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
10757ee9 427
428 TIterator* iter = list->MakeIterator();
429 AliTPCcalibTracksGain* cal = 0;
430
431 while ((cal = (AliTPCcalibTracksGain*)iter->Next())) {
432 if (!cal->InheritsFrom(AliTPCcalibTracksGain::Class())) {
433 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
434 return -1;
435 }
0d3279d4 436
10757ee9 437 Add(cal);
438 }
439 return 0;
440}
441
442void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
443 //
444 // Adds another AliTPCcalibTracksGain object to this object.
445 //
446
447 fSimpleFitter->Add(cal->fSimpleFitter);
448 fSqrtFitter->Add(cal->fSqrtFitter);
449 fLogFitter->Add(cal->fLogFitter);
450 fSingleSectorFitter->Add(cal->fSingleSectorFitter);
0d3279d4 451 //
452 //
453 //
454 fFitter0M->Add(cal->fFitter0M);
455 fFitter1M->Add(cal->fFitter1M);
456 fFitter2M->Add(cal->fFitter2M);
457 fFitter0T->Add(cal->fFitter0T);
458 fFitter1T->Add(cal->fFitter1T);
459 fFitter2T->Add(cal->fFitter2T);
460 //
461 //
462 // histograms
463 //
464 for (Int_t i=0; i<73; i++){
465 TH1F *his,*hism;
466 his = (TH1F*)fArrayQM->At(i);
467 hism = (TH1F*)cal->fArrayQM->At(i);
468 if (his && hism) his->Add(hism);
469 his = (TH1F*)fArrayQT->At(i);
470 hism = (TH1F*)cal->fArrayQT->At(i);
471 if (his && hism) his->Add(hism);
472 }
473 //
474 //
475 for (Int_t i=0; i<37; i++){
476 TProfile *his,*hism;
477 his = (TProfile*)fProfileArrayQM->At(i);
478 hism = (TProfile*)cal->fProfileArrayQM->At(i);
479 if (his && hism) his->Add(hism);
480 his = (TProfile*)fProfileArrayQT->At(i);
481 hism = (TProfile*)cal->fProfileArrayQT->At(i);
482 if (his && hism) his->Add(hism);
483 }
484 //
485 //
486 for (Int_t i=0; i<37; i++){
487 TProfile2D *his,*hism;
488 his = (TProfile2D*)fProfileArrayQM2D->At(i);
489 hism = (TProfile2D*)cal->fProfileArrayQM2D->At(i);
490 if (his && hism) his->Add(hism);
491 his = (TProfile2D*)fProfileArrayQT2D->At(i);
492 hism = (TProfile2D*)cal->fProfileArrayQT2D->At(i);
493 if (his && hism) his->Add(hism);
494 }
495 //
10757ee9 496 // this will be gone for the a new ROOT version > v5-17-05
497 for (UInt_t iSegment = 0; iSegment < 36; iSegment++) {
498 fNShortClusters[iSegment] += cal->fNShortClusters[iSegment];
499 fNMediumClusters[iSegment] += cal->fNMediumClusters[iSegment];
500 fNLongClusters[iSegment] += cal->fNLongClusters[iSegment];
501 }
502
503 // just for debugging, remove me
504 fTotalTracks += cal->fTotalTracks;
505 fAcceptedTracks += cal->fAcceptedTracks;
506 fDebugCalPadRaw->Add(cal->fDebugCalPadRaw);
507 fDebugCalPadCorr->Add(cal->fDebugCalPadCorr);
508
10757ee9 509}
510
511void AliTPCcalibTracksGain::AddTrack(AliTPCseed* seed) {
512 //
513 // The clusters making up the track (seed) are added to various fit functions.
514 // See AddCluster(...) for more detail.
515 //
516
0d3279d4 517 DumpTrack(seed);
518 //
519 // simple histograming part
520 for (Int_t i=0; i<159; i++){
521 AliTPCclusterMI* cluster = seed->GetClusterPointer(i);
522 if (cluster) AddCluster(cluster);
523 }
10757ee9 524}
525
0d3279d4 526void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster){
527 //
528 // Adding cluster information to the simple histograms
529 // No correction, fittings are applied
530 //
531 Float_t kThreshold=5;
532 if (cluster->GetX()<=0) return;
533 if (cluster->GetQ()<=kThreshold) return;
534 //
535
536
537 Int_t sector = cluster->GetDetector();
538 TH1F * his;
539 his = GetQT(sector);
540 if (his) his->Fill(cluster->GetQ());
541 his = GetQT(-1);
542 if (his) his->Fill(cluster->GetQ());
543 his = GetQM(sector);
544 if (his) his->Fill(cluster->GetMax());
545 his = GetQM(-1);
546 if (his) his->Fill(cluster->GetMax());
547 //
548 sector = sector%36;
549 TProfile *prof;
550 prof = GetProfileQT(sector);
551 if (prof) prof->Fill(cluster->GetX(),cluster->GetQ());
552 prof = GetProfileQT(-1);
553 if (prof) prof->Fill(cluster->GetX(),cluster->GetQ());
554 prof = GetProfileQM(sector);
555 if (prof) prof->Fill(cluster->GetX(),cluster->GetMax());
556 prof = GetProfileQM(-1);
557 if (prof) prof->Fill(cluster->GetX(),cluster->GetMax());
558 //
559 Float_t phi = cluster->GetY()/cluster->GetX();
560 TProfile2D *prof2;
561 prof2 = GetProfileQT2D(sector);
562 if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetQ());
563 prof2 = GetProfileQT2D(-1);
564 if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetQ());
565 prof2 = GetProfileQM2D(sector);
566 if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetMax());
567 prof2 = GetProfileQM2D(-1);
568 if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetMax());
569
570 //
571}
572
573
574
f1c2a4a3 575void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t /*momenta*/, Float_t/* mdedx*/, Int_t padType,
576 Float_t xcenter, TVectorD& dedxQ, TVectorD& /*dedxM*/, Float_t /*fraction*/, Float_t fraction2, Float_t dedge,
577 TVectorD& /*parY*/, TVectorD& /*parZ*/, TVectorD& meanPos) {
10757ee9 578 //
579 // Adds cluster to the appropriate fitter for later analysis.
580 // The charge used for the fit is the maximum charge for this specific cluster or the
581 // accumulated charge per cluster, depending on the value of fgkUseTotalCharge.
582 // Depending on the pad size where the cluster is registered, the value will be put in
583 // the appropriate fitter. Furthermore, for each pad size three different types of fitters
584 // are used. The fit functions are the same for all fitters (parabolic functions), but the value
585 // added to each fitter is different. The simple fitter gets the charge plugged in as is, the sqrt fitter
586 // gets the square root of the charge, and the log fitter gets fgkM*(1+q/fgkM), where q is the original charge
587 // and fgkM==25.
588 //
589
590 if (!cluster) {
591 Error("AddCluster", "Cluster not valid.");
592 return;
593 }
594
595 if (dedge < 3.) return;
596 if (fraction2 > 0.7) return;
597
598 //Int_t padType = GetPadType(cluster->GetX());
599 Double_t xx[7];
600 //Double_t centerPad[2] = {0};
601 //AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
602 //xx[0] = cluster->GetX() - centerPad[0];
603 //xx[1] = cluster->GetY() - centerPad[1];
604 xx[0] = cluster->GetX() - xcenter;
605 xx[1] = cluster->GetY();
606 xx[2] = xx[0] * xx[0];
607 xx[3] = xx[1] * xx[1];
608 xx[4] = xx[0] * xx[1];
609 xx[5] = TMath::Abs(cluster->GetZ()) - TMath::Abs(meanPos[4]);
610 xx[6] = xx[5] * xx[5];
0d3279d4 611 //
612 // Update profile histograms
613 //
10757ee9 614
0d3279d4 615 //
616 // Update fitters
617 //
10757ee9 618 Int_t segment = cluster->GetDetector() % 36;
619 Double_t q = fgkUseTotalCharge ? ((Double_t)(cluster->GetQ())) : ((Double_t)(cluster->GetMax())); // note: no normalization to pad size!
620
621 // just for debugging
622 Int_t row = 0;
623 Int_t pad = 0;
624 GetRowPad(cluster->GetX(), cluster->GetY(), row, pad);
625 fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
626
627 // correct charge by normalising to mean charge per track
628 q /= dedxQ[2];
629
630 // just for debugging
631 fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
632
633 Double_t sqrtQ = TMath::Sqrt(q);
634 Double_t logQ = fgkM * TMath::Log(1 + q / fgkM);
635 fSimpleFitter->GetFitter(segment, padType)->AddPoint(xx, q);
636 fSqrtFitter->GetFitter(segment, padType)->AddPoint(xx, sqrtQ);
637 fLogFitter->GetFitter(segment, padType)->AddPoint(xx, logQ);
638 fSingleSectorFitter->GetFitter(0, padType)->AddPoint(xx, q);
639
640 // this will be gone for the a new ROOT version > v5-17-05
641 if (padType == kShortPads)
642 fNShortClusters[segment]++;
0d3279d4 643 if (padType == kMediumPads)
10757ee9 644 fNMediumClusters[segment]++;
0d3279d4 645 if (padType == kLongPads)
10757ee9 646 fNLongClusters[segment]++;
647}
648
649void AliTPCcalibTracksGain::Evaluate(Bool_t robust, Double_t frac) {
650 //
651 // Evaluates all fitters contained in this object.
652 // If the robust option is set to kTRUE a robust fit is performed with frac as
653 // the minimal fraction of good points (see TLinearFitter::EvalRobust for details).
654 // Beware: Robust fitting is much slower!
655 //
656
657 fSimpleFitter->Evaluate(robust, frac);
658 fSqrtFitter->Evaluate(robust, frac);
659 fLogFitter->Evaluate(robust, frac);
660 fSingleSectorFitter->Evaluate(robust, frac);
0d3279d4 661 fFitter0M->Eval();
662 fFitter1M->Eval();
663 fFitter2M->Eval();
664 fFitter0T->Eval();
665 fFitter1T->Eval();
666 fFitter2T->Eval();
10757ee9 667}
668
669AliTPCCalPad* AliTPCcalibTracksGain::CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
670 //
671 // Creates the calibration object AliTPCcalPad using fitted parameterization
672 //
673 TObjArray tpc(72);
674 for (UInt_t iSector = 0; iSector < 72; iSector++)
675 tpc.Add(CreateFitCalROC(iSector, fitType, undoTransformation, normalizeToPadSize));
676 return new AliTPCCalPad(&tpc);
677}
678
679AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
680 //
681 // Create the AliTPCCalROC with the values per pad
682 // sector - sector of interest
683 // fitType -
684 //
685
686 TVectorD par(8);
687 if (sector < 36) {
688 GetParameters(sector % 36, 0, fitType, par);
689 return CreateFitCalROC(sector, 0, par, fitType, undoTransformation, normalizeToPadSize);
690 }
691 else {
692 GetParameters(sector % 36, 1, fitType, par);
693 AliTPCCalROC* roc1 = CreateFitCalROC(sector, 1, par, fitType, undoTransformation, normalizeToPadSize);
694 GetParameters(sector % 36, 2, fitType, par);
695 AliTPCCalROC* roc2 = CreateFitCalROC(sector, 2, par, fitType, undoTransformation, normalizeToPadSize);
696 AliTPCCalROC* roc3 = CreateCombinedCalROC(roc1, roc2);
697 delete roc1;
698 delete roc2;
699 return roc3;
700 }
701}
702
703AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t padType, TVectorD &fitParam, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
704 //
705 // This function is essentially a copy of AliTPCCalROC::CreateGlobalFitCalROC(...), with the
706 // modifications, that the center of the region of same pad size is used as the origin
707 // of the fit function instead of the center of the ROC.
708 // The possibility of a linear fit is removed as well because it is not needed.
709 // Only values for pads with the given pad size are calculated, the rest is 0.
710 // Set undoTransformation for undoing the transformation that was applied to the
711 // charge values before they were put into the fitter (thus allowing comparison to the original
712 // charge values). For fitType use 0 for the simple fitter, 1 for the sqrt fitter, 2 for the log fitter.
713 // If normalizeToPadSize is true, the values are normalized to the pad size.
714 // Please be aware, that you even need to specify the fitType if you want to normalize to the pad size without
715 // undoing the transformation (because normalizing involves undoing the trafo first, then normalizing, then
716 // applying the trafo again).
717 // Please note: The normalization to the pad size is a simple linear scaling with the pad length, which
718 // actually doesn't describe reality!
719 //
720
721 Float_t dlx, dly;
722 Double_t centerPad[2] = {0};
723 Float_t localXY[3] = {0};
724 AliTPCROC* tpcROC = AliTPCROC::Instance();
725 if ((padType == 0 && sector >= tpcROC->GetNInnerSector()) || (padType > 0 && sector < tpcROC->GetNInnerSector()) || sector >= tpcROC->GetNSector())
726 return 0;
727 AliTPCCalROC* lROCfitted = new AliTPCCalROC(sector);
728 //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
729 UInt_t startRow = 0;
730 UInt_t endRow = 0;
731 switch (padType) {
732 case kShortPads:
733 startRow = 0;
734 endRow = lROCfitted->GetNrows();
735 break;
736 case kMediumPads:
737 startRow = 0;
738 endRow = 64;
739 break;
740 case kLongPads:
741 startRow = 64;
742 endRow = lROCfitted->GetNrows();
743 break;
744 }
745
746 AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
747 Double_t value = 0;
748 for (UInt_t irow = startRow; irow < endRow; irow++) {
749 for (UInt_t ipad = 0; ipad < lROCfitted->GetNPads(irow); ipad++) {
750 tpcROC->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number
751 dlx = localXY[0] - centerPad[0];
752 dly = localXY[1] - centerPad[1];
753 value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
754
755 // Let q' = value be the transformed value without any pad size corrections,
756 // let T be the transformation and let l be the pad size
757 // 1) don't undo transformation, don't normalize: return q'
758 // 2) undo transformation, don't normalize: return T^{-1} q'
759 // 3) undo transformation, normalize: return (T^{-1} q') / l
760 // 4) don't undo transformation, normalize: return T((T^{-1} q') / l)
761 if (!undoTransformation && !normalizeToPadSize) {/* value remains unchanged */} // (1)
762 else { // (2), (3), (4)
763 //calculate T^{-1}
764 switch (fitType) {
765 case 0: /* value remains unchanged */ break;
766 case 1: value = value * value; break;
767 case 2: value = (TMath::Exp(value / fgkM) - 1) * fgkM; break;
768 default: Error("CreateFitCalROC", "Wrong fit type."); break;
769 }
770 if (normalizeToPadSize) value /= GetPadLength(localXY[0]); // (3)
771 }
772 if (!undoTransformation && normalizeToPadSize) { // (4)
773 // calculate T
774 switch (fitType) {
775 case 0: /* value remains unchanged */ break;
776 case 1: value = TMath::Sqrt(value); break;
777 case 2: value = fgkM * TMath::Log(1 + value / fgkM); break;
778 default: Error("CreateFitCalROC", "Wrong fit type."); break;
779 }
780 }
781 lROCfitted->SetValue(irow, ipad, value);
782 }
783 }
784 return lROCfitted;
785}
786
787AliTPCCalROC* AliTPCcalibTracksGain::CreateCombinedCalROC(const AliTPCCalROC* roc1, const AliTPCCalROC* roc2) {
788 //
789 // Combines the medium pad size values of roc1 with the long pad size values of roc2 into a new
790 // AliTPCCalROC. Returns a null pointer if any one of the ROCs is an IROC; issues a warning message
791 // if the sectors of roc1 and roc2 don't match, but still continue and use the sector of roc1 as the
792 // sector of the new ROC.
793 //
794
795 if (!roc1 || !roc2) return 0;
796 if ((Int_t)(roc1->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
797 if ((Int_t)(roc2->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
798 if (roc1->GetSector() != roc2->GetSector()) Warning("CreateCombinedCalROC", "Sector number mismatch.");
799 AliTPCCalROC* roc = new AliTPCCalROC(roc1->GetSector());
800
801 for (UInt_t iRow = 0; iRow < 64; iRow++) {
802 for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
803 roc->SetValue(iRow, iPad, roc1->GetValue(iRow, iPad));
804 }
805 for (UInt_t iRow = 64; iRow < roc->GetNrows(); iRow++) {
806 for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
807 roc->SetValue(iRow, iPad, roc2->GetValue(iRow, iPad));
808 }
809 return roc;
810}
811
812void AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) {
813 //
814 // Puts the fit parameters for the specified segment (IROC & OROC), padType and fitType
815 // into the fitParam TVectorD (which should contain 8 elements).
816 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
817 // Note: The fitter has to be evaluated first!
818 //
819
820 GetFitter(segment, padType, fitType)->GetParameters(fitParam);
821}
822
823void AliTPCcalibTracksGain::GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError) {
824 //
825 // Puts the fit parameter errors for the specified segment (IROC & OROC), padType and fitType
826 // into the fitError TVectorD (which should contain 8 elements).
827 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
828 // Note: The fitter has to be evaluated first!
829 //
830
831 GetFitter(segment, padType, fitType)->GetErrors(fitError);
832 fitError *= TMath::Sqrt(GetRedChi2(segment, padType, fitType));
833}
834
835Double_t AliTPCcalibTracksGain::GetRedChi2(UInt_t segment, UInt_t padType, UInt_t fitType) {
836 //
837 // Returns the reduced chi^2 value for the specified segment, padType and fitType.
838 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
839 // Note: The fitter has to be evaluated first!
840 //
841
842 // this will be gone for the a new ROOT version > v5-17-05
843 Int_t lNClusters = 0;
844 switch (padType) {
845 case kShortPads:
846 lNClusters = fNShortClusters[segment];
847 break;
848 case kMediumPads:
849 lNClusters = fNMediumClusters[segment];
850 break;
851 case kLongPads:
852 lNClusters = fNLongClusters[segment];
853 break;
854 }
855 return GetFitter(segment, padType, fitType)->GetChisquare()/(lNClusters - 8);
856}
857
858void AliTPCcalibTracksGain::GetCovarianceMatrix(UInt_t segment, UInt_t padType, UInt_t fitType, TMatrixD& covMatrix) {
859 //
860 // Returns the covariance matrix for the specified segment, padType, fitType.
861 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
862 //
863
864 GetFitter(segment, padType, fitType)->GetCovarianceMatrix(covMatrix);
865}
866
867TLinearFitter* AliTPCcalibTracksGain::GetFitter(UInt_t segment, UInt_t padType, UInt_t fitType) {
868 //
869 // Returns the TLinearFitter object for the specified segment, padType, fitType.
870 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
871 //
872
873 switch (fitType) {
874 case kSimpleFitter:
875 return fSimpleFitter->GetFitter(segment, padType);
876 case kSqrtFitter:
877 return fSqrtFitter->GetFitter(segment, padType);
878 case kLogFitter:
879 return fLogFitter->GetFitter(segment, padType);
880 case 3:
881 return fSingleSectorFitter->GetFitter(0, padType);
882 }
883 return 0;
884}
885
886Double_t AliTPCcalibTracksGain::GetPadLength(Double_t lx) {
887 //
888 // The function returns 0.75 for an IROC, 1. for an OROC at medium pad size position,
889 // 1.5 for an OROC at long pad size position, -1 if out of bounds.
890 //
891
892 Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
893 Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
894 Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
895 Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
896 Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
897 Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
898
899 // if IROC
900 if (lx >= irocLow && lx <= irocUp) return 0.75;
901 // if OROC medium pads
902 if (lx >= orocLow1 && lx <= orocUp1) return 1.;
903 // if OROC long pads
904 if (lx >= orocLow2 && lx <= orocUp2) return 1.5;
905 // if out of bounds
906 return -1;
907}
908
909Int_t AliTPCcalibTracksGain::GetPadType(Double_t lx) {
910 //
911 // The function returns 0 for an IROC, 1 for an OROC at medium pad size position,
912 // 2 for an OROC at long pad size position, -1 if out of bounds.
913 //
914
915 if (GetPadLength(lx) == 0.75) return 0;
916 else if (GetPadLength(lx) == 1.) return 1;
917 else if (GetPadLength(lx) == 1.5) return 2;
918 return -1;
919}
920
921// ONLY FOR DEBUGGING PURPOSES - REMOVE ME WHEN NOT NEEDED ANYMORE
922Bool_t AliTPCcalibTracksGain::GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad) {
923 //
924 // Calculate the row and pad number when the local coordinates are given.
925 // Returns kFALSE if the position is out of range, otherwise return kTRUE.
926 // WARNING: This function is preliminary and probably isn't very accurate!!
927 //
928
929 Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
930 //Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
931 Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
932 //Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
933 Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
934 //Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
935
936 if (GetPadType(lx) == 0) {
937 row = (Int_t)((lx - irocLow) / fgTPCparam->GetInnerPadPitchLength());
938 pad = (Int_t)((ly + fgTPCparam->GetYInner(row)) / fgTPCparam->GetInnerPadPitchWidth());
939 } else if (GetPadType(lx) == 1) {
940 row = (Int_t)((lx - orocLow1) / fgTPCparam->GetOuter1PadPitchLength());
941 pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
942 } else if (GetPadType(lx) == 2) {
943 row = fgTPCparam->GetNRowUp1() + (Int_t)((lx - orocLow2) / fgTPCparam->GetOuter2PadPitchLength());
944 pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
945 }
946 else return kFALSE;
947 return kTRUE;
948}
949
950void AliTPCcalibTracksGain::DumpTrack(AliTPCseed* track) {
951 //
952 // Dump track information to the debug stream
953 //
954
10757ee9 955 Int_t rows[200];
0d3279d4 956 Int_t sector[3];
957 Int_t npoints[3];
958 TVectorD dedxM[3];
959 TVectorD dedxQ[3];
960 TVectorD parY[3];
961 TVectorD parZ[3];
962 TVectorD meanPos[3];
963 //
964 Int_t count=0;
10757ee9 965 for (Int_t ipad = 0; ipad < 3; ipad++) {
0d3279d4 966 dedxM[ipad].ResizeTo(5);
967 dedxQ[ipad].ResizeTo(5);
968 parY[ipad].ResizeTo(3);
969 parZ[ipad].ResizeTo(3);
970 meanPos[ipad].ResizeTo(6);
971 Bool_t isOK = GetDedx(track, ipad, rows, sector[ipad], npoints[ipad], dedxM[ipad], dedxQ[ipad], parY[ipad], parZ[ipad], meanPos[ipad]);
972 if (isOK)
973 AddTracklet(sector[ipad],ipad, dedxQ[ipad], dedxM[ipad], parY[ipad], parZ[ipad], meanPos[ipad] );
974 if (isOK) count++;
10757ee9 975 }
ae28e92e 976
977 TTreeSRedirector * cstream = GetDebugStreamer();
978 if (cstream){
979 (*cstream) << "Track" <<
0d3279d4 980 "Track.=" << track << // track information
0d3279d4 981 "\n";
ae28e92e 982 //
983 //
984 //
985 if ( GetStreamLevel()>1 && count>1){
986 (*cstream) << "TrackG" <<
987 "Track.=" << track << // track information
988 "ncount="<<count<<
989 // info for pad type 0
990 "sector0="<<sector[0]<<
991 "npoints0="<<npoints[0]<<
992 "dedxM0.="<<&dedxM[0]<<
993 "dedxQ0.="<<&dedxQ[0]<<
994 "parY0.="<<&parY[0]<<
995 "parZ0.="<<&parZ[0]<<
996 "meanPos0.="<<&meanPos[0]<<
997 //
998 // info for pad type 1
999 "sector1="<<sector[1]<<
1000 "npoints1="<<npoints[1]<<
1001 "dedxM1.="<<&dedxM[1]<<
1002 "dedxQ1.="<<&dedxQ[1]<<
1003 "parY1.="<<&parY[1]<<
1004 "parZ1.="<<&parZ[1]<<
1005 "meanPos1.="<<&meanPos[1]<<
1006 //
1007 // info for pad type 2
1008 "sector2="<<sector[2]<<
1009 "npoints2="<<npoints[2]<<
1010 "dedxM2.="<<&dedxM[2]<<
1011 "dedxQ2.="<<&dedxQ[2]<<
1012 "parY2.="<<&parY[2]<<
1013 "parZ2.="<<&parZ[2]<<
1014 "meanPos2.="<<&meanPos[2]<<
1015 //
1016 "\n";
1017 }
0d3279d4 1018 }
0d3279d4 1019}
1020
f1c2a4a3 1021Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /*rows*/,
0d3279d4 1022 Int_t &sector, Int_t& npoints,
1023 TVectorD &dedxM, TVectorD &dedxQ,
1024 TVectorD &parY, TVectorD &parZ, TVectorD&meanPos)
1025{
1026 //
1027 // GetDedx for given sector for given track
1028 // padType - type of pads
1029 //
1030
1031 static TLinearFitter fitY(2, "pol1");
1032 static TLinearFitter fitZ(2, "pol1");
1033 //
1034 //
10757ee9 1035 Int_t firstRow = 0, lastRow = 0;
1036 Int_t minRow = 100;
1037 Float_t xcenter = 0;
1038 const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
1039 const Float_t kedgey = 4.;
1040 if (padType == 0) {
1041 firstRow = 0;
1042 lastRow = fgTPCparam->GetNRowLow();
1043 xcenter = 108.47;
1044 }
1045 if (padType == 1) {
1046 firstRow = fgTPCparam->GetNRowLow();
1047 lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1048 xcenter = 166.60;
1049 }
1050 if (padType == 2) {
1051 firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1052 lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
1053 xcenter = 222.6;
1054 }
1055 minRow = (lastRow - firstRow) / 2;
1056 //
1057 //
1058 Int_t nclusters = 0;
1059 Int_t nclustersNE = 0; // number of not edge clusters
1060 Int_t lastSector = -1;
1061 Float_t amplitudeQ[100];
1062 Float_t amplitudeM[100];
1063 Int_t rowIn[100];
1064 Int_t index[100];
1065 //
0d3279d4 1066 //
10757ee9 1067 fitY.ClearPoints();
1068 fitZ.ClearPoints();
10757ee9 1069
1070 for (Int_t iCluster = firstRow; iCluster < lastRow; iCluster++) {
1071 AliTPCclusterMI* cluster = track->GetClusterPointer(iCluster);
1072 if (cluster) {
1073 Int_t detector = cluster->GetDetector() ;
1074 if (lastSector == -1) lastSector = detector;
1075 if (lastSector != detector) continue;
1076 amplitudeQ[nclusters] = cluster->GetQ();
1077 amplitudeM[nclusters] = cluster->GetMax();
1078 rowIn[nclusters] = iCluster;
1079 nclusters++;
1080 Double_t dx = cluster->GetX() - xcenter;
1081 Double_t y = cluster->GetY();
1082 Double_t z = cluster->GetZ();
1083 fitY.AddPoint(&dx, y);
1084 fitZ.AddPoint(&dx, z);
1085 meanPos[0] += dx;
1086 meanPos[1] += dx;
1087 meanPos[2] += y;
1088 meanPos[3] += y*y;
1089 meanPos[4] += z;
1090 meanPos[5] += z*z;
1091 if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) nclustersNE++;
1092 }
1093 }
1094
1095 if (nclusters < minRow / 2) return kFALSE;
1096 if (nclustersNE < minRow / 2) return kFALSE;
1097 for (Int_t i = 0; i < 6; i++) meanPos[i] /= Double_t(nclusters);
1098 fitY.Eval();
1099 fitZ.Eval();
1100 fitY.GetParameters(parY);
1101 fitZ.GetParameters(parZ);
1102 //
1103 // calculate truncated mean
1104 //
1105 TMath::Sort(nclusters, amplitudeQ, index, kFALSE);
0d3279d4 1106 //
1107 //
1108 //
10757ee9 1109 Float_t ndedx[5];
1110 for (Int_t i = 0; i < 5; i++) {
1111 dedxQ[i] = 0;
1112 dedxM[i] = 0;
1113 ndedx[i] = 0;
1114 }
1115 //
1116 // dedx calculation
1117 //
1118 Int_t inonEdge = 0;
1119 for (Int_t i = 0; i < nclusters; i++) {
1120 Int_t rowSorted = rowIn[index[i]];
1121 AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1122
1123 if (TMath::Abs(cluster->GetY()) > cluster->GetX()*ktany - kedgey) continue; //don't take edge clusters
1124 inonEdge++;
1125 if (inonEdge < nclustersNE * 0.5) {
1126 ndedx[0]++;
1127 dedxQ[0] += amplitudeQ[index[i]];
1128 dedxM[0] += amplitudeM[index[i]];
1129 }
1130 if (inonEdge < nclustersNE * 0.6) {
1131 ndedx[1]++;
1132 dedxQ[1] += amplitudeQ[index[i]];
1133 dedxM[1] += amplitudeM[index[i]];
1134 }
1135 if (inonEdge < nclustersNE * 0.7) {
1136 ndedx[2]++;
1137 dedxQ[2] += amplitudeQ[index[i]];
1138 dedxM[2] += amplitudeM[index[i]];
1139 }
1140 if (inonEdge < nclustersNE * 0.8) {
1141 ndedx[3]++;
1142 dedxQ[3] += amplitudeQ[index[i]];
1143 dedxM[3] += amplitudeM[index[i]];
1144 }
1145 if (inonEdge < nclustersNE * 0.9) {
1146 ndedx[4]++;
1147 dedxQ[4] += amplitudeQ[index[i]];
1148 dedxM[4] += amplitudeM[index[i]];
1149 }
1150 }
1151 for (Int_t i = 0; i < 5; i++) {
1152 dedxQ[i] /= ndedx[i];
1153 dedxM[i] /= ndedx[i];
1154 }
ae28e92e 1155 TTreeSRedirector * cstream = GetDebugStreamer();
10757ee9 1156 inonEdge = 0;
0d3279d4 1157 Float_t momenta = track->GetP();
1158 Float_t mdedx = track->GetdEdx();
10757ee9 1159 for (Int_t i = 0; i < nclusters; i++) {
1160 Int_t rowSorted = rowIn[index[i]];
1161 AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1162 if (!cluster) {
1163 printf("Problem\n");
1164 continue;
1165 }
1166 if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) inonEdge++;
1167 Float_t dedge = cluster->GetX()*ktany - TMath::Abs(cluster->GetY());
1168 Float_t fraction = Float_t(i) / Float_t(nclusters);
1169 Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
10757ee9 1170
1171 AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
1172
ae28e92e 1173 if (cstream) (*cstream) << "dEdx" <<
10757ee9 1174 "Cl.=" << cluster << // cluster of interest
1175 "P=" << momenta << // track momenta
1176 "dedx=" << mdedx << // mean dedx - corrected for angle
1177 "IPad=" << padType << // pad type 0..2
1178 "xc=" << xcenter << // x center of chamber
1179 "dedxQ.=" << &dedxQ << // dedxQ - total charge
1180 "dedxM.=" << &dedxM << // dedxM - maximal charge
1181 "fraction=" << fraction << // fraction - order in statistic (0,1)
1182 "fraction2=" << fraction2 << // fraction - order in statistic (0,1)
1183 "dedge=" << dedge << // distance to the edge
1184 "parY.=" << &parY << // line fit
1185 "parZ.=" << &parZ << // line fit
1186 "meanPos.=" << &meanPos << // mean position (dx, dx^2, y,y^2, z, z^2)
1187 "\n";
1188 }
0d3279d4 1189
ae28e92e 1190 if (cstream) (*cstream) << "dEdxT" <<
0d3279d4 1191 "P=" << momenta << // track momenta
1192 "npoints="<<inonEdge<< // number of points
1193 "sector="<<lastSector<< // sector number
1194 "dedx=" << mdedx << // mean dedx - corrected for angle
1195 "IPad=" << padType << // pad type 0..2
1196 "xc=" << xcenter << // x center of chamber
1197 "dedxQ.=" << &dedxQ << // dedxQ - total charge
1198 "dedxM.=" << &dedxM << // dedxM - maximal charge
1199 "parY.=" << &parY << // line fit
1200 "parZ.=" << &parZ << // line fit
1201 "meanPos.=" << &meanPos << // mean position (dx, dx^2, y,y^2, z, z^2)
1202 "\n";
1203
1204 sector = lastSector;
1205 npoints = inonEdge;
10757ee9 1206 return kTRUE;
1207}
0d3279d4 1208
1209void AliTPCcalibTracksGain::AddTracklet(UInt_t sector, UInt_t padType,TVectorD &dedxQ, TVectorD &dedxM,TVectorD& parY, TVectorD& parZ, TVectorD& meanPos){
1210 //
1211 // Add measured point - dedx to the fitter
1212 //
1213 //
8076baa0 1214 //chain->SetAlias("dr","(250-abs(meanPos.fElements[4]))/250");
0d3279d4 1215 //chain->SetAlias("tz","(0+abs(parZ.fElements[1]))");
1216 //chain->SetAlias("ty","(0+abs(parY.fElements[1]))");
1217 //chain->SetAlias("corrg","sqrt((1+ty^2)*(1+tz^2))");
8076baa0 1218 //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 1219
1220 Double_t xxx[100];
1221 //
1222 // z and angular part
1223 //
8076baa0 1224
1225 xxx[0] = (250.-TMath::Abs(meanPos[4]))/250.;
0d3279d4 1226 xxx[1] = TMath::Abs(parY[1]);
1227 xxx[2] = TMath::Abs(parZ[1]);
1228 xxx[3] = xxx[0]*xxx[1];
1229 xxx[4] = xxx[0]*xxx[2];
1230 xxx[5] = xxx[1]*xxx[2];
1231 xxx[6] = xxx[0]*xxx[0];
8076baa0 1232 xxx[7] = xxx[1]*xxx[1];
1233 xxx[8] = xxx[2]*xxx[2];
0d3279d4 1234 //
1235 // chamber part
1236 //
1237 Int_t tsector = sector%36;
1238 for (Int_t i=0;i<35;i++){
8076baa0 1239 xxx[9+i]=(i==tsector)?1:0;
0d3279d4 1240 }
1241 TLinearFitter *fitterM = fFitter0M;
1242 if (padType==1) fitterM=fFitter1M;
1243 if (padType==2) fitterM=fFitter2M;
1244 fitterM->AddPoint(xxx,dedxM[1]);
1245 //
1246 TLinearFitter *fitterT = fFitter0T;
1247 if (padType==1) fitterT = fFitter1T;
1248 if (padType==2) fitterT = fFitter2T;
1249 fitterT->AddPoint(xxx,dedxQ[1]);
1250}
8076baa0 1251
1252
1253TGraph *AliTPCcalibTracksGain::CreateAmpGraph(Int_t ipad, Bool_t qmax){
1254 //
1255 // create the amplitude graph
1256 // The normalized amplitudes are extrapolated to the 0 angle (y,z) and 0 drift length
1257 //
1258
1259 TVectorD vec;
1260 if (qmax){
1261 if (ipad==0) fFitter0M->GetParameters(vec);
1262 if (ipad==1) fFitter1M->GetParameters(vec);
1263 if (ipad==2) fFitter2M->GetParameters(vec);
1264 }else{
1265 if (ipad==0) fFitter0T->GetParameters(vec);
1266 if (ipad==1) fFitter1T->GetParameters(vec);
1267 if (ipad==2) fFitter2T->GetParameters(vec);
1268 }
1269
1270 Float_t amp[36];
1271 Float_t sec[36];
1272 for (Int_t i=0;i<35;i++){
1273 sec[i]=i;
1274 amp[i]=vec[10+i]+vec[0];
1275 }
1276 amp[35]=vec[0];
1277 Float_t mean = TMath::Mean(36,amp);
1278 for (Int_t i=0;i<36;i++){
1279 sec[i]=i;
1280 amp[i]=(amp[i]-mean)/mean;
1281 }
1282 TGraph *gr = new TGraph(36,sec,amp);
1283 return gr;
1284}
1285
1286
1287void AliTPCcalibTracksGain::UpdateClusterParam(AliTPCClusterParam* clparam){
1288 //
1289 // SetQ normalization parameters
1290 //
1291 // void SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm);
1292
1293 TVectorD vec;
1294 //
1295 fFitter0T->GetParameters(vec);
1296 clparam->SetQnorm(0,0,&vec);
1297 fFitter1T->GetParameters(vec);
1298 clparam->SetQnorm(1,0,&vec);
1299 fFitter2T->GetParameters(vec);
1300 clparam->SetQnorm(2,0,&vec);
1301 //
1302 fFitter0M->GetParameters(vec);
1303 clparam->SetQnorm(0,1,&vec);
1304 fFitter1M->GetParameters(vec);
1305 clparam->SetQnorm(1,1,&vec);
1306 fFitter2M->GetParameters(vec);
1307 clparam->SetQnorm(2,1,&vec);
1308 //
1309
1310}
b8601924 1311
1312
1313void AliTPCcalibTracksGain::Analyze(){
1314
1315 Evaluate();
1316
1317}
1318
1319