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