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