TPC/AliTPCcalibTracksGain.cxx

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