]> 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
Additional protection (Yves)
[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//
684602c8 99//
100// In order to debug some numerical algorithm all data data which are used for
101// fitters can be stored in the debug streamers. In case of fitting roblems the
102// errors and tendencies can be checked.
103//
104// Debug Streamers:
105//
106//
107//
108//
109//
10757ee9 110////////////////////////////////////////////////////////////////////////////
111
684602c8 112/*
113 .x ~/UliStyle.C
114 gSystem->Load("libANALYSIS");
da7d274e 115 gSystem->Load("libSTAT");
684602c8 116 gSystem->Load("libTPCcalib");
da7d274e 117
684602c8 118 TFile fcalib("CalibObjects.root");
119 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
120 AliTPCcalibTracksGain * gain = ( AliTPCcalibTracksGain *)array->FindObject("calibTracksGain");
121
2acad464 122 //
123 // Angular and drift correction
124 //
125 AliTPCClusterParam *param = new AliTPCClusterParam;param->SetInstance(param);
126 gain->UpdateClusterParam(param);
127 TF2 fdrty("fdrty","AliTPCClusterParam::SQnorm(0,0,x,y,0)",0,1,0,1)
128
684602c8 129 //
130 // Make visual Tree - compare with Kr calibration
131 //
132 AliTPCCalPad * gain010 = gain->CreateFitCalPad(0,kTRUE,0); gain010->SetName("CGain010");
133 AliTPCCalPad * gain110 = gain->CreateFitCalPad(1,kTRUE,0); gain110->SetName("CGain110");
134 AliTPCCalPad * gain210 = gain->CreateFitCalPad(2,kTRUE,0); gain210->SetName("CGain210");
2acad464 135 TFile fkr("/u/miranov/GainMap.root");
684602c8 136 AliTPCCalPad *gainKr = fkr.Get("GainMap"); fkr->SetName("KrGain");
137 //
138 AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
139 preprocesor->AddComponent(gain010);
140 preprocesor->AddComponent(gain110);
141 preprocesor->AddComponent(gain210);
142 preprocesor->AddComponent(gainKr);
143 preprocesor->DumpToFile("cosmicGain.root");
144 //
145 //
146 //
147 // Simple session using the debug streamers
148 //
149
150 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
151 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
152 AliXRDPROOFtoolkit tool;
153
154 TChain * chain0 = tool.MakeChain("gain.txt","dEdx",0,1000000);
155 TChain * chain1 = tool.MakeChain("gain.txt","Track",0,1000000);
156 TChain * chain2 = tool.MakeChain("gain.txt","TrackG",0,1000000);
157 chain0->Lookup();
158 chain1->Lookup();
159 chain2->Lookup();
160
161 chain2->SetAlias("k1","1/0.855");
162 chain2->SetAlias("k0","1/0.9928");
163 chain2->SetAlias("k2","1/1.152");
164
684602c8 165*/
166
167
168
169
170
10757ee9 171#include <TPDGCode.h>
172#include <TStyle.h>
173#include "TSystem.h"
174#include "TMatrixD.h"
175#include "TTreeStream.h"
176#include "TF1.h"
177#include "AliTPCParamSR.h"
178#include "AliTPCClusterParam.h"
179#include "AliTrackPointArray.h"
180#include "TCint.h"
181#include "AliTPCcalibTracksGain.h"
182#include <TH1.h>
183#include <TH3F.h>
184#include <TLinearFitter.h>
185#include <TTreeStream.h>
186#include <TFile.h>
187#include <TCollection.h>
188#include <TIterator.h>
5b00528f 189#include <TProfile.h>
190#include <TProfile2D.h>
c1418a4c 191#include <TProof.h>
10757ee9 192
193//
194// AliRoot includes
195//
196#include "AliMagF.h"
197#include "AliMathBase.h"
198//
199#include "AliTPCROC.h"
200#include "AliTPCParamSR.h"
201#include "AliTPCCalROC.h"
202#include "AliTPCCalPad.h"
8076baa0 203#include "AliTPCClusterParam.h"
10757ee9 204//
205#include "AliTracker.h"
206#include "AliESD.h"
207#include "AliESDtrack.h"
208#include "AliESDfriend.h"
209#include "AliESDfriendTrack.h"
210#include "AliTPCseed.h"
211#include "AliTPCclusterMI.h"
212#include "AliTPCcalibTracksCuts.h"
213#include "AliTPCFitPad.h"
da7d274e 214#include "TStatToolkit.h"
215#include "TString.h"
216#include "TCut.h"
10757ee9 217
da7d274e 218//
10757ee9 219#include <TTree.h>
220#include "AliESDEvent.h"
221
222/*
223
224TFile f("TPCCalibTracksGain.root")
225
226gSystem->Load("libPWG1.so")
227AliTreeDraw comp
228comp.SetTree(dEdx)
229Double_t chi2;
230TVectorD vec(3)
231TMatrixD mat(3,3)
232TString * str = comp.FitPlane("Cl.fQ/dedxQ.fElements[0]","Cl.fY++Cl.fX","Cl.fDetector<36",chi2,vec,mat)
233
234*/
235
236ClassImp(AliTPCcalibTracksGain)
237
238const Bool_t AliTPCcalibTracksGain::fgkUseTotalCharge = kTRUE;
239const Double_t AliTPCcalibTracksGain::fgkM = 25.;
240const char* AliTPCcalibTracksGain::fgkDebugStreamFileName = "TPCCalibTracksGain.root";
241AliTPCParamSR* AliTPCcalibTracksGain::fgTPCparam = new AliTPCParamSR();
242
243AliTPCcalibTracksGain::AliTPCcalibTracksGain() :
b8601924 244 AliTPCcalibBase(),
0d3279d4 245 fCuts(0), // cuts that are used for sieving the tracks used for calibration
2acad464 246 fGainMap(0), // gain map to be applied
0d3279d4 247 //
248 // Simple Array of histograms
249 //
250 fArrayQM(0), // Qmax normalized
251 fArrayQT(0), // Qtot normalized
252 fProfileArrayQM(0), // Qmax normalized versus local X
253 fProfileArrayQT(0), // Qtot normalized versus local X
254 fProfileArrayQM2D(0), // Qmax normalized versus local X and phi
255 fProfileArrayQT2D(0), // Qtot normalized versus local X and phi
256 //
257 // Fitters
258 //
259 fSimpleFitter(0), // simple fitter for short pads
260 fSqrtFitter(0), // sqrt fitter for medium pads
261 fLogFitter(0), // log fitter for long pads
684602c8 262
0d3279d4 263 fFitter0M(0), // fitting of the atenuation, angular correction, and mean chamber gain
264 fFitter1M(0), // fitting of the atenuation, angular correction, and mean chamber gain
265 fFitter2M(0), // fitting of the atenuation, angular correction, and mean chamber gain
266 fFitter0T(0), // fitting of the atenuation, angular correction, and mean chamber gain
267 fFitter1T(0), // fitting of the atenuation, angular correction, and mean chamber gain
268 fFitter2T(0), // fitting of the atenuation, angular correction, and mean chamber gain
684602c8 269 //
270 fDFitter0M(0), // fitting of the atenuation, angular correction
271 fDFitter1M(0), // fitting of the atenuation, angular correction
272 fDFitter2M(0), // fitting of the atenuation, angular correction
273 fDFitter0T(0), // fitting of the atenuation, angular correction
274 fDFitter1T(0), // fitting of the atenuation, angular correction
275 fDFitter2T(0), // fitting of the atenuation, angular correction
276 //
0d3279d4 277 fSingleSectorFitter(0), // just for debugging
278 //
279 // Counters
280 //
281 fTotalTracks(0), // just for debugging
282 fAcceptedTracks(0), // just for debugging
283 fDebugCalPadRaw(0), // just for debugging
2acad464 284 fDebugCalPadCorr(0) // just for debugging
0d3279d4 285
10757ee9 286{
287 //
288 // Default constructor.
289 //
290}
291
292AliTPCcalibTracksGain::AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj) :
b8601924 293 AliTPCcalibBase(obj),
0d3279d4 294 fCuts(obj.fCuts), // cuts that are used for sieving the tracks used for calibration
2acad464 295 fGainMap(new AliTPCCalPad(*(obj.fGainMap))), // gain map to be applied
f1c2a4a3 296 fArrayQM(0), // Qmax normalized
297 fArrayQT(0), // Qtot normalized
0d3279d4 298 //
299 // Simple Histograms
300 //
301 fProfileArrayQM(obj.fProfileArrayQM), // Qmax normalized versus local X
302 fProfileArrayQT(obj.fProfileArrayQT), // Qtot normalized versus local X
303 fProfileArrayQM2D(obj.fProfileArrayQM2D), // Qmax normalized versus local X and phi
304 fProfileArrayQT2D(obj.fProfileArrayQT2D), // Qtot normalized versus local X and phi
305 //
306 // Fitters
307 //
308 fSimpleFitter(obj.fSimpleFitter), // simple fitter for short pads
309 fSqrtFitter(obj.fSqrtFitter), // sqrt fitter for medium pads
310 fLogFitter(obj.fLogFitter), // log fitter for long pads
311 fFitter0M(obj.fFitter0M),
312 fFitter1M(obj.fFitter1M),
313 fFitter2M(obj.fFitter2M),
314 fFitter0T(obj.fFitter0T),
315 fFitter1T(obj.fFitter1T),
316 fFitter2T(obj.fFitter2T),
684602c8 317 //
318 fDFitter0M(obj.fDFitter0M),
319 fDFitter1M(obj.fDFitter1M),
320 fDFitter2M(obj.fDFitter2M),
321 fDFitter0T(obj.fDFitter0T),
322 fDFitter1T(obj.fDFitter1T),
323 fDFitter2T(obj.fDFitter2T),
0d3279d4 324 fSingleSectorFitter(obj.fSingleSectorFitter), // just for debugging
325 //
326 // Counters
327 //
328 fTotalTracks(obj.fTotalTracks), // just for debugging
329 fAcceptedTracks(obj.fAcceptedTracks), // just for debugging
330 fDebugCalPadRaw(obj.fDebugCalPadRaw), // just for debugging
2acad464 331 fDebugCalPadCorr(obj.fDebugCalPadCorr) // just for debugging
0d3279d4 332
10757ee9 333{
334 //
335 // Copy constructor.
336 //
10757ee9 337}
338
339AliTPCcalibTracksGain& AliTPCcalibTracksGain::operator=(const AliTPCcalibTracksGain& rhs) {
340 //
341 // Assignment operator.
342 //
343
344 if (this != &rhs) {
345 TNamed::operator=(rhs);
346 fDebugCalPadRaw = new AliTPCCalPad(*(rhs.fDebugCalPadRaw));
347 fDebugCalPadCorr = new AliTPCCalPad(*(rhs.fDebugCalPadCorr));
348 fSimpleFitter = new AliTPCFitPad(*(rhs.fSimpleFitter));
349 fSqrtFitter = new AliTPCFitPad(*(rhs.fSqrtFitter));
350 fLogFitter = new AliTPCFitPad(*(rhs.fLogFitter));
351 fSingleSectorFitter = new AliTPCFitPad(*(rhs.fSingleSectorFitter));
10757ee9 352 fCuts = new AliTPCcalibTracksCuts(*(rhs.fCuts));
2acad464 353 fGainMap = new AliTPCCalPad(*(rhs.fGainMap));
10757ee9 354 }
355 return *this;
356}
357
f1c2a4a3 358AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts, TNamed* /*debugStreamPrefix*/, AliTPCcalibTracksGain* prevIter) :
b8601924 359 AliTPCcalibBase(),
0d3279d4 360 fCuts(0), // cuts that are used for sieving the tracks used for calibration
2acad464 361 fGainMap(0), // gain map to be applied
f1c2a4a3 362 fArrayQM(0), // Qmax normalized
363 fArrayQT(0), // Qtot normalized
0d3279d4 364 //
365 // Simple Histograms
366 //
367 fProfileArrayQM(0), // Qmax normalized versus local X
368 fProfileArrayQT(0), // Qtot normalized versus local X
369 fProfileArrayQM2D(0), // Qmax normalized versus local X and phi
370 fProfileArrayQT2D(0), // Qtot normalized versus local X and phi
371 //
372 // Fitters
373 //
374 fSimpleFitter(0), // simple fitter for short pads
375 fSqrtFitter(0), // sqrt fitter for medium pads
376 fLogFitter(0), // log fitter for long pads
377 fFitter0M(0), // fitting of the atenuation, angular correction, and mean chamber gain
378 fFitter1M(0), // fitting of the atenuation, angular correction, and mean chamber gain
379 fFitter2M(0), // fitting of the atenuation, angular correction, and mean chamber gain
380 fFitter0T(0), // fitting of the atenuation, angular correction, and mean chamber gain
381 fFitter1T(0), // fitting of the atenuation, angular correction, and mean chamber gain
382 fFitter2T(0), // fitting of the atenuation, angular correction, and mean chamber gain
684602c8 383 //
384 fDFitter0M(0), // fitting of the atenuation, angular correction
385 fDFitter1M(0), // fitting of the atenuation, angular correction
386 fDFitter2M(0), // fitting of the atenuation, angular correction
387 fDFitter0T(0), // fitting of the atenuation, angular correction
388 fDFitter1T(0), // fitting of the atenuation, angular correction
389 fDFitter2T(0), // fitting of the atenuation, angular correction
0d3279d4 390 fSingleSectorFitter(0), // just for debugging
391 //
392 // Counters
393 //
394 fTotalTracks(0), // just for debugging
395 fAcceptedTracks(0), // just for debugging
396 fDebugCalPadRaw(0), // just for debugging
2acad464 397 fDebugCalPadCorr(0) // just for debugging
0d3279d4 398
10757ee9 399{
400 //
401 // Constructor.
0d3279d4 402 //
b8601924 403 this->SetNameTitle(name, title);
10757ee9 404 fCuts = cuts;
0d3279d4 405 //
406 // Fitter initialization
407 //
10757ee9 408 fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
409 fSqrtFitter = new AliTPCFitPad(8, "hyp7", "");
410 fLogFitter = new AliTPCFitPad(8, "hyp7", "");
411 fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
8076baa0 412 //
413 fFitter0M = new TLinearFitter(45,"hyp44");
414 fFitter1M = new TLinearFitter(45,"hyp44");
415 fFitter2M = new TLinearFitter(45,"hyp44");
416 fFitter0T = new TLinearFitter(45,"hyp44");
417 fFitter1T = new TLinearFitter(45,"hyp44");
418 fFitter2T = new TLinearFitter(45,"hyp44");
0d3279d4 419 //
684602c8 420 fDFitter0M = new TLinearFitter(10,"hyp9");
421 fDFitter1M = new TLinearFitter(10,"hyp9");
422 fDFitter2M = new TLinearFitter(10,"hyp9");
423 fDFitter0T = new TLinearFitter(10,"hyp9");
424 fDFitter1T = new TLinearFitter(10,"hyp9");
425 fDFitter2T = new TLinearFitter(10,"hyp9");
426 //
cbc19295 427 //
428 fFitter0M->StoreData(kFALSE);
429 fFitter1M->StoreData(kFALSE);
430 fFitter2M->StoreData(kFALSE);
431 fFitter0T->StoreData(kFALSE);
432 fFitter1T->StoreData(kFALSE);
433 fFitter2T->StoreData(kFALSE);
434 //
684602c8 435 fDFitter0M->StoreData(kFALSE);
436 fDFitter1M->StoreData(kFALSE);
437 fDFitter2M->StoreData(kFALSE);
438 fDFitter0T->StoreData(kFALSE);
439 fDFitter1T->StoreData(kFALSE);
440 fDFitter2T->StoreData(kFALSE);
441 //
0d3279d4 442 //
443 // Add profile histograms -JUST for visualization - Not used for real calibration
444 //
445 //
446 fArrayQM=new TObjArray(73); // Qmax normalized
447 fArrayQT=new TObjArray(73); // Qtot normalized
448 fProfileArrayQM = new TObjArray(37); // Qmax normalized versus local X
449 fProfileArrayQT = new TObjArray(37); // Qtot normalized versus local X
450 fProfileArrayQM2D = new TObjArray(37); // Qmax normalized versus local X and phi
451 fProfileArrayQT2D = new TObjArray(37); // Qtot normalized versus local X and phi
452 char hname[1000];
453 for (Int_t i=0; i<73; i++){
454 sprintf(hname,"QM_%d",i);
455 fArrayQM->AddAt(new TH1F(hname,hname,200,0,1000),i);
456 sprintf(hname,"QT_%d",i);
457 fArrayQT->AddAt(new TH1F(hname,hname,200,0,1000),i);
458 }
10757ee9 459
0d3279d4 460 for (Int_t i=0; i<37;i++){
461 sprintf(hname,"QMvsx_%d",i);
462 fProfileArrayQM->AddAt(new TProfile(hname,hname,50,89,250),i);
463 sprintf(hname,"QTvsx_%d",i);
464 fProfileArrayQT->AddAt(new TProfile(hname,hname,50,89,250),i);
465 sprintf(hname,"QM2D_%d",i);
466 fProfileArrayQM2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
467 sprintf(hname,"QT2D_%d",i);
468 fProfileArrayQT2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
469 }
470 //
471 // just for debugging -counters
472 //
10757ee9 473 fTotalTracks = 0;
474 fAcceptedTracks = 0;
475 fDebugCalPadRaw = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
0d3279d4 476 fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
10757ee9 477 // this will be gone for the a new ROOT version > v5-17-05
478 for (UInt_t i = 0; i < 36; i++) {
479 fNShortClusters[i] = 0;
480 fNMediumClusters[i] = 0;
481 fNLongClusters[i] = 0;
482 }
483 }
484
485AliTPCcalibTracksGain::~AliTPCcalibTracksGain() {
486 //
487 // Destructor.
488 //
489
490 Info("Destructor","");
491 if (fSimpleFitter) delete fSimpleFitter;
492 if (fSqrtFitter) delete fSqrtFitter;
493 if (fLogFitter) delete fLogFitter;
494 if (fSingleSectorFitter) delete fSingleSectorFitter;
495
10757ee9 496 if (fDebugCalPadRaw) delete fDebugCalPadRaw;
497 if (fDebugCalPadCorr) delete fDebugCalPadCorr;
498}
499
500void AliTPCcalibTracksGain::Terminate(){
501 //
502 // Evaluate fitters and close the debug stream.
503 // Also move or copy the debug stream, if a debugStreamPrefix is provided.
504 //
505
506 Evaluate();
ae28e92e 507 AliTPCcalibBase::Terminate();
10757ee9 508}
509
10757ee9 510
10757ee9 511
512void AliTPCcalibTracksGain::Process(AliTPCseed* seed) {
513 //
514 // Main method to be called when a new seed is supposed to be processed
515 // and be used for gain calibration. Its quality is checked before it
516 // is added.
517 //
518
c1418a4c 519
10757ee9 520 fTotalTracks++;
b8601924 521 if (!fCuts->AcceptTrack(seed)) return;
c1418a4c 522 //
523 // reinint on proof
524 // if (gProof){
525 static Bool_t doinit= kTRUE;
526 if (doinit){
527 fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
528 fSqrtFitter = new AliTPCFitPad(8, "hyp7", "");
529 fLogFitter = new AliTPCFitPad(8, "hyp7", "");
530 fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
531 //
532 fFitter0M = new TLinearFitter(45,"hyp44");
533 fFitter1M = new TLinearFitter(45,"hyp44");
534 fFitter2M = new TLinearFitter(45,"hyp44");
535 fFitter0T = new TLinearFitter(45,"hyp44");
536 fFitter1T = new TLinearFitter(45,"hyp44");
537 fFitter2T = new TLinearFitter(45,"hyp44");
684602c8 538 //
539 fDFitter0M = new TLinearFitter(10,"hyp9");
540 fDFitter1M = new TLinearFitter(10,"hyp9");
541 fDFitter2M = new TLinearFitter(10,"hyp9");
542 fDFitter0T = new TLinearFitter(10,"hyp9");
543 fDFitter1T = new TLinearFitter(10,"hyp9");
544 fDFitter2T = new TLinearFitter(10,"hyp9");
545 doinit=kFALSE;
c1418a4c 546 }
547 //}
548
10757ee9 549 fAcceptedTracks++;
550 AddTrack(seed);
551}
552
553Long64_t AliTPCcalibTracksGain::Merge(TCollection *list) {
554 //
555 // Merge() merges the results of all AliTPCcalibTracksGain objects contained in
556 // list, thus allowing a distributed computation of several files, e.g. on PROOF.
557 // The merged results are merged with the data members of the AliTPCcalibTracksGain
558 // object used for calling the Merge method.
559 // The return value is 0 /*the total number of tracks used for calibration*/ if the merge
560 // is successful, otherwise it is -1.
561 //
562
563 if (!list || list->IsEmpty()) return -1;
564
565 if (!fSimpleFitter) fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
566 if (!fSqrtFitter) fSqrtFitter = new AliTPCFitPad(8, "hyp7", "");
567 if (!fLogFitter) fLogFitter = new AliTPCFitPad(8, "hyp7", "");
568 if (!fSingleSectorFitter) fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
569
10757ee9 570
571 // just for debugging
0d3279d4 572 if (!fDebugCalPadRaw) fDebugCalPadRaw = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
573 if (!fDebugCalPadCorr) fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
10757ee9 574
575 TIterator* iter = list->MakeIterator();
576 AliTPCcalibTracksGain* cal = 0;
577
578 while ((cal = (AliTPCcalibTracksGain*)iter->Next())) {
579 if (!cal->InheritsFrom(AliTPCcalibTracksGain::Class())) {
580 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
581 return -1;
582 }
0d3279d4 583
10757ee9 584 Add(cal);
585 }
586 return 0;
587}
588
2acad464 589Float_t AliTPCcalibTracksGain::GetGain(AliTPCclusterMI* cl){
590 //
591 // Return local gain at cluster position
592 //
593 Float_t factor = 1;
09bed5bf 594 if(!fGainMap) return factor;
595
2acad464 596 AliTPCCalROC * roc = fGainMap->GetCalROC(cl->GetDetector());
597 Int_t irow = cl->GetRow();
598 if (roc){
599 if (irow < 63) { // IROC
600 factor = roc->GetValue(irow, TMath::Nint(cl->GetPad()));
601 } else { // OROC
602 factor = roc->GetValue(irow - 63, TMath::Nint(cl->GetPad()));
603 }
604 }
605 if (factor<0.1) factor=0.1;
606 return factor;
607}
608
609
610Float_t AliTPCcalibTracksGain::GetMaxNorm(AliTPCclusterMI * cl){
611 //
612 // Get normalized amplituded if the gain map provided
613 //
614 return cl->GetMax()/GetGain(cl);
615}
616
617
618Float_t AliTPCcalibTracksGain::GetQNorm(AliTPCclusterMI * cl){
619 //
620 // Get normalized amplituded if the gain map provided
621 //
622 return cl->GetQ()/GetGain(cl);
623}
624
625
626
10757ee9 627void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
628 //
629 // Adds another AliTPCcalibTracksGain object to this object.
630 //
631
632 fSimpleFitter->Add(cal->fSimpleFitter);
633 fSqrtFitter->Add(cal->fSqrtFitter);
634 fLogFitter->Add(cal->fLogFitter);
635 fSingleSectorFitter->Add(cal->fSingleSectorFitter);
0d3279d4 636 //
637 //
638 //
639 fFitter0M->Add(cal->fFitter0M);
640 fFitter1M->Add(cal->fFitter1M);
641 fFitter2M->Add(cal->fFitter2M);
642 fFitter0T->Add(cal->fFitter0T);
643 fFitter1T->Add(cal->fFitter1T);
644 fFitter2T->Add(cal->fFitter2T);
645 //
684602c8 646 fDFitter0M->Add(cal->fDFitter0M);
647 fDFitter1M->Add(cal->fDFitter1M);
648 fDFitter2M->Add(cal->fDFitter2M);
649 fDFitter0T->Add(cal->fDFitter0T);
650 fDFitter1T->Add(cal->fDFitter1T);
651 fDFitter2T->Add(cal->fDFitter2T);
652 //
0d3279d4 653 //
654 // histograms
655 //
656 for (Int_t i=0; i<73; i++){
657 TH1F *his,*hism;
658 his = (TH1F*)fArrayQM->At(i);
659 hism = (TH1F*)cal->fArrayQM->At(i);
660 if (his && hism) his->Add(hism);
661 his = (TH1F*)fArrayQT->At(i);
662 hism = (TH1F*)cal->fArrayQT->At(i);
663 if (his && hism) his->Add(hism);
664 }
665 //
666 //
667 for (Int_t i=0; i<37; i++){
668 TProfile *his,*hism;
669 his = (TProfile*)fProfileArrayQM->At(i);
670 hism = (TProfile*)cal->fProfileArrayQM->At(i);
671 if (his && hism) his->Add(hism);
672 his = (TProfile*)fProfileArrayQT->At(i);
673 hism = (TProfile*)cal->fProfileArrayQT->At(i);
674 if (his && hism) his->Add(hism);
675 }
676 //
677 //
678 for (Int_t i=0; i<37; i++){
679 TProfile2D *his,*hism;
680 his = (TProfile2D*)fProfileArrayQM2D->At(i);
681 hism = (TProfile2D*)cal->fProfileArrayQM2D->At(i);
682 if (his && hism) his->Add(hism);
683 his = (TProfile2D*)fProfileArrayQT2D->At(i);
684 hism = (TProfile2D*)cal->fProfileArrayQT2D->At(i);
685 if (his && hism) his->Add(hism);
686 }
687 //
10757ee9 688 // this will be gone for the a new ROOT version > v5-17-05
689 for (UInt_t iSegment = 0; iSegment < 36; iSegment++) {
690 fNShortClusters[iSegment] += cal->fNShortClusters[iSegment];
691 fNMediumClusters[iSegment] += cal->fNMediumClusters[iSegment];
692 fNLongClusters[iSegment] += cal->fNLongClusters[iSegment];
693 }
694
695 // just for debugging, remove me
696 fTotalTracks += cal->fTotalTracks;
697 fAcceptedTracks += cal->fAcceptedTracks;
698 fDebugCalPadRaw->Add(cal->fDebugCalPadRaw);
699 fDebugCalPadCorr->Add(cal->fDebugCalPadCorr);
700
10757ee9 701}
702
703void AliTPCcalibTracksGain::AddTrack(AliTPCseed* seed) {
704 //
705 // The clusters making up the track (seed) are added to various fit functions.
706 // See AddCluster(...) for more detail.
707 //
708
0d3279d4 709 DumpTrack(seed);
710 //
711 // simple histograming part
712 for (Int_t i=0; i<159; i++){
713 AliTPCclusterMI* cluster = seed->GetClusterPointer(i);
714 if (cluster) AddCluster(cluster);
715 }
10757ee9 716}
717
0d3279d4 718void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster){
719 //
720 // Adding cluster information to the simple histograms
2acad464 721 // No correction, fittings are applied
722
723
0d3279d4 724 Float_t kThreshold=5;
725 if (cluster->GetX()<=0) return;
726 if (cluster->GetQ()<=kThreshold) return;
727 //
0d3279d4 728
729 Int_t sector = cluster->GetDetector();
730 TH1F * his;
731 his = GetQT(sector);
2acad464 732 if (his) his->Fill(GetQNorm(cluster));
0d3279d4 733 his = GetQT(-1);
2acad464 734 if (his) his->Fill(GetQNorm(cluster));
0d3279d4 735 his = GetQM(sector);
2acad464 736 if (his) his->Fill(GetMaxNorm(cluster));
0d3279d4 737 his = GetQM(-1);
2acad464 738 if (his) his->Fill(GetMaxNorm(cluster));
0d3279d4 739 //
740 sector = sector%36;
741 TProfile *prof;
742 prof = GetProfileQT(sector);
2acad464 743 if (prof) prof->Fill(cluster->GetX(),GetQNorm(cluster));
0d3279d4 744 prof = GetProfileQT(-1);
2acad464 745 if (prof) prof->Fill(cluster->GetX(),GetQNorm(cluster));
0d3279d4 746 prof = GetProfileQM(sector);
2acad464 747 if (prof) prof->Fill(cluster->GetX(),GetMaxNorm(cluster));
0d3279d4 748 prof = GetProfileQM(-1);
2acad464 749 if (prof) prof->Fill(cluster->GetX(),GetMaxNorm(cluster));
0d3279d4 750 //
751 Float_t phi = cluster->GetY()/cluster->GetX();
752 TProfile2D *prof2;
753 prof2 = GetProfileQT2D(sector);
2acad464 754 if (prof2) prof2->Fill(cluster->GetX(),phi,GetQNorm(cluster));
0d3279d4 755 prof2 = GetProfileQT2D(-1);
2acad464 756 if (prof2) prof2->Fill(cluster->GetX(),phi,GetQNorm(cluster));
0d3279d4 757 prof2 = GetProfileQM2D(sector);
2acad464 758 if (prof2) prof2->Fill(cluster->GetX(),phi,GetMaxNorm(cluster));
0d3279d4 759 prof2 = GetProfileQM2D(-1);
2acad464 760 if (prof2) prof2->Fill(cluster->GetX(),phi,GetMaxNorm(cluster));
0d3279d4 761
762 //
763}
764
765
766
f1c2a4a3 767void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t /*momenta*/, Float_t/* mdedx*/, Int_t padType,
768 Float_t xcenter, TVectorD& dedxQ, TVectorD& /*dedxM*/, Float_t /*fraction*/, Float_t fraction2, Float_t dedge,
769 TVectorD& /*parY*/, TVectorD& /*parZ*/, TVectorD& meanPos) {
10757ee9 770 //
771 // Adds cluster to the appropriate fitter for later analysis.
772 // The charge used for the fit is the maximum charge for this specific cluster or the
773 // accumulated charge per cluster, depending on the value of fgkUseTotalCharge.
774 // Depending on the pad size where the cluster is registered, the value will be put in
775 // the appropriate fitter. Furthermore, for each pad size three different types of fitters
776 // are used. The fit functions are the same for all fitters (parabolic functions), but the value
777 // added to each fitter is different. The simple fitter gets the charge plugged in as is, the sqrt fitter
778 // gets the square root of the charge, and the log fitter gets fgkM*(1+q/fgkM), where q is the original charge
779 // and fgkM==25.
780 //
2acad464 781 Float_t kedge = 3;
782 Float_t kfraction = 0.7;
783 Int_t kinorm = 2;
784
785
786 // Where to put selection on threshold?
787 // Defined by the Q/dEdxT variable - see debug streamer:
788 //
789 // Debug stream variables: (Where tu cut ?)
790 // chain0->Draw("Cl.fQ/dedxQ.fElements[1]>>his(100,0,3)","fraction2<0.6&&dedge>3","",1000000);
791 // mean 1 sigma 0.25
792 // chain0->Draw("Cl.fMax/dedxM.fElements[1]>>his(100,0,3)","fraction2<0.6&&dedge>3","",1000000)
793 // mean 1 sigma 0.25
794 // chain0->Draw("Cl.fQ/dedxQ.fElements[2]>>his(100,0,3)","fraction2<0.7&&dedge>3","",1000000)
795 // mean 1 sigma 0.29
796 // chain0->Draw("Cl.fMax/dedxM.fElements[2]>>his(100,0,3)","fraction2<0.7&&dedge>3","",1000000)
797 // mean 1 sigma 0.27
798 // chain0->Draw("Cl.fQ/dedxQ.fElements[3]>>his(100,0,3)","fraction2<0.8&&dedge>3","",1000000)
799 // mean 1 sigma 0.32
800 //
801 // chain0->Draw("Cl.fQ/dedxQ.fElements[4]>>his(100,0,3)","fraction2<0.9&&dedge>3","",1000000)
802 // mean 1 sigma 0.4
803
804 // Fraction choosen 0.7
10757ee9 805
806 if (!cluster) {
807 Error("AddCluster", "Cluster not valid.");
808 return;
809 }
810
2acad464 811 if (dedge < kedge) return;
812 if (fraction2 > kfraction) return;
10757ee9 813
814 //Int_t padType = GetPadType(cluster->GetX());
815 Double_t xx[7];
816 //Double_t centerPad[2] = {0};
817 //AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
818 //xx[0] = cluster->GetX() - centerPad[0];
819 //xx[1] = cluster->GetY() - centerPad[1];
820 xx[0] = cluster->GetX() - xcenter;
821 xx[1] = cluster->GetY();
822 xx[2] = xx[0] * xx[0];
823 xx[3] = xx[1] * xx[1];
824 xx[4] = xx[0] * xx[1];
825 xx[5] = TMath::Abs(cluster->GetZ()) - TMath::Abs(meanPos[4]);
826 xx[6] = xx[5] * xx[5];
0d3279d4 827 //
828 // Update profile histograms
829 //
10757ee9 830
0d3279d4 831 //
832 // Update fitters
833 //
10757ee9 834 Int_t segment = cluster->GetDetector() % 36;
2acad464 835 Double_t q = fgkUseTotalCharge ? ((Double_t)(GetQNorm(cluster))) : ((Double_t)(GetMaxNorm(cluster))); // note: no normalization to pad size!
10757ee9 836
837 // just for debugging
838 Int_t row = 0;
839 Int_t pad = 0;
840 GetRowPad(cluster->GetX(), cluster->GetY(), row, pad);
841 fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
842
843 // correct charge by normalising to mean charge per track
2acad464 844 q /= dedxQ[kinorm];
10757ee9 845
846 // just for debugging
847 fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
848
849 Double_t sqrtQ = TMath::Sqrt(q);
850 Double_t logQ = fgkM * TMath::Log(1 + q / fgkM);
684602c8 851 TLinearFitter * fitter =0;
852 //
853 fitter = fSimpleFitter->GetFitter(segment, padType);
854 fitter->AddPoint(xx, q);
855 //
856 fitter = fSqrtFitter->GetFitter(segment, padType);
857 fitter->AddPoint(xx, sqrtQ);
858 //
859 fitter = fLogFitter->GetFitter(segment, padType);
860 fitter->AddPoint(xx, logQ);
861 //
862 fitter=fSingleSectorFitter->GetFitter(0, padType);
863 fitter->AddPoint(xx, q);
10757ee9 864
865 // this will be gone for the a new ROOT version > v5-17-05
866 if (padType == kShortPads)
867 fNShortClusters[segment]++;
0d3279d4 868 if (padType == kMediumPads)
10757ee9 869 fNMediumClusters[segment]++;
0d3279d4 870 if (padType == kLongPads)
10757ee9 871 fNLongClusters[segment]++;
872}
873
874void AliTPCcalibTracksGain::Evaluate(Bool_t robust, Double_t frac) {
875 //
876 // Evaluates all fitters contained in this object.
877 // If the robust option is set to kTRUE a robust fit is performed with frac as
878 // the minimal fraction of good points (see TLinearFitter::EvalRobust for details).
879 // Beware: Robust fitting is much slower!
880 //
881
882 fSimpleFitter->Evaluate(robust, frac);
883 fSqrtFitter->Evaluate(robust, frac);
884 fLogFitter->Evaluate(robust, frac);
885 fSingleSectorFitter->Evaluate(robust, frac);
0d3279d4 886 fFitter0M->Eval();
887 fFitter1M->Eval();
888 fFitter2M->Eval();
889 fFitter0T->Eval();
890 fFitter1T->Eval();
891 fFitter2T->Eval();
684602c8 892 //
893 fDFitter0M->Eval();
894 fDFitter1M->Eval();
895 fDFitter2M->Eval();
896 fDFitter0T->Eval();
897 fDFitter1T->Eval();
898 fDFitter2T->Eval();
10757ee9 899}
900
901AliTPCCalPad* AliTPCcalibTracksGain::CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
902 //
903 // Creates the calibration object AliTPCcalPad using fitted parameterization
904 //
905 TObjArray tpc(72);
906 for (UInt_t iSector = 0; iSector < 72; iSector++)
907 tpc.Add(CreateFitCalROC(iSector, fitType, undoTransformation, normalizeToPadSize));
908 return new AliTPCCalPad(&tpc);
909}
910
911AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
912 //
913 // Create the AliTPCCalROC with the values per pad
914 // sector - sector of interest
915 // fitType -
916 //
917
918 TVectorD par(8);
919 if (sector < 36) {
920 GetParameters(sector % 36, 0, fitType, par);
921 return CreateFitCalROC(sector, 0, par, fitType, undoTransformation, normalizeToPadSize);
922 }
923 else {
924 GetParameters(sector % 36, 1, fitType, par);
925 AliTPCCalROC* roc1 = CreateFitCalROC(sector, 1, par, fitType, undoTransformation, normalizeToPadSize);
926 GetParameters(sector % 36, 2, fitType, par);
927 AliTPCCalROC* roc2 = CreateFitCalROC(sector, 2, par, fitType, undoTransformation, normalizeToPadSize);
928 AliTPCCalROC* roc3 = CreateCombinedCalROC(roc1, roc2);
929 delete roc1;
930 delete roc2;
931 return roc3;
932 }
933}
934
935AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t padType, TVectorD &fitParam, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
936 //
937 // This function is essentially a copy of AliTPCCalROC::CreateGlobalFitCalROC(...), with the
938 // modifications, that the center of the region of same pad size is used as the origin
939 // of the fit function instead of the center of the ROC.
940 // The possibility of a linear fit is removed as well because it is not needed.
941 // Only values for pads with the given pad size are calculated, the rest is 0.
942 // Set undoTransformation for undoing the transformation that was applied to the
943 // charge values before they were put into the fitter (thus allowing comparison to the original
944 // charge values). For fitType use 0 for the simple fitter, 1 for the sqrt fitter, 2 for the log fitter.
945 // If normalizeToPadSize is true, the values are normalized to the pad size.
946 // Please be aware, that you even need to specify the fitType if you want to normalize to the pad size without
947 // undoing the transformation (because normalizing involves undoing the trafo first, then normalizing, then
948 // applying the trafo again).
949 // Please note: The normalization to the pad size is a simple linear scaling with the pad length, which
950 // actually doesn't describe reality!
951 //
952
953 Float_t dlx, dly;
954 Double_t centerPad[2] = {0};
955 Float_t localXY[3] = {0};
956 AliTPCROC* tpcROC = AliTPCROC::Instance();
957 if ((padType == 0 && sector >= tpcROC->GetNInnerSector()) || (padType > 0 && sector < tpcROC->GetNInnerSector()) || sector >= tpcROC->GetNSector())
958 return 0;
959 AliTPCCalROC* lROCfitted = new AliTPCCalROC(sector);
960 //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
961 UInt_t startRow = 0;
962 UInt_t endRow = 0;
963 switch (padType) {
964 case kShortPads:
965 startRow = 0;
966 endRow = lROCfitted->GetNrows();
967 break;
968 case kMediumPads:
969 startRow = 0;
970 endRow = 64;
971 break;
972 case kLongPads:
973 startRow = 64;
974 endRow = lROCfitted->GetNrows();
975 break;
976 }
977
978 AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
979 Double_t value = 0;
980 for (UInt_t irow = startRow; irow < endRow; irow++) {
981 for (UInt_t ipad = 0; ipad < lROCfitted->GetNPads(irow); ipad++) {
982 tpcROC->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number
983 dlx = localXY[0] - centerPad[0];
984 dly = localXY[1] - centerPad[1];
985 value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
986
987 // Let q' = value be the transformed value without any pad size corrections,
988 // let T be the transformation and let l be the pad size
989 // 1) don't undo transformation, don't normalize: return q'
990 // 2) undo transformation, don't normalize: return T^{-1} q'
991 // 3) undo transformation, normalize: return (T^{-1} q') / l
992 // 4) don't undo transformation, normalize: return T((T^{-1} q') / l)
993 if (!undoTransformation && !normalizeToPadSize) {/* value remains unchanged */} // (1)
994 else { // (2), (3), (4)
995 //calculate T^{-1}
996 switch (fitType) {
997 case 0: /* value remains unchanged */ break;
998 case 1: value = value * value; break;
999 case 2: value = (TMath::Exp(value / fgkM) - 1) * fgkM; break;
1000 default: Error("CreateFitCalROC", "Wrong fit type."); break;
1001 }
1002 if (normalizeToPadSize) value /= GetPadLength(localXY[0]); // (3)
1003 }
1004 if (!undoTransformation && normalizeToPadSize) { // (4)
1005 // calculate T
1006 switch (fitType) {
1007 case 0: /* value remains unchanged */ break;
1008 case 1: value = TMath::Sqrt(value); break;
1009 case 2: value = fgkM * TMath::Log(1 + value / fgkM); break;
1010 default: Error("CreateFitCalROC", "Wrong fit type."); break;
1011 }
1012 }
1013 lROCfitted->SetValue(irow, ipad, value);
1014 }
1015 }
1016 return lROCfitted;
1017}
1018
1019AliTPCCalROC* AliTPCcalibTracksGain::CreateCombinedCalROC(const AliTPCCalROC* roc1, const AliTPCCalROC* roc2) {
1020 //
1021 // Combines the medium pad size values of roc1 with the long pad size values of roc2 into a new
1022 // AliTPCCalROC. Returns a null pointer if any one of the ROCs is an IROC; issues a warning message
1023 // if the sectors of roc1 and roc2 don't match, but still continue and use the sector of roc1 as the
1024 // sector of the new ROC.
1025 //
1026
1027 if (!roc1 || !roc2) return 0;
1028 if ((Int_t)(roc1->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
1029 if ((Int_t)(roc2->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
1030 if (roc1->GetSector() != roc2->GetSector()) Warning("CreateCombinedCalROC", "Sector number mismatch.");
1031 AliTPCCalROC* roc = new AliTPCCalROC(roc1->GetSector());
1032
1033 for (UInt_t iRow = 0; iRow < 64; iRow++) {
1034 for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
1035 roc->SetValue(iRow, iPad, roc1->GetValue(iRow, iPad));
1036 }
1037 for (UInt_t iRow = 64; iRow < roc->GetNrows(); iRow++) {
1038 for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
1039 roc->SetValue(iRow, iPad, roc2->GetValue(iRow, iPad));
1040 }
1041 return roc;
1042}
1043
684602c8 1044Bool_t AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) {
10757ee9 1045 //
1046 // Puts the fit parameters for the specified segment (IROC & OROC), padType and fitType
1047 // into the fitParam TVectorD (which should contain 8 elements).
1048 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1049 // Note: The fitter has to be evaluated first!
1050 //
684602c8 1051 TLinearFitter * fitter = GetFitter(segment, padType, fitType);
1052 if (fitter){
1053 fitter->Eval();
1054 fitter->GetParameters(fitParam);
1055 return kTRUE;
1056 }else{
1057 Error("AliTPCcalibTracksGain::GetParameters",
1058 Form("Fitter%d_%d_%d not availble", segment, padType, fitType));
1059 return kFALSE;
1060 }
1061 return kFALSE;
10757ee9 1062}
1063
1064void AliTPCcalibTracksGain::GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError) {
1065 //
1066 // Puts the fit parameter errors for the specified segment (IROC & OROC), padType and fitType
1067 // into the fitError TVectorD (which should contain 8 elements).
1068 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1069 // Note: The fitter has to be evaluated first!
1070 //
1071
1072 GetFitter(segment, padType, fitType)->GetErrors(fitError);
1073 fitError *= TMath::Sqrt(GetRedChi2(segment, padType, fitType));
1074}
1075
1076Double_t AliTPCcalibTracksGain::GetRedChi2(UInt_t segment, UInt_t padType, UInt_t fitType) {
1077 //
1078 // Returns the reduced chi^2 value for the specified segment, padType and fitType.
1079 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1080 // Note: The fitter has to be evaluated first!
1081 //
1082
1083 // this will be gone for the a new ROOT version > v5-17-05
1084 Int_t lNClusters = 0;
1085 switch (padType) {
1086 case kShortPads:
1087 lNClusters = fNShortClusters[segment];
1088 break;
1089 case kMediumPads:
1090 lNClusters = fNMediumClusters[segment];
1091 break;
1092 case kLongPads:
1093 lNClusters = fNLongClusters[segment];
1094 break;
1095 }
1096 return GetFitter(segment, padType, fitType)->GetChisquare()/(lNClusters - 8);
1097}
1098
1099void AliTPCcalibTracksGain::GetCovarianceMatrix(UInt_t segment, UInt_t padType, UInt_t fitType, TMatrixD& covMatrix) {
1100 //
1101 // Returns the covariance matrix for the specified segment, padType, fitType.
1102 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1103 //
1104
1105 GetFitter(segment, padType, fitType)->GetCovarianceMatrix(covMatrix);
1106}
1107
1108TLinearFitter* AliTPCcalibTracksGain::GetFitter(UInt_t segment, UInt_t padType, UInt_t fitType) {
1109 //
1110 // Returns the TLinearFitter object for the specified segment, padType, fitType.
1111 // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1112 //
1113
1114 switch (fitType) {
1115 case kSimpleFitter:
1116 return fSimpleFitter->GetFitter(segment, padType);
1117 case kSqrtFitter:
1118 return fSqrtFitter->GetFitter(segment, padType);
1119 case kLogFitter:
1120 return fLogFitter->GetFitter(segment, padType);
1121 case 3:
1122 return fSingleSectorFitter->GetFitter(0, padType);
1123 }
1124 return 0;
1125}
1126
1127Double_t AliTPCcalibTracksGain::GetPadLength(Double_t lx) {
1128 //
1129 // The function returns 0.75 for an IROC, 1. for an OROC at medium pad size position,
1130 // 1.5 for an OROC at long pad size position, -1 if out of bounds.
1131 //
1132
1133 Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
1134 Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
1135 Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
1136 Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
1137 Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
1138 Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
1139
1140 // if IROC
1141 if (lx >= irocLow && lx <= irocUp) return 0.75;
1142 // if OROC medium pads
1143 if (lx >= orocLow1 && lx <= orocUp1) return 1.;
1144 // if OROC long pads
1145 if (lx >= orocLow2 && lx <= orocUp2) return 1.5;
1146 // if out of bounds
1147 return -1;
1148}
1149
1150Int_t AliTPCcalibTracksGain::GetPadType(Double_t lx) {
1151 //
1152 // The function returns 0 for an IROC, 1 for an OROC at medium pad size position,
1153 // 2 for an OROC at long pad size position, -1 if out of bounds.
1154 //
1155
1156 if (GetPadLength(lx) == 0.75) return 0;
1157 else if (GetPadLength(lx) == 1.) return 1;
1158 else if (GetPadLength(lx) == 1.5) return 2;
1159 return -1;
1160}
1161
1162// ONLY FOR DEBUGGING PURPOSES - REMOVE ME WHEN NOT NEEDED ANYMORE
1163Bool_t AliTPCcalibTracksGain::GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad) {
1164 //
1165 // Calculate the row and pad number when the local coordinates are given.
1166 // Returns kFALSE if the position is out of range, otherwise return kTRUE.
1167 // WARNING: This function is preliminary and probably isn't very accurate!!
1168 //
1169
1170 Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
1171 //Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
1172 Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
1173 //Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
1174 Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
1175 //Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
1176
1177 if (GetPadType(lx) == 0) {
1178 row = (Int_t)((lx - irocLow) / fgTPCparam->GetInnerPadPitchLength());
1179 pad = (Int_t)((ly + fgTPCparam->GetYInner(row)) / fgTPCparam->GetInnerPadPitchWidth());
1180 } else if (GetPadType(lx) == 1) {
1181 row = (Int_t)((lx - orocLow1) / fgTPCparam->GetOuter1PadPitchLength());
1182 pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
1183 } else if (GetPadType(lx) == 2) {
1184 row = fgTPCparam->GetNRowUp1() + (Int_t)((lx - orocLow2) / fgTPCparam->GetOuter2PadPitchLength());
1185 pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
1186 }
1187 else return kFALSE;
1188 return kTRUE;
1189}
1190
1191void AliTPCcalibTracksGain::DumpTrack(AliTPCseed* track) {
1192 //
1193 // Dump track information to the debug stream
1194 //
1195
10757ee9 1196 Int_t rows[200];
0d3279d4 1197 Int_t sector[3];
1198 Int_t npoints[3];
1199 TVectorD dedxM[3];
1200 TVectorD dedxQ[3];
1201 TVectorD parY[3];
1202 TVectorD parZ[3];
1203 TVectorD meanPos[3];
1204 //
1205 Int_t count=0;
10757ee9 1206 for (Int_t ipad = 0; ipad < 3; ipad++) {
0d3279d4 1207 dedxM[ipad].ResizeTo(5);
1208 dedxQ[ipad].ResizeTo(5);
1209 parY[ipad].ResizeTo(3);
1210 parZ[ipad].ResizeTo(3);
1211 meanPos[ipad].ResizeTo(6);
1212 Bool_t isOK = GetDedx(track, ipad, rows, sector[ipad], npoints[ipad], dedxM[ipad], dedxQ[ipad], parY[ipad], parZ[ipad], meanPos[ipad]);
1213 if (isOK)
1214 AddTracklet(sector[ipad],ipad, dedxQ[ipad], dedxM[ipad], parY[ipad], parZ[ipad], meanPos[ipad] );
1215 if (isOK) count++;
10757ee9 1216 }
ae28e92e 1217
1218 TTreeSRedirector * cstream = GetDebugStreamer();
1219 if (cstream){
1220 (*cstream) << "Track" <<
0d3279d4 1221 "Track.=" << track << // track information
0d3279d4 1222 "\n";
ae28e92e 1223 //
1224 //
1225 //
1226 if ( GetStreamLevel()>1 && count>1){
1227 (*cstream) << "TrackG" <<
1228 "Track.=" << track << // track information
1229 "ncount="<<count<<
1230 // info for pad type 0
1231 "sector0="<<sector[0]<<
1232 "npoints0="<<npoints[0]<<
1233 "dedxM0.="<<&dedxM[0]<<
1234 "dedxQ0.="<<&dedxQ[0]<<
1235 "parY0.="<<&parY[0]<<
1236 "parZ0.="<<&parZ[0]<<
1237 "meanPos0.="<<&meanPos[0]<<
1238 //
1239 // info for pad type 1
1240 "sector1="<<sector[1]<<
1241 "npoints1="<<npoints[1]<<
1242 "dedxM1.="<<&dedxM[1]<<
1243 "dedxQ1.="<<&dedxQ[1]<<
1244 "parY1.="<<&parY[1]<<
1245 "parZ1.="<<&parZ[1]<<
1246 "meanPos1.="<<&meanPos[1]<<
1247 //
1248 // info for pad type 2
1249 "sector2="<<sector[2]<<
1250 "npoints2="<<npoints[2]<<
1251 "dedxM2.="<<&dedxM[2]<<
1252 "dedxQ2.="<<&dedxQ[2]<<
1253 "parY2.="<<&parY[2]<<
1254 "parZ2.="<<&parZ[2]<<
1255 "meanPos2.="<<&meanPos[2]<<
1256 //
1257 "\n";
1258 }
0d3279d4 1259 }
0d3279d4 1260}
1261
f1c2a4a3 1262Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /*rows*/,
0d3279d4 1263 Int_t &sector, Int_t& npoints,
1264 TVectorD &dedxM, TVectorD &dedxQ,
1265 TVectorD &parY, TVectorD &parZ, TVectorD&meanPos)
1266{
1267 //
1268 // GetDedx for given sector for given track
1269 // padType - type of pads
1270 //
1271
1272 static TLinearFitter fitY(2, "pol1");
1273 static TLinearFitter fitZ(2, "pol1");
684602c8 1274 fitY.StoreData(kFALSE);
1275 fitZ.StoreData(kFALSE);
0d3279d4 1276 //
1277 //
10757ee9 1278 Int_t firstRow = 0, lastRow = 0;
1279 Int_t minRow = 100;
1280 Float_t xcenter = 0;
1281 const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
1282 const Float_t kedgey = 4.;
1283 if (padType == 0) {
1284 firstRow = 0;
1285 lastRow = fgTPCparam->GetNRowLow();
1286 xcenter = 108.47;
1287 }
1288 if (padType == 1) {
1289 firstRow = fgTPCparam->GetNRowLow();
1290 lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1291 xcenter = 166.60;
1292 }
1293 if (padType == 2) {
1294 firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1295 lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
1296 xcenter = 222.6;
1297 }
1298 minRow = (lastRow - firstRow) / 2;
1299 //
1300 //
1301 Int_t nclusters = 0;
1302 Int_t nclustersNE = 0; // number of not edge clusters
1303 Int_t lastSector = -1;
1304 Float_t amplitudeQ[100];
1305 Float_t amplitudeM[100];
1306 Int_t rowIn[100];
1307 Int_t index[100];
1308 //
0d3279d4 1309 //
10757ee9 1310 fitY.ClearPoints();
1311 fitZ.ClearPoints();
10757ee9 1312
1313 for (Int_t iCluster = firstRow; iCluster < lastRow; iCluster++) {
1314 AliTPCclusterMI* cluster = track->GetClusterPointer(iCluster);
1315 if (cluster) {
1316 Int_t detector = cluster->GetDetector() ;
1317 if (lastSector == -1) lastSector = detector;
1318 if (lastSector != detector) continue;
2acad464 1319 amplitudeQ[nclusters] = GetQNorm(cluster);
1320 amplitudeM[nclusters] = GetMaxNorm(cluster);
10757ee9 1321 rowIn[nclusters] = iCluster;
1322 nclusters++;
1323 Double_t dx = cluster->GetX() - xcenter;
1324 Double_t y = cluster->GetY();
1325 Double_t z = cluster->GetZ();
1326 fitY.AddPoint(&dx, y);
1327 fitZ.AddPoint(&dx, z);
1328 meanPos[0] += dx;
1329 meanPos[1] += dx;
1330 meanPos[2] += y;
1331 meanPos[3] += y*y;
1332 meanPos[4] += z;
1333 meanPos[5] += z*z;
1334 if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) nclustersNE++;
1335 }
1336 }
1337
1338 if (nclusters < minRow / 2) return kFALSE;
1339 if (nclustersNE < minRow / 2) return kFALSE;
1340 for (Int_t i = 0; i < 6; i++) meanPos[i] /= Double_t(nclusters);
1341 fitY.Eval();
1342 fitZ.Eval();
1343 fitY.GetParameters(parY);
1344 fitZ.GetParameters(parZ);
1345 //
1346 // calculate truncated mean
1347 //
1348 TMath::Sort(nclusters, amplitudeQ, index, kFALSE);
0d3279d4 1349 //
1350 //
1351 //
10757ee9 1352 Float_t ndedx[5];
1353 for (Int_t i = 0; i < 5; i++) {
1354 dedxQ[i] = 0;
1355 dedxM[i] = 0;
1356 ndedx[i] = 0;
1357 }
1358 //
1359 // dedx calculation
1360 //
1361 Int_t inonEdge = 0;
1362 for (Int_t i = 0; i < nclusters; i++) {
1363 Int_t rowSorted = rowIn[index[i]];
1364 AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1365
1366 if (TMath::Abs(cluster->GetY()) > cluster->GetX()*ktany - kedgey) continue; //don't take edge clusters
1367 inonEdge++;
1368 if (inonEdge < nclustersNE * 0.5) {
1369 ndedx[0]++;
1370 dedxQ[0] += amplitudeQ[index[i]];
1371 dedxM[0] += amplitudeM[index[i]];
1372 }
1373 if (inonEdge < nclustersNE * 0.6) {
1374 ndedx[1]++;
1375 dedxQ[1] += amplitudeQ[index[i]];
1376 dedxM[1] += amplitudeM[index[i]];
1377 }
1378 if (inonEdge < nclustersNE * 0.7) {
1379 ndedx[2]++;
1380 dedxQ[2] += amplitudeQ[index[i]];
1381 dedxM[2] += amplitudeM[index[i]];
1382 }
1383 if (inonEdge < nclustersNE * 0.8) {
1384 ndedx[3]++;
1385 dedxQ[3] += amplitudeQ[index[i]];
1386 dedxM[3] += amplitudeM[index[i]];
1387 }
1388 if (inonEdge < nclustersNE * 0.9) {
1389 ndedx[4]++;
1390 dedxQ[4] += amplitudeQ[index[i]];
1391 dedxM[4] += amplitudeM[index[i]];
1392 }
1393 }
1394 for (Int_t i = 0; i < 5; i++) {
1395 dedxQ[i] /= ndedx[i];
1396 dedxM[i] /= ndedx[i];
1397 }
ae28e92e 1398 TTreeSRedirector * cstream = GetDebugStreamer();
10757ee9 1399 inonEdge = 0;
0d3279d4 1400 Float_t momenta = track->GetP();
1401 Float_t mdedx = track->GetdEdx();
10757ee9 1402 for (Int_t i = 0; i < nclusters; i++) {
1403 Int_t rowSorted = rowIn[index[i]];
1404 AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1405 if (!cluster) {
1406 printf("Problem\n");
1407 continue;
1408 }
1409 if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) inonEdge++;
1410 Float_t dedge = cluster->GetX()*ktany - TMath::Abs(cluster->GetY());
1411 Float_t fraction = Float_t(i) / Float_t(nclusters);
1412 Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
10757ee9 1413
1414 AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
2acad464 1415 Float_t gain = GetGain(cluster);
ae28e92e 1416 if (cstream) (*cstream) << "dEdx" <<
2acad464 1417 "Cl.=" << cluster << // cluster of interest
1418 "gain="<<gain<< // gain at cluster position
1419 "P=" << momenta << // track momenta
1420 "dedx=" << mdedx << // mean dedx - corrected for angle
1421 "IPad=" << padType << // pad type 0..2
1422 "xc=" << xcenter << // x center of chamber
1423 "dedxQ.=" << &dedxQ << // dedxQ - total charge
1424 "dedxM.=" << &dedxM << // dedxM - maximal charge
1425 "fraction=" << fraction << // fraction - order in statistic (0,1)
1426 "fraction2=" << fraction2 << // fraction - order in statistic (0,1)
1427 "dedge=" << dedge << // distance to the edge
1428 "parY.=" << &parY << // line fit
1429 "parZ.=" << &parZ << // line fit
1430 "meanPos.=" << &meanPos << // mean position (dx, dx^2, y,y^2, z, z^2)
1431 "\n";
10757ee9 1432 }
0d3279d4 1433
ae28e92e 1434 if (cstream) (*cstream) << "dEdxT" <<
0d3279d4 1435 "P=" << momenta << // track momenta
1436 "npoints="<<inonEdge<< // number of points
1437 "sector="<<lastSector<< // sector number
1438 "dedx=" << mdedx << // mean dedx - corrected for angle
1439 "IPad=" << padType << // pad type 0..2
1440 "xc=" << xcenter << // x center of chamber
1441 "dedxQ.=" << &dedxQ << // dedxQ - total charge
1442 "dedxM.=" << &dedxM << // dedxM - maximal charge
1443 "parY.=" << &parY << // line fit
1444 "parZ.=" << &parZ << // line fit
1445 "meanPos.=" << &meanPos << // mean position (dx, dx^2, y,y^2, z, z^2)
1446 "\n";
1447
1448 sector = lastSector;
1449 npoints = inonEdge;
10757ee9 1450 return kTRUE;
1451}
0d3279d4 1452
1453void AliTPCcalibTracksGain::AddTracklet(UInt_t sector, UInt_t padType,TVectorD &dedxQ, TVectorD &dedxM,TVectorD& parY, TVectorD& parZ, TVectorD& meanPos){
1454 //
1455 // Add measured point - dedx to the fitter
1456 //
1457 //
8076baa0 1458 //chain->SetAlias("dr","(250-abs(meanPos.fElements[4]))/250");
0d3279d4 1459 //chain->SetAlias("tz","(0+abs(parZ.fElements[1]))");
1460 //chain->SetAlias("ty","(0+abs(parY.fElements[1]))");
1461 //chain->SetAlias("corrg","sqrt((1+ty^2)*(1+tz^2))");
8076baa0 1462 //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 1463
1464 Double_t xxx[100];
1465 //
1466 // z and angular part
1467 //
8076baa0 1468
1469 xxx[0] = (250.-TMath::Abs(meanPos[4]))/250.;
0d3279d4 1470 xxx[1] = TMath::Abs(parY[1]);
1471 xxx[2] = TMath::Abs(parZ[1]);
1472 xxx[3] = xxx[0]*xxx[1];
1473 xxx[4] = xxx[0]*xxx[2];
1474 xxx[5] = xxx[1]*xxx[2];
1475 xxx[6] = xxx[0]*xxx[0];
8076baa0 1476 xxx[7] = xxx[1]*xxx[1];
1477 xxx[8] = xxx[2]*xxx[2];
0d3279d4 1478 //
1479 // chamber part
1480 //
1481 Int_t tsector = sector%36;
1482 for (Int_t i=0;i<35;i++){
8076baa0 1483 xxx[9+i]=(i==tsector)?1:0;
0d3279d4 1484 }
1485 TLinearFitter *fitterM = fFitter0M;
1486 if (padType==1) fitterM=fFitter1M;
1487 if (padType==2) fitterM=fFitter2M;
1488 fitterM->AddPoint(xxx,dedxM[1]);
1489 //
1490 TLinearFitter *fitterT = fFitter0T;
1491 if (padType==1) fitterT = fFitter1T;
1492 if (padType==2) fitterT = fFitter2T;
1493 fitterT->AddPoint(xxx,dedxQ[1]);
684602c8 1494 //
1495 TLinearFitter *dfitterM = fDFitter0M;
1496 if (padType==1) dfitterM=fDFitter1M;
1497 if (padType==2) dfitterM=fDFitter2M;
1498 dfitterM->AddPoint(xxx,dedxM[1]);
1499 //
1500 TLinearFitter *dfitterT = fDFitter0T;
1501 if (padType==1) dfitterT = fDFitter1T;
1502 if (padType==2) dfitterT = fDFitter2T;
1503 dfitterT->AddPoint(xxx,dedxQ[1]);
0d3279d4 1504}
8076baa0 1505
1506
1507TGraph *AliTPCcalibTracksGain::CreateAmpGraph(Int_t ipad, Bool_t qmax){
1508 //
1509 // create the amplitude graph
1510 // The normalized amplitudes are extrapolated to the 0 angle (y,z) and 0 drift length
1511 //
1512
1513 TVectorD vec;
1514 if (qmax){
1515 if (ipad==0) fFitter0M->GetParameters(vec);
1516 if (ipad==1) fFitter1M->GetParameters(vec);
1517 if (ipad==2) fFitter2M->GetParameters(vec);
1518 }else{
1519 if (ipad==0) fFitter0T->GetParameters(vec);
1520 if (ipad==1) fFitter1T->GetParameters(vec);
1521 if (ipad==2) fFitter2T->GetParameters(vec);
1522 }
1523
1524 Float_t amp[36];
1525 Float_t sec[36];
1526 for (Int_t i=0;i<35;i++){
1527 sec[i]=i;
1528 amp[i]=vec[10+i]+vec[0];
1529 }
1530 amp[35]=vec[0];
1531 Float_t mean = TMath::Mean(36,amp);
1532 for (Int_t i=0;i<36;i++){
1533 sec[i]=i;
1534 amp[i]=(amp[i]-mean)/mean;
1535 }
1536 TGraph *gr = new TGraph(36,sec,amp);
1537 return gr;
1538}
1539
1540
1541void AliTPCcalibTracksGain::UpdateClusterParam(AliTPCClusterParam* clparam){
1542 //
1543 // SetQ normalization parameters
1544 //
1545 // void SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm);
1546
1547 TVectorD vec;
684602c8 1548
8076baa0 1549 //
684602c8 1550 fDFitter0T->Eval();
1551 fDFitter1T->Eval();
1552 fDFitter2T->Eval();
1553 fDFitter0M->Eval();
1554 fDFitter1M->Eval();
1555 fDFitter2M->Eval();
1556 fDFitter0T->GetParameters(vec);
8076baa0 1557 clparam->SetQnorm(0,0,&vec);
684602c8 1558 fDFitter1T->GetParameters(vec);
8076baa0 1559 clparam->SetQnorm(1,0,&vec);
684602c8 1560 fDFitter2T->GetParameters(vec);
8076baa0 1561 clparam->SetQnorm(2,0,&vec);
1562 //
684602c8 1563 fDFitter0M->GetParameters(vec);
8076baa0 1564 clparam->SetQnorm(0,1,&vec);
684602c8 1565 fDFitter1M->GetParameters(vec);
8076baa0 1566 clparam->SetQnorm(1,1,&vec);
684602c8 1567 fDFitter2M->GetParameters(vec);
8076baa0 1568 clparam->SetQnorm(2,1,&vec);
1569 //
1570
1571}
b8601924 1572
1573
1574void AliTPCcalibTracksGain::Analyze(){
1575
1576 Evaluate();
1577
1578}
1579
1580
da7d274e 1581
1582TVectorD * AliTPCcalibTracksGain::MakeQPosNorm(TTree * chain0, Int_t ipad, Bool_t isMax, Int_t maxPoints, Int_t verbose){
1583 //
1584 // Input parameters
1585 // chain0 - the tree with information -Debug stream
1586 // ipad - 0 IROC
1587 // - 1 OROC medium
1588 // - 2 OROC LONG
1589 // isMax - kFALSE - total charge param
1590 // kTRUE - Max charge param
1591 //
1592 // maxPoints - number of points for fit
1593 //
1594 // verbose -
1595 //
1596 /* e.g
1597 ipad=0
1598 isMax=kTRUE;
1599 maxPoints=1000000;
1600 */
1601 // Make Q normalization as function of following parameters
1602 // 1 - dp - relative pad position
1603 // 2 - dt - relative time position
1604 // 3 - di - drift length (norm to 1);
1605 // 4 - dq0 - Tot/Max charge
1606 // 5 - dq1 - Max/Tot charge
1607 // 6 - sy - sigma y - shape
1608 // 7 - sz - sigma z - shape
1609 //
1610 // Coeficient of Taylor expansion fitted
1611 // Fit parameters returned as TVectorD
1612 // Fit parameters to be used in corresponding correction function
1613 // in AliTPCclusterParam
1614 //
1615 //
1616 TStatToolkit toolkit;
1617 Double_t chi2;
1618 TVectorD fitParam;
1619 TMatrixD covMatrix;
1620 Int_t npoints;
1621 TCut cutA("dedge>3&&fraction2<0.7");
1622 chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1623 chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1624 chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1625 chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1626 chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1627 chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1628 chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1629 //
1630 TString fstring="";
1631 fstring+="dp++"; //1
1632 fstring+="dt++"; //2
1633 fstring+="dp*dp++"; //3
1634 fstring+="dt*dt++"; //4
1635 fstring+="dt*dt*dt++"; //5
1636 fstring+="dp*dt++"; //6
1637 fstring+="dp*dt*dt++"; //7
1638 fstring+="(dq0)++"; //8
1639 fstring+="(dq1)++"; //9
1640 //
1641 //
1642 fstring+="dp*dp*(di)++"; //10
1643 fstring+="dt*dt*(di)++"; //11
1644 fstring+="dp*dp*sy++"; //12
1645 fstring+="dt*sz++"; //13
1646 fstring+="dt*dt*sz++"; //14
1647 fstring+="dt*dt*dt*sz++"; //15
1648 //
1649 fstring+="dp*dp*1*sy*sz++"; //16
1650 fstring+="dt*sy*sz++"; //17
1651 fstring+="dt*dt*sy*sz++"; //18
1652 fstring+="dt*dt*dt*sy*sz++"; //19
1653 //
1654 fstring+="dp*dp*(dq0)++"; //20
1655 fstring+="dt*1*(dq0)++"; //21
1656 fstring+="dt*dt*(dq0)++"; //22
1657 fstring+="dt*dt*dt*(dq0)++"; //23
1658 //
1659 fstring+="dp*dp*(dq1)++"; //24
1660 fstring+="dt*(dq1)++"; //25
1661 fstring+="dt*dt*(dq1)++"; //26
1662 fstring+="dt*dt*dt*(dq1)++"; //27
1663
1664 TString var;
1665 if (isMax) var = "Cl.fMax/gain/dedxM.fElements[2]";
1666 if (!isMax) var = "Cl.fQ/gain/dedxQ.fElements[2]";
1667 TString cutP="IPad==";
1668 cutP+=ipad;
1669 //
1670 TString *strq0 = toolkit.FitPlane(chain0,var.Data(),fstring.Data(), cutP.Data()+cutA, chi2,npoints,fitParam,covMatrix,-1,0,maxPoints);
1671 //
1672 //
1673 if (verbose){
1674 printf("Chi2/npoints = %f",TMath::Sqrt(chi2/npoints));
1675 printf("\nFit function\n:%s\n",strq0->Data());
1676 }
1677 TVectorD *vec = new TVectorD(fitParam);
1678 return vec;
1679}
1680
1681void AliTPCcalibTracksGain::MakeQPosNormAll(TTree * chain, AliTPCClusterParam * param, Int_t maxPoints, Int_t verbose){
1682 //
1683 // Fill the content of the of the AliTPCclusterParam
1684 // with fitted values of corrections
1685 //
1686
1687}
1688
1689
1690
1691/*
1692
1693 Position correction fit:
1694 //
1695TStatToolkit toolkit;
1696Double_t chi2;
1697TVectorD fitParam;
1698TMatrixD covMatrix;
1699Int_t npoints;
1700//
1701TCut cutA("dedge>3&&fraction2<0.7");
1702chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1703chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1704chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1705chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1706chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1707chain0->SetAlias("sy","(0.2/sqrt(0.01^2+Cl.fSigmaY2))");
1708chain0->SetAlias("sz","(0.2/sqrt(0.01^2+Cl.fSigmaZ2))");
1709
1710TString fstring="";
1711
1712fstring+="dp++"; //1
1713fstring+="dt++"; //2
1714fstring+="dp*dp++"; //3
1715fstring+="dt*dt++"; //4
1716fstring+="dt*dt*dt++"; //5
1717fstring+="dp*dt++"; //6
1718fstring+="dp*dt*dt++"; //7
1719fstring+="(dq0)++"; //8
1720fstring+="(dq1)++"; //9
1721//
1722//
1723fstring+="dp*dp*(di)++"; //10
1724fstring+="dt*dt*(di)++"; //11
1725fstring+="dp*dp*sy++"; //12
1726fstring+="dt*sz++"; //13
1727fstring+="dt*dt*sz++"; //14
1728fstring+="dt*dt*dt*sz++"; //15
1729//
1730fstring+="dp*dp*1*sy*sz++"; //16
1731fstring+="dt*sy*sz++"; //17
1732fstring+="dt*dt*sy*sz++"; //18
1733fstring+="dt*dt*dt*sy*sz++"; //19
1734//
1735fstring+="dp*dp*(dq0)++"; //20
1736fstring+="dt*1*(dq0)++"; //21
1737fstring+="dt*dt*(dq0)++"; //22
1738fstring+="dt*dt*dt*(dq0)++"; //23
1739
1740fstring+="dp*dp*(dq1)++"; //24
1741fstring+="dt*(dq1)++"; //25
1742fstring+="dt*dt*(dq1)++"; //26
1743fstring+="dt*dt*dt*(dq1)++"; //27
1744
1745
1746 TString *strq0 = toolkit.FitPlane(chain0,"Cl.fMax/gain/dedxM.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
1747 TString *strqt0 = toolkit.FitPlane(chain0,"Cl.fQ/gain/dedxQ.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
1748
1749 chain0->SetAlias("qcorM0",strq0->Data());
1750 chain0->SetAlias("qcorT0",strqt0->Data());
1751//chain0->SetAlias("mmqcorM0","min(max(qcorM0,0.75),1.15)");
1752 chain0->Draw("(Cl.fMax/gain/dedxM.fElements[2]):min(max(qcorM0,0.75),1.15)","IPad==0"+cutA,"prof",100000)
1753
1754 fraction05 -
1755 sigma 0.2419
1756 sigma fit 0.2302
1757 sigma fit with shape 0.2257
1758 fraction 07
1759 qtot sigma 0.322
1760 qmax sigma 0.292
1761 qmax sigma fit 0.2702
1762 qmax sigma fit+ratio 0.2638
1763
1764*/
1765