TPC/AliTPCcalibTracksGain.cxx

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