TPC/AliTPCcalibTracksGain.cxx

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