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 //       === Calibration class for gain calibration using tracks ===
  19 //
  20 //    1) Genereal idea
  21 //    ================
  22 //       A 6-parametric parabolic function
  23 //
  24 //          G(x, y) = p0 + p1*x + p2*y + p3*x^2 + p4*y^2 + p5 * x*y
  25 //
  26 //       is fitted to the maximum charge values or total charge values of
  27 //       all the clusters contained in the tracks that are added to this
  28 //       object. This fit is performed for each read out chamber, in fact even
  29 //       for each type of pad sizes (thus for one segment, which consists of
  30 //       an IROC and an OROC, there are three fitters used, corresponding to
  31 //       the three pad sizes). The coordinate origin is at the center of the
  32 //       particular pad size region on each ROC.
  33 //
  34 //       Because of the Landau nature of the charge deposition we use
  35 //       different "types" of fitters instead of one to minimize the effect
  36 //       of the long Landau tail. The difference between the fitters is only
  37 //       the charge value, that is put into them, i.e. the charge is subject
  38 //       to a transformation. At this point we use three different fit types:
  39 //
  40 //          a) simple: the charge is put in as it is
  41 //          b) sqrt:   the square root of the charge is put into the fitter
  42 //          c) log:    fgkM * Log(1+q/fgkM) is put into the fitter, with
  43 //                     q being the untransformed charge and fgkM=25
  44 //
  45 //       The results of the fits may be visualized and further used by
  46 //       creating an AliTPCCalROC or AliTPCCalPad. You may specify to undo
  47 //       the transformation and/or to normalize to the pad size.
  48 //
  49 //       Not every track you add to this object is actually used for
  50 //       calibration. There are some cuts and conditions to exclude bad
  51 //       tracks, e.g. a pt cut to cut out tracks with too much charge
  52 //       deposition or a cut on edge clusters which are not fully
  53 //       registered and don't give a usable signal.
  54 //
  55 //    2) Interface / usage
  56 //    ====================
  57 //       For each track to be added you need to call Process().
  58 //       This method expects an AliTPCseed, which contains the necessary
  59 //       cluster information. At the moment of writing this information
  60 //       is stored in an AliESDfriend corresponding to an AliESD.
  61 //       You may also call AddTrack() if you don't want the cuts and
  62 //       other quality conditions to kick in (thus forcing the object to
  63 //       accept the track) or AddCluster() for adding single clusters.
  64 //       Call one of the Evaluate functions to evaluate the fitter(s) and
  65 //       to retrieve the fit parameters, erros and so on. You can also
  66 //       do this later on by using the different Getters.
  67 //
  68 //       The visualization methods CreateFitCalPad() and CreateFitCalROC()
  69 //       are straight forward to use.
  70 //
  71 //       Note: If you plan to write this object to a ROOT file, make sure
  72 //             you evaluate all the fitters *before* writing, because due
  73 //             to a bug in the fitter component writing fitters doesn't
  74 //             work properly (yet). Be aware that you cannot re-evaluate
  75 //             the fitters after loading this object from file.
  76 //             (This will be gone for a new ROOT version > v5-17-05)
  77 //
  78 ////////////////////////////////////////////////////////////////////////////
  79
  80 #include <TPDGCode.h>
  81 #include <TStyle.h>
  82 #include "TSystem.h"
  83 #include "TMatrixD.h"
  84 #include "TTreeStream.h"
  85 #include "TF1.h"
  86 #include "AliTPCParamSR.h"
  87 #include "AliTPCClusterParam.h"
  88 #include "AliTrackPointArray.h"
  89 #include "TCint.h"
  90 #include "AliTPCcalibTracksGain.h"
  91 #include <TH1.h>
  92 #include <TH3F.h>
  93 #include <TLinearFitter.h>
  94 #include <TTreeStream.h>
  95 #include <TFile.h>
  96 #include <TCollection.h>
  97 #include <TIterator.h>
  98
  99 //
 100 // AliRoot includes
 101 //
 102 #include "AliMagF.h"
 103 #include "AliMathBase.h"
 104 //
 105 #include "AliTPCROC.h"
 106 #include "AliTPCParamSR.h"
 107 #include "AliTPCCalROC.h"
 108 #include "AliTPCCalPad.h"
 109 //
 110 #include "AliTracker.h"
 111 #include "AliESD.h"
 112 #include "AliESDtrack.h"
 113 #include "AliESDfriend.h"
 114 #include "AliESDfriendTrack.h"
 115 #include "AliTPCseed.h"
 116 #include "AliTPCclusterMI.h"
 117 #include "AliTPCcalibTracksCuts.h"
 118 #include "AliTPCFitPad.h"
 119
 120 // REMOVE ALL OF THIS
 121 #include <TTree.h>
 122 #include "AliESDEvent.h"
 123
 124 /*
 125
 126 TFile f("TPCCalibTracksGain.root")
 127
 128 gSystem->Load("libPWG1.so")
 129 AliTreeDraw comp
 130 comp.SetTree(dEdx)
 131 Double_t chi2;
 132 TVectorD vec(3)
 133 TMatrixD mat(3,3)
 134 TString * str = comp.FitPlane("Cl.fQ/dedxQ.fElements[0]","Cl.fY++Cl.fX","Cl.fDetector<36",chi2,vec,mat)
 135
 136 */
 137
 138 ClassImp(AliTPCcalibTracksGain)
 139
 140 const Bool_t   AliTPCcalibTracksGain::fgkUseTotalCharge = kTRUE;
 141 const Double_t AliTPCcalibTracksGain::fgkM = 25.;
 142 const char*    AliTPCcalibTracksGain::fgkDebugStreamFileName = "TPCCalibTracksGain.root";
 143 AliTPCParamSR* AliTPCcalibTracksGain::fgTPCparam = new AliTPCParamSR();
 144
 145 AliTPCcalibTracksGain::AliTPCcalibTracksGain() :
 146    TNamed(),
 147    fDebugCalPadRaw(0),
 148    fDebugCalPadCorr(0),
 149    fDebugStream(0),
 150    fSimpleFitter(0),
 151    fSqrtFitter(0),
 152    fLogFitter(0),
 153    fSingleSectorFitter(0),
 154    fPrevIter(0),
 155    fDebugStreamPrefix(0),
 156    fCuts(0)
 157 {
 158    //
 159    // Default constructor.
 160    //
 161 }
 162
 163 AliTPCcalibTracksGain::AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj) :
 164    TNamed(obj)
 165 {
 166    //
 167    // Copy constructor.
 168    //
 169
 170    fDebugCalPadRaw = new AliTPCCalPad(*(obj.fDebugCalPadRaw));
 171    fDebugCalPadCorr = new AliTPCCalPad(*(obj.fDebugCalPadCorr));
 172    fSimpleFitter = new AliTPCFitPad(*(obj.fSimpleFitter));
 173    fSqrtFitter = new AliTPCFitPad(*(obj.fSqrtFitter));
 174    fLogFitter = new AliTPCFitPad(*(obj.fLogFitter));
 175    fSingleSectorFitter = new AliTPCFitPad(*(obj.fSingleSectorFitter));
 176    fPrevIter = new AliTPCcalibTracksGain(*(obj.fPrevIter));
 177    fDebugStreamPrefix = new TObjString(*(obj.fDebugStreamPrefix));
 178    fCuts = new AliTPCcalibTracksCuts(*(obj.fCuts));
 179 }
 180
 181 AliTPCcalibTracksGain& AliTPCcalibTracksGain::operator=(const AliTPCcalibTracksGain& rhs) {
 182    //
 183    // Assignment operator.
 184    //
 185
 186    if (this != &rhs) {
 187       TNamed::operator=(rhs);
 188       fDebugCalPadRaw = new AliTPCCalPad(*(rhs.fDebugCalPadRaw));
 189       fDebugCalPadCorr = new AliTPCCalPad(*(rhs.fDebugCalPadCorr));
 190       fSimpleFitter = new AliTPCFitPad(*(rhs.fSimpleFitter));
 191       fSqrtFitter = new AliTPCFitPad(*(rhs.fSqrtFitter));
 192       fLogFitter = new AliTPCFitPad(*(rhs.fLogFitter));
 193       fSingleSectorFitter = new AliTPCFitPad(*(rhs.fSingleSectorFitter));
 194       fPrevIter = new AliTPCcalibTracksGain(*(rhs.fPrevIter));
 195       fDebugStreamPrefix = new TObjString(*(rhs.fDebugStreamPrefix));
 196       fCuts = new AliTPCcalibTracksCuts(*(rhs.fCuts));
 197    }
 198    return *this;
 199 }
 200
 201 AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts, TNamed* debugStreamPrefix, AliTPCcalibTracksGain* prevIter) :
 202    TNamed(name, title),
 203    fDebugCalPadRaw(0),
 204    fDebugCalPadCorr(0),
 205    fDebugStream(0),
 206    fSimpleFitter(0),
 207    fSqrtFitter(0),
 208    fLogFitter(0),
 209    fSingleSectorFitter(0),
 210    fPrevIter(0),
 211    fDebugStreamPrefix(0),
 212    fCuts(0)
 213 {
 214    //
 215    // Constructor.
 216    //
 217
 218    G__SetCatchException(0);
 219
 220    fCuts = cuts;
 221    if (debugStreamPrefix) fDebugStreamPrefix = new TObjString(debugStreamPrefix->GetTitle());
 222    fPrevIter = prevIter;
 223
 224    fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
 225    fSqrtFitter   = new AliTPCFitPad(8, "hyp7", "");
 226    fLogFitter    = new AliTPCFitPad(8, "hyp7", "");
 227    fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
 228
 229    // just for debugging
 230    fTotalTracks     = 0;
 231    fAcceptedTracks  = 0;
 232    fDebugCalPadRaw  = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
 233    fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
 234
 235    // this will be gone for the a new ROOT version > v5-17-05
 236    for (UInt_t i = 0; i < 36; i++) {
 237       fNShortClusters[i]  = 0;
 238       fNMediumClusters[i] = 0;
 239       fNLongClusters[i]   = 0;
 240    }
 241  }
 242
 243 AliTPCcalibTracksGain::~AliTPCcalibTracksGain() {
 244    //
 245    // Destructor.
 246    //
 247
 248    Info("Destructor","");
 249    if (fSimpleFitter) delete fSimpleFitter;
 250    if (fSqrtFitter) delete fSqrtFitter;
 251    if (fLogFitter) delete fLogFitter;
 252    if (fSingleSectorFitter) delete fSingleSectorFitter;
 253
 254
 255    if (fDebugStream) {
 256       //fDebugStream->GetFile()->Close();
 257       printf("Deleting debug stream object\n");
 258       delete fDebugStream;
 259    }
 260
 261    if (fDebugStreamPrefix) delete fDebugStreamPrefix;
 262
 263    if (fDebugCalPadRaw) delete fDebugCalPadRaw;
 264    if (fDebugCalPadCorr) delete fDebugCalPadCorr;
 265 }
 266
 267 void AliTPCcalibTracksGain::Terminate(){
 268    //
 269    // Evaluate fitters and close the debug stream.
 270    // Also move or copy the debug stream, if a debugStreamPrefix is provided.
 271    //
 272
 273    Evaluate();
 274    if (fDebugStream) {
 275      delete fDebugStream;
 276      fDebugStream = 0;
 277    }
 278
 279    if (fDebugStreamPrefix) {
 280       TString debugStreamPrefix = fDebugStreamPrefix->GetString();
 281       TString destFile("");
 282       destFile += debugStreamPrefix;
 283       destFile += "/";
 284       destFile += gSystem->HostName();
 285       destFile += "_TPCCalibTracksGain.root";
 286       if (debugStreamPrefix.BeginsWith("root://")) {
 287          TFile::Cp("TPCCalibTracksGain.root", destFile.Data());
 288       } else {
 289          TString command("mv TPCCalibTracksGain.root ");
 290          command += destFile;
 291          gSystem->Exec(command.Data());
 292       }
 293    }
 294 }
 295
 296 void AliTPCcalibTracksGain::AddInfo(TChain* chain, char* debugStreamPrefix, char* prevIterFileName) {
 297    //
 298    // Add some parameters to the chain.
 299    // debugStreamPrefix: If specified, contains the location (either normal or xrootd directory)
 300    //                    where the debug stream is moved (normal directory) or copied to (xrootd).
 301    // prevIterFileName: If specified, contains an AliTPCcalibTracksGain object from a previous run
 302    //                   for doing an iterative calibration procedure (right now unused).
 303    // Note: The parameters are *not* added to this class, you need to do it later by retrieving
 304    // the parameters from the chain and passing them to the constructor!
 305    //
 306
 307    if (debugStreamPrefix) {
 308       TNamed* objDebugStreamPrefix = new TNamed("debugStreamPrefix", debugStreamPrefix);
 309       chain->GetUserInfo()->AddLast((TObject*)objDebugStreamPrefix);
 310    }
 311
 312    if (prevIterFileName) {
 313       TFile paramFile(prevIterFileName);
 314       if (paramFile.IsZombie()) {
 315          printf("File %s not found. Continuing without previous iteration.\n", prevIterFileName);
 316          return;
 317       }
 318
 319       AliTPCcalibTracksGain *prevIter = (AliTPCcalibTracksGain*)paramFile.Get("calibTracksGain");
 320       if (prevIter) {
 321          chain->GetUserInfo()->AddLast((TObject*)prevIter);
 322       } else
 323          printf("No calibTracksGain object found. Continuing without previous iteration.\n");
 324    }
 325 }
 326
 327 Bool_t AliTPCcalibTracksGain::AcceptTrack(AliTPCseed* track) {
 328    //
 329    // Decides whether to accept a track or not depending on track parameters and cuts
 330    // contained as AliTPCcalibTracksCuts object fCuts.
 331    // Tracks are discarded if the number of clusters is too low or the transverse
 332    // momentum is too low.
 333    // The corresponding cut values are specified in the fCuts member.
 334    //
 335
 336    if (track->GetNumberOfClusters() < fCuts->GetMinClusters()) return kFALSE;
 337    //if ((TMath::Abs(track->GetY() / track->GetX()) > fCuts->GetEdgeYXCutNoise())
 338    //   && (TMath::Abs(track->GetTgl()) < fCuts->GetEdgeThetaCutNoise())) return kFALSE;
 339    //if (track->GetNumberOfClusters() / (track->GetNFoundable()+1.) < fCuts->GetMinRatio()) return kFALSE;
 340    if (TMath::Abs(track->GetSigned1Pt()) > fCuts->GetMax1pt()) return kFALSE;
 341
 342    //if (track->GetPt() < 50.) return kFALSE;
 343    return kTRUE;
 344 }
 345
 346 void AliTPCcalibTracksGain::Process(AliTPCseed* seed) {
 347    //
 348    // Main method to be called when a new seed is supposed to be processed
 349    // and be used for gain calibration. Its quality is checked before it
 350    // is added.
 351    //
 352
 353    fTotalTracks++;
 354    if (!AcceptTrack(seed)) return;
 355    fAcceptedTracks++;
 356    AddTrack(seed);
 357 }
 358
 359 Long64_t AliTPCcalibTracksGain::Merge(TCollection *list) {
 360    //
 361    // Merge() merges the results of all AliTPCcalibTracksGain objects contained in
 362    // list, thus allowing a distributed computation of several files, e.g. on PROOF.
 363    // The merged results are merged with the data members of the AliTPCcalibTracksGain
 364    // object used for calling the Merge method.
 365    // The return value is 0 /*the total number of tracks used for calibration*/ if the merge
 366    // is successful, otherwise it is -1.
 367    //
 368
 369    if (!list || list->IsEmpty()) return -1;
 370
 371    if (!fSimpleFitter) fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
 372    if (!fSqrtFitter)   fSqrtFitter   = new AliTPCFitPad(8, "hyp7", "");
 373    if (!fLogFitter)    fLogFitter    = new AliTPCFitPad(8, "hyp7", "");
 374    if (!fSingleSectorFitter) fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
 375
 376    // this will be gone for the a new ROOT version > v5-17-05
 377    for (UInt_t i = 0; i < 36; i++) {
 378       fNShortClusters[i]  = 0;
 379       fNMediumClusters[i] = 0;
 380       fNLongClusters[i]   = 0;
 381    }
 382
 383    // just for debugging
 384    if (fDebugCalPadRaw)  delete fDebugCalPadRaw;
 385    if (fDebugCalPadCorr) delete fDebugCalPadCorr;
 386    fDebugCalPadRaw  = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
 387    fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
 388    fTotalTracks     = 0;
 389    fAcceptedTracks  = 0;
 390
 391    TIterator* iter = list->MakeIterator();
 392    AliTPCcalibTracksGain* cal = 0;
 393
 394    while ((cal = (AliTPCcalibTracksGain*)iter->Next())) {
 395       if (!cal->InheritsFrom(AliTPCcalibTracksGain::Class())) {
 396          Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 397          return -1;
 398       }
 399       Add(cal);
 400    }
 401    return 0;
 402 }
 403
 404 void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
 405    //
 406    // Adds another AliTPCcalibTracksGain object to this object.
 407    //
 408
 409    fSimpleFitter->Add(cal->fSimpleFitter);
 410    fSqrtFitter->Add(cal->fSqrtFitter);
 411    fLogFitter->Add(cal->fLogFitter);
 412    fSingleSectorFitter->Add(cal->fSingleSectorFitter);
 413
 414    // this will be gone for the a new ROOT version > v5-17-05
 415    for (UInt_t iSegment = 0; iSegment < 36; iSegment++) {
 416       fNShortClusters[iSegment] += cal->fNShortClusters[iSegment];
 417       fNMediumClusters[iSegment] += cal->fNMediumClusters[iSegment];
 418       fNLongClusters[iSegment] += cal->fNLongClusters[iSegment];
 419    }
 420
 421    // just for debugging, remove me
 422    fTotalTracks += cal->fTotalTracks;
 423    fAcceptedTracks += cal->fAcceptedTracks;
 424    fDebugCalPadRaw->Add(cal->fDebugCalPadRaw);
 425    fDebugCalPadCorr->Add(cal->fDebugCalPadCorr);
 426
 427    // Let's see later what to do with fCuts, fDebugStream and fDebugStreamPrefix,
 428    // we'll probably won't merge them.
 429 }
 430
 431 void AliTPCcalibTracksGain::AddTrack(AliTPCseed* seed) {
 432    //
 433    // The clusters making up the track (seed) are added to various fit functions.
 434    // See AddCluster(...) for more detail.
 435    //
 436
 437    if (!fDebugStream) fDebugStream = new TTreeSRedirector(fgkDebugStreamFileName);
 438    DumpTrack(seed);
 439 }
 440
 441 void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t momenta, Float_t mdedx, Int_t padType,
 442             Float_t xcenter, TVectorD dedxQ, TVectorD dedxM, Float_t fraction, Float_t fraction2, Float_t dedge,
 443             TVectorD parY, TVectorD parZ, TVectorD meanPos) {
 444    //
 445    // Adds cluster to the appropriate fitter for later analysis.
 446    // The charge used for the fit is the maximum charge for this specific cluster or the
 447    // accumulated charge per cluster, depending on the value of fgkUseTotalCharge.
 448    // Depending on the pad size where the cluster is registered, the value will be put in
 449    // the appropriate fitter. Furthermore, for each pad size three different types of fitters
 450    // are used. The fit functions are the same for all fitters (parabolic functions), but the value
 451    // added to each fitter is different. The simple fitter gets the charge plugged in as is, the sqrt fitter
 452    // gets the square root of the charge, and the log fitter gets fgkM*(1+q/fgkM), where q is the original charge
 453    // and fgkM==25.
 454    //
 455
 456    if (!cluster) {
 457       Error("AddCluster", "Cluster not valid.");
 458       return;
 459    }
 460
 461    if (dedge < 3.) return;
 462    if (fraction2 > 0.7) return;
 463
 464    //Int_t padType = GetPadType(cluster->GetX());
 465    Double_t xx[7];
 466    //Double_t centerPad[2] = {0};
 467    //AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
 468    //xx[0] = cluster->GetX() - centerPad[0];
 469    //xx[1] = cluster->GetY() - centerPad[1];
 470    xx[0] = cluster->GetX() - xcenter;
 471    xx[1] = cluster->GetY();
 472    xx[2] = xx[0] * xx[0];
 473    xx[3] = xx[1] * xx[1];
 474    xx[4] = xx[0] * xx[1];
 475    xx[5] = TMath::Abs(cluster->GetZ()) - TMath::Abs(meanPos[4]);
 476    xx[6] = xx[5] * xx[5];
 477
 478    Int_t segment = cluster->GetDetector() % 36;
 479    Double_t q = fgkUseTotalCharge ? ((Double_t)(cluster->GetQ())) : ((Double_t)(cluster->GetMax()));  // note: no normalization to pad size!
 480
 481    // just for debugging
 482    Int_t row = 0;
 483    Int_t pad = 0;
 484    GetRowPad(cluster->GetX(), cluster->GetY(), row, pad);
 485    fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
 486
 487    // correct charge by normalising to mean charge per track
 488    q /= dedxQ[2];
 489
 490    // just for debugging
 491    fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
 492
 493    Double_t sqrtQ = TMath::Sqrt(q);
 494    Double_t logQ = fgkM * TMath::Log(1 + q / fgkM);
 495    fSimpleFitter->GetFitter(segment, padType)->AddPoint(xx, q);
 496    fSqrtFitter->GetFitter(segment, padType)->AddPoint(xx, sqrtQ);
 497    fLogFitter->GetFitter(segment, padType)->AddPoint(xx, logQ);
 498    fSingleSectorFitter->GetFitter(0, padType)->AddPoint(xx, q);
 499
 500    // this will be gone for the a new ROOT version > v5-17-05
 501    if (padType == kShortPads)
 502       fNShortClusters[segment]++;
 503    else if (padType == kMediumPads)
 504       fNMediumClusters[segment]++;
 505    else if (padType == kLongPads)
 506       fNLongClusters[segment]++;
 507 }
 508
 509 void AliTPCcalibTracksGain::Evaluate(Bool_t robust, Double_t frac) {
 510    //
 511    // Evaluates all fitters contained in this object.
 512    // If the robust option is set to kTRUE a robust fit is performed with frac as
 513    // the minimal fraction of good points (see TLinearFitter::EvalRobust for details).
 514    // Beware: Robust fitting is much slower!
 515    //
 516
 517    fSimpleFitter->Evaluate(robust, frac);
 518    fSqrtFitter->Evaluate(robust, frac);
 519    fLogFitter->Evaluate(robust, frac);
 520    fSingleSectorFitter->Evaluate(robust, frac);
 521 }
 522
 523 AliTPCCalPad* AliTPCcalibTracksGain::CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 524   //
 525   // Creates the calibration object AliTPCcalPad using fitted parameterization
 526   //
 527    TObjArray tpc(72);
 528    for (UInt_t iSector = 0; iSector < 72; iSector++)
 529       tpc.Add(CreateFitCalROC(iSector, fitType, undoTransformation, normalizeToPadSize));
 530    return new AliTPCCalPad(&tpc);
 531 }
 532
 533 AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 534   //
 535   // Create the AliTPCCalROC with the values per pad
 536   // sector  - sector of interest
 537   // fitType -
 538   //
 539
 540    TVectorD par(8);
 541    if (sector < 36) {
 542       GetParameters(sector % 36, 0, fitType, par);
 543       return CreateFitCalROC(sector, 0, par, fitType, undoTransformation, normalizeToPadSize);
 544    }
 545    else {
 546       GetParameters(sector % 36, 1, fitType, par);
 547       AliTPCCalROC* roc1 = CreateFitCalROC(sector, 1, par, fitType, undoTransformation, normalizeToPadSize);
 548       GetParameters(sector % 36, 2, fitType, par);
 549       AliTPCCalROC* roc2 = CreateFitCalROC(sector, 2, par, fitType, undoTransformation, normalizeToPadSize);
 550       AliTPCCalROC* roc3 = CreateCombinedCalROC(roc1, roc2);
 551       delete roc1;
 552       delete roc2;
 553       return roc3;
 554    }
 555 }
 556
 557 AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t padType, TVectorD &fitParam, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 558    //
 559    // This function is essentially a copy of AliTPCCalROC::CreateGlobalFitCalROC(...), with the
 560    // modifications, that the center of the region of same pad size is used as the origin
 561    // of the fit function instead of the center of the ROC.
 562    // The possibility of a linear fit is removed as well because it is not needed.
 563    // Only values for pads with the given pad size are calculated, the rest is 0.
 564    // Set undoTransformation for undoing the transformation that was applied to the
 565    // charge values before they were put into the fitter (thus allowing comparison to the original
 566    // charge values). For fitType use 0 for the simple fitter, 1 for the sqrt fitter, 2 for the log fitter.
 567    // If normalizeToPadSize is true, the values are normalized to the pad size.
 568    // Please be aware, that you even need to specify the fitType if you want to normalize to the pad size without
 569    // undoing the transformation (because normalizing involves undoing the trafo first, then normalizing, then
 570    // applying the trafo again).
 571    // Please note: The normalization to the pad size is a simple linear scaling with the pad length, which
 572    //              actually doesn't describe reality!
 573    //
 574
 575    Float_t dlx, dly;
 576    Double_t centerPad[2] = {0};
 577    Float_t localXY[3] = {0};
 578    AliTPCROC* tpcROC = AliTPCROC::Instance();
 579    if ((padType == 0 && sector >= tpcROC->GetNInnerSector()) || (padType > 0 && sector < tpcROC->GetNInnerSector()) || sector >= tpcROC->GetNSector())
 580       return 0;
 581    AliTPCCalROC* lROCfitted = new AliTPCCalROC(sector);
 582    //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
 583    UInt_t startRow = 0;
 584    UInt_t endRow = 0;
 585    switch (padType) {
 586       case kShortPads:
 587          startRow = 0;
 588          endRow = lROCfitted->GetNrows();
 589          break;
 590       case kMediumPads:
 591          startRow = 0;
 592          endRow = 64;
 593          break;
 594       case kLongPads:
 595          startRow = 64;
 596          endRow = lROCfitted->GetNrows();
 597          break;
 598    }
 599
 600    AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
 601    Double_t value = 0;
 602    for (UInt_t irow = startRow; irow < endRow; irow++) {
 603       for (UInt_t ipad = 0; ipad < lROCfitted->GetNPads(irow); ipad++) {
 604          tpcROC->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 605          dlx = localXY[0] - centerPad[0];
 606          dly = localXY[1] - centerPad[1];
 607          value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
 608
 609          // Let q' = value be the transformed value without any pad size corrections,
 610          // let T be the transformation and let l be the pad size
 611          //    1) don't undo transformation, don't normalize: return q'
 612          //    2) undo transformation,       don't normalize: return T^{-1} q'
 613          //    3) undo transformation,       normalize:       return (T^{-1} q') / l
 614          //    4) don't undo transformation, normalize:       return T((T^{-1} q') / l)
 615          if (!undoTransformation && !normalizeToPadSize) {/* value remains unchanged */}  // (1)
 616          else {                                                                           // (2), (3), (4)
 617             //calculate T^{-1}
 618             switch (fitType) {
 619                case  0: /* value remains unchanged */ break;
 620                case  1: value = value * value; break;
 621                case  2: value = (TMath::Exp(value / fgkM) - 1) * fgkM; break;
 622                default: Error("CreateFitCalROC", "Wrong fit type."); break;
 623             }
 624             if (normalizeToPadSize) value /= GetPadLength(localXY[0]);                    // (3)
 625          }
 626          if (!undoTransformation && normalizeToPadSize) {                                 // (4)
 627             // calculate T
 628             switch (fitType) {
 629                case  0: /* value remains unchanged */ break;
 630                case  1: value = TMath::Sqrt(value); break;
 631                case  2: value = fgkM * TMath::Log(1 + value / fgkM); break;
 632                default: Error("CreateFitCalROC", "Wrong fit type."); break;
 633             }
 634          }
 635          lROCfitted->SetValue(irow, ipad, value);
 636       }
 637    }
 638    return lROCfitted;
 639 }
 640
 641 AliTPCCalROC* AliTPCcalibTracksGain::CreateCombinedCalROC(const AliTPCCalROC* roc1, const AliTPCCalROC* roc2) {
 642    //
 643    // Combines the medium pad size values of roc1 with the long pad size values of roc2 into a new
 644    // AliTPCCalROC. Returns a null pointer if any one of the ROCs is an IROC; issues a warning message
 645    // if the sectors of roc1 and roc2 don't match, but still continue and use the sector of roc1 as the
 646    // sector of the new ROC.
 647    //
 648
 649    if (!roc1 || !roc2) return 0;
 650    if ((Int_t)(roc1->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
 651    if ((Int_t)(roc2->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
 652    if (roc1->GetSector() != roc2->GetSector()) Warning("CreateCombinedCalROC", "Sector number mismatch.");
 653    AliTPCCalROC* roc = new AliTPCCalROC(roc1->GetSector());
 654
 655    for (UInt_t iRow = 0; iRow < 64; iRow++) {
 656       for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
 657          roc->SetValue(iRow, iPad, roc1->GetValue(iRow, iPad));
 658    }
 659    for (UInt_t iRow = 64; iRow < roc->GetNrows(); iRow++) {
 660       for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
 661          roc->SetValue(iRow, iPad, roc2->GetValue(iRow, iPad));
 662    }
 663    return roc;
 664 }
 665
 666 void AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) {
 667    //
 668    // Puts the fit parameters for the specified segment (IROC & OROC), padType and fitType
 669    // into the fitParam TVectorD (which should contain 8 elements).
 670    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 671    // Note: The fitter has to be evaluated first!
 672    //
 673
 674    GetFitter(segment, padType, fitType)->GetParameters(fitParam);
 675 }
 676
 677 void AliTPCcalibTracksGain::GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError) {
 678    //
 679    // Puts the fit parameter errors for the specified segment (IROC & OROC), padType and fitType
 680    // into the fitError TVectorD (which should contain 8 elements).
 681    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 682    // Note: The fitter has to be evaluated first!
 683    //
 684
 685    GetFitter(segment, padType, fitType)->GetErrors(fitError);
 686    fitError *= TMath::Sqrt(GetRedChi2(segment, padType, fitType));
 687 }
 688
 689 Double_t AliTPCcalibTracksGain::GetRedChi2(UInt_t segment, UInt_t padType, UInt_t fitType) {
 690    //
 691    // Returns the reduced chi^2 value for the specified segment, padType and fitType.
 692    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 693    // Note: The fitter has to be evaluated first!
 694    //
 695
 696    // this will be gone for the a new ROOT version > v5-17-05
 697    Int_t lNClusters = 0;
 698    switch (padType) {
 699       case kShortPads:
 700          lNClusters = fNShortClusters[segment];
 701          break;
 702       case kMediumPads:
 703          lNClusters = fNMediumClusters[segment];
 704          break;
 705       case kLongPads:
 706          lNClusters = fNLongClusters[segment];
 707          break;
 708    }
 709    return GetFitter(segment, padType, fitType)->GetChisquare()/(lNClusters - 8);
 710 }
 711
 712 void AliTPCcalibTracksGain::GetCovarianceMatrix(UInt_t segment, UInt_t padType, UInt_t fitType, TMatrixD& covMatrix) {
 713    //
 714    // Returns the covariance matrix for the specified segment, padType, fitType.
 715    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 716    //
 717
 718    GetFitter(segment, padType, fitType)->GetCovarianceMatrix(covMatrix);
 719 }
 720
 721 TLinearFitter* AliTPCcalibTracksGain::GetFitter(UInt_t segment, UInt_t padType, UInt_t fitType) {
 722    //
 723    // Returns the TLinearFitter object for the specified segment, padType, fitType.
 724    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 725    //
 726
 727    switch (fitType) {
 728       case kSimpleFitter:
 729          return fSimpleFitter->GetFitter(segment, padType);
 730       case kSqrtFitter:
 731          return fSqrtFitter->GetFitter(segment, padType);
 732       case kLogFitter:
 733          return fLogFitter->GetFitter(segment, padType);
 734       case 3:
 735          return fSingleSectorFitter->GetFitter(0, padType);
 736    }
 737    return 0;
 738 }
 739
 740 Double_t AliTPCcalibTracksGain::GetPadLength(Double_t lx) {
 741    //
 742    // The function returns 0.75 for an IROC, 1. for an OROC at medium pad size position,
 743    // 1.5 for an OROC at long pad size position, -1 if out of bounds.
 744    //
 745
 746    Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
 747    Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
 748    Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
 749    Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
 750    Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
 751    Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
 752
 753    // if IROC
 754    if (lx >= irocLow && lx <= irocUp) return 0.75;
 755    // if OROC medium pads
 756    if (lx >= orocLow1 && lx <= orocUp1) return 1.;
 757    // if OROC long pads
 758    if (lx >= orocLow2 && lx <= orocUp2) return 1.5;
 759    // if out of bounds
 760    return -1;
 761 }
 762
 763 Int_t AliTPCcalibTracksGain::GetPadType(Double_t lx) {
 764    //
 765    // The function returns 0 for an IROC, 1 for an OROC at medium pad size position,
 766    // 2 for an OROC at long pad size position, -1 if out of bounds.
 767    //
 768
 769    if (GetPadLength(lx) == 0.75) return 0;
 770    else if (GetPadLength(lx) == 1.) return 1;
 771    else if (GetPadLength(lx) == 1.5) return 2;
 772    return -1;
 773 }
 774
 775 // ONLY FOR DEBUGGING PURPOSES - REMOVE ME WHEN NOT NEEDED ANYMORE
 776 Bool_t AliTPCcalibTracksGain::GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad) {
 777    //
 778    // Calculate the row and pad number when the local coordinates are given.
 779    // Returns kFALSE if the position is out of range, otherwise return kTRUE.
 780    // WARNING: This function is preliminary and probably isn't very accurate!!
 781    //
 782
 783    Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
 784    //Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
 785    Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
 786    //Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
 787    Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
 788    //Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
 789
 790    if (GetPadType(lx) == 0) {
 791       row = (Int_t)((lx - irocLow) / fgTPCparam->GetInnerPadPitchLength());
 792       pad = (Int_t)((ly + fgTPCparam->GetYInner(row)) / fgTPCparam->GetInnerPadPitchWidth());
 793    } else if (GetPadType(lx) == 1) {
 794       row = (Int_t)((lx - orocLow1) / fgTPCparam->GetOuter1PadPitchLength());
 795       pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
 796    } else if (GetPadType(lx) == 2) {
 797       row = fgTPCparam->GetNRowUp1() + (Int_t)((lx - orocLow2) / fgTPCparam->GetOuter2PadPitchLength());
 798       pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
 799    }
 800    else return kFALSE;
 801    return kTRUE;
 802 }
 803
 804 void AliTPCcalibTracksGain::DumpTrack(AliTPCseed* track) {
 805    //
 806    //  Dump track information to the debug stream
 807    //
 808
 809    (*fDebugStream) << "Track" <<
 810       "Track.=" << track <<        // track information
 811       "\n";
 812    Int_t rows[200];
 813    for (Int_t ipad = 0; ipad < 3; ipad++) {
 814       GetDedx(track, ipad, rows);
 815    }
 816 }
 817
 818 Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* rows) {
 819    //
 820    // GetDedx for given sector for given track
 821    // padType - type of pads
 822    //
 823
 824    Int_t firstRow = 0, lastRow = 0;
 825    Int_t minRow = 100;
 826    Float_t xcenter = 0;
 827    const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
 828    const Float_t kedgey = 4.;
 829    if (padType == 0) {
 830       firstRow = 0;
 831       lastRow = fgTPCparam->GetNRowLow();
 832       xcenter = 108.47;
 833    }
 834    if (padType == 1) {
 835       firstRow = fgTPCparam->GetNRowLow();
 836       lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
 837       xcenter = 166.60;
 838    }
 839    if (padType == 2) {
 840       firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
 841       lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
 842       xcenter = 222.6;
 843    }
 844    minRow = (lastRow - firstRow) / 2;
 845    //
 846    //
 847    Int_t nclusters = 0;
 848    Int_t nclustersNE = 0; // number of not edge clusters
 849    Int_t lastSector = -1;
 850    Float_t amplitudeQ[100];
 851    Float_t amplitudeM[100];
 852    Int_t rowIn[100];
 853    Int_t index[100];
 854    //
 855    static TLinearFitter fitY(2, "pol1");
 856    static TLinearFitter fitZ(2, "pol1");
 857    static TVectorD parY(2);
 858    static TVectorD parZ(2);
 859    fitY.ClearPoints();
 860    fitZ.ClearPoints();
 861    TVectorD meanPos(6);
 862
 863    for (Int_t iCluster = firstRow; iCluster < lastRow; iCluster++) {
 864       AliTPCclusterMI* cluster = track->GetClusterPointer(iCluster);
 865       if (cluster) {
 866          Int_t detector = cluster->GetDetector() ;
 867          if (lastSector == -1) lastSector = detector;
 868          if (lastSector != detector) continue;
 869          amplitudeQ[nclusters] = cluster->GetQ();
 870          amplitudeM[nclusters] = cluster->GetMax();
 871          rowIn[nclusters] = iCluster;
 872          nclusters++;
 873          Double_t dx = cluster->GetX() - xcenter;
 874          Double_t y = cluster->GetY();
 875          Double_t z = cluster->GetZ();
 876          fitY.AddPoint(&dx, y);
 877          fitZ.AddPoint(&dx, z);
 878          meanPos[0] += dx;
 879          meanPos[1] += dx;
 880          meanPos[2] += y;
 881          meanPos[3] += y*y;
 882          meanPos[4] += z;
 883          meanPos[5] += z*z;
 884          if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) nclustersNE++;
 885       }
 886    }
 887
 888    if (nclusters < minRow / 2) return kFALSE;
 889    if (nclustersNE < minRow / 2) return kFALSE;
 890    for (Int_t i = 0; i < 6; i++) meanPos[i] /= Double_t(nclusters);
 891    fitY.Eval();
 892    fitZ.Eval();
 893    fitY.GetParameters(parY);
 894    fitZ.GetParameters(parZ);
 895    //
 896    // calculate truncated mean
 897    //
 898    TMath::Sort(nclusters, amplitudeQ, index, kFALSE);
 899
 900    TVectorD dedxQ(5);
 901    TVectorD dedxM(5);
 902    Float_t ndedx[5];
 903    for (Int_t i = 0; i < 5; i++) {
 904       dedxQ[i] = 0;
 905       dedxM[i] = 0;
 906       ndedx[i] = 0;
 907    }
 908    //
 909    // dedx calculation
 910    //
 911    Int_t inonEdge = 0;
 912    for (Int_t i = 0; i < nclusters; i++) {
 913       Int_t rowSorted = rowIn[index[i]];
 914       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
 915
 916       if (TMath::Abs(cluster->GetY()) > cluster->GetX()*ktany - kedgey) continue;  //don't take edge clusters
 917       inonEdge++;
 918       if (inonEdge < nclustersNE * 0.5) {
 919          ndedx[0]++;
 920          dedxQ[0] += amplitudeQ[index[i]];
 921          dedxM[0] += amplitudeM[index[i]];
 922       }
 923       if (inonEdge < nclustersNE * 0.6) {
 924          ndedx[1]++;
 925          dedxQ[1] += amplitudeQ[index[i]];
 926          dedxM[1] += amplitudeM[index[i]];
 927       }
 928       if (inonEdge < nclustersNE * 0.7) {
 929          ndedx[2]++;
 930          dedxQ[2] += amplitudeQ[index[i]];
 931          dedxM[2] += amplitudeM[index[i]];
 932       }
 933       if (inonEdge < nclustersNE * 0.8) {
 934          ndedx[3]++;
 935          dedxQ[3] += amplitudeQ[index[i]];
 936          dedxM[3] += amplitudeM[index[i]];
 937       }
 938       if (inonEdge < nclustersNE * 0.9) {
 939          ndedx[4]++;
 940          dedxQ[4] += amplitudeQ[index[i]];
 941          dedxM[4] += amplitudeM[index[i]];
 942       }
 943    }
 944    for (Int_t i = 0; i < 5; i++) {
 945       dedxQ[i] /= ndedx[i];
 946       dedxM[i] /= ndedx[i];
 947    }
 948
 949    inonEdge = 0;
 950    for (Int_t i = 0; i < nclusters; i++) {
 951       Int_t rowSorted = rowIn[index[i]];
 952       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
 953       if (!cluster) {
 954          printf("Problem\n");
 955          continue;
 956       }
 957       if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) inonEdge++;
 958       Float_t dedge = cluster->GetX()*ktany - TMath::Abs(cluster->GetY());
 959       Float_t fraction = Float_t(i) / Float_t(nclusters);
 960       Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
 961       Float_t momenta = track->GetP();
 962       Float_t mdedx = track->GetdEdx();
 963
 964       AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
 965
 966       (*fDebugStream) << "dEdx" <<
 967          "Cl.=" << cluster <<           // cluster of interest
 968          "P=" << momenta <<             // track momenta
 969          "dedx=" << mdedx <<            // mean dedx - corrected for angle
 970          "IPad=" << padType <<          // pad type 0..2
 971          "xc=" << xcenter <<            // x center of chamber
 972          "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
 973          "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
 974          "fraction=" << fraction <<     // fraction - order in statistic (0,1)
 975          "fraction2=" << fraction2 <<   // fraction - order in statistic (0,1)
 976          "dedge=" << dedge <<           // distance to the edge
 977          "parY.=" << &parY <<           // line fit
 978          "parZ.=" << &parZ <<           // line fit
 979          "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
 980          "\n";
 981    }
 982    return kTRUE;
 983 }