TPC/AliTPCcalibTracksGain.cxx

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