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