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