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