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