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