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