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