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