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