TPC/AliTPCcalibTracksGain.cxx

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