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