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