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 //    In order to debug some numerical algorithm all data data which are used for
 101 //    fitters can be stored in the debug streamers. In case of fitting roblems the
 102 //    errors and tendencies can be checked.
 103 //
 104 //    Debug Streamers:
 105 //
 106 //
 107 //
 108 //
 109 //
 110 ////////////////////////////////////////////////////////////////////////////
 111
 112 /*
 113   .x ~/UliStyle.C
 114   gSystem->Load("libANALYSIS");
 115   gSystem->Load("libTPCcalib");
 116   TFile fcalib("CalibObjects.root");
 117   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
 118   AliTPCcalibTracksGain * gain = ( AliTPCcalibTracksGain *)array->FindObject("calibTracksGain");
 119
 120   //
 121   // Make visual Tree - compare with Kr calibration
 122   //
 123   AliTPCCalPad * gain010 = gain->CreateFitCalPad(0,kTRUE,0); gain010->SetName("CGain010");
 124   AliTPCCalPad * gain110 = gain->CreateFitCalPad(1,kTRUE,0); gain110->SetName("CGain110");
 125   AliTPCCalPad * gain210 = gain->CreateFitCalPad(2,kTRUE,0); gain210->SetName("CGain210");
 126   TFile fkr("/u/akalweit/macros/GainCalibration/GainMap.root");
 127   AliTPCCalPad *gainKr = fkr.Get("GainMap"); fkr->SetName("KrGain");
 128   //
 129   AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
 130   preprocesor->AddComponent(gain010);
 131   preprocesor->AddComponent(gain110);
 132   preprocesor->AddComponent(gain210);
 133   preprocesor->AddComponent(gainKr);
 134   preprocesor->DumpToFile("cosmicGain.root");
 135   //
 136   //
 137   //
 138   // Simple session using the debug streamers
 139   //
 140
 141   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 142   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
 143   AliXRDPROOFtoolkit tool;
 144
 145   TChain * chain0 = tool.MakeChain("gain.txt","dEdx",0,1000000);
 146   TChain * chain1 = tool.MakeChain("gain.txt","Track",0,1000000);
 147   TChain * chain2 = tool.MakeChain("gain.txt","TrackG",0,1000000);
 148   chain0->Lookup();
 149   chain1->Lookup();
 150   chain2->Lookup();
 151
 152   chain2->SetAlias("k1","1/0.855");
 153   chain2->SetAlias("k0","1/0.9928");
 154   chain2->SetAlias("k2","1/1.152");
 155
 156
 157
 158
 159 */
 160
 161
 162
 163
 164
 165 #include <TPDGCode.h>
 166 #include <TStyle.h>
 167 #include "TSystem.h"
 168 #include "TMatrixD.h"
 169 #include "TTreeStream.h"
 170 #include "TF1.h"
 171 #include "AliTPCParamSR.h"
 172 #include "AliTPCClusterParam.h"
 173 #include "AliTrackPointArray.h"
 174 #include "TCint.h"
 175 #include "AliTPCcalibTracksGain.h"
 176 #include <TH1.h>
 177 #include <TH3F.h>
 178 #include <TLinearFitter.h>
 179 #include <TTreeStream.h>
 180 #include <TFile.h>
 181 #include <TCollection.h>
 182 #include <TIterator.h>
 183 #include <TProfile.h>
 184 #include <TProfile2D.h>
 185 #include <TProof.h>
 186
 187 //
 188 // AliRoot includes
 189 //
 190 #include "AliMagF.h"
 191 #include "AliMathBase.h"
 192 //
 193 #include "AliTPCROC.h"
 194 #include "AliTPCParamSR.h"
 195 #include "AliTPCCalROC.h"
 196 #include "AliTPCCalPad.h"
 197 #include "AliTPCClusterParam.h"
 198 //
 199 #include "AliTracker.h"
 200 #include "AliESD.h"
 201 #include "AliESDtrack.h"
 202 #include "AliESDfriend.h"
 203 #include "AliESDfriendTrack.h"
 204 #include "AliTPCseed.h"
 205 #include "AliTPCclusterMI.h"
 206 #include "AliTPCcalibTracksCuts.h"
 207 #include "AliTPCFitPad.h"
 208
 209 // REMOVE ALL OF THIS
 210 #include <TTree.h>
 211 #include "AliESDEvent.h"
 212
 213 /*
 214
 215 TFile f("TPCCalibTracksGain.root")
 216
 217 gSystem->Load("libPWG1.so")
 218 AliTreeDraw comp
 219 comp.SetTree(dEdx)
 220 Double_t chi2;
 221 TVectorD vec(3)
 222 TMatrixD mat(3,3)
 223 TString * str = comp.FitPlane("Cl.fQ/dedxQ.fElements[0]","Cl.fY++Cl.fX","Cl.fDetector<36",chi2,vec,mat)
 224
 225 */
 226
 227 ClassImp(AliTPCcalibTracksGain)
 228
 229 const Bool_t   AliTPCcalibTracksGain::fgkUseTotalCharge = kTRUE;
 230 const Double_t AliTPCcalibTracksGain::fgkM = 25.;
 231 const char*    AliTPCcalibTracksGain::fgkDebugStreamFileName = "TPCCalibTracksGain.root";
 232 AliTPCParamSR* AliTPCcalibTracksGain::fgTPCparam = new AliTPCParamSR();
 233
 234 AliTPCcalibTracksGain::AliTPCcalibTracksGain() :
 235   AliTPCcalibBase(),
 236   fCuts(0),            // cuts that are used for sieving the tracks used for calibration
 237   //
 238   // Simple Array of histograms
 239   //
 240   fArrayQM(0),                // Qmax normalized
 241   fArrayQT(0),                // Qtot normalized
 242   fProfileArrayQM(0),         // Qmax normalized  versus local X
 243   fProfileArrayQT(0),         // Qtot normalized  versus local X
 244   fProfileArrayQM2D(0),       // Qmax normalized  versus local X and phi
 245   fProfileArrayQT2D(0),       // Qtot normalized  versus local X and phi
 246   //
 247   // Fitters
 248   //
 249   fSimpleFitter(0),         // simple fitter for short pads
 250   fSqrtFitter(0),           // sqrt fitter for medium pads
 251   fLogFitter(0),            // log fitter for long pads
 252
 253   fFitter0M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 254   fFitter1M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 255   fFitter2M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 256   fFitter0T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 257   fFitter1T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 258   fFitter2T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 259   //
 260   fDFitter0M(0),              // fitting of the atenuation, angular correction
 261   fDFitter1M(0),              // fitting of the atenuation, angular correction
 262   fDFitter2M(0),              // fitting of the atenuation, angular correction
 263   fDFitter0T(0),              // fitting of the atenuation, angular correction
 264   fDFitter1T(0),              // fitting of the atenuation, angular correction
 265   fDFitter2T(0),              // fitting of the atenuation, angular correction
 266   //
 267   fSingleSectorFitter(0),   // just for debugging
 268   //
 269   // Counters
 270   //
 271   fTotalTracks(0),         // just for debugging
 272   fAcceptedTracks(0),      // just for debugging
 273   fDebugCalPadRaw(0),      // just for debugging
 274   fDebugCalPadCorr(0),     // just for debugging
 275   fPrevIter(0)        // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!)
 276
 277 {
 278    //
 279    // Default constructor.
 280    //
 281 }
 282
 283 AliTPCcalibTracksGain::AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj) :
 284   AliTPCcalibBase(obj),
 285   fCuts(obj.fCuts),            // cuts that are used for sieving the tracks used for calibration
 286   fArrayQM(0),                // Qmax normalized
 287   fArrayQT(0),                // Qtot normalized
 288   //
 289   // Simple Histograms
 290   //
 291   fProfileArrayQM(obj.fProfileArrayQM),         // Qmax normalized  versus local X
 292   fProfileArrayQT(obj.fProfileArrayQT),         // Qtot normalized  versus local X
 293   fProfileArrayQM2D(obj.fProfileArrayQM2D),       // Qmax normalized  versus local X and phi
 294   fProfileArrayQT2D(obj.fProfileArrayQT2D),       // Qtot normalized  versus local X and phi
 295   //
 296   // Fitters
 297   //
 298   fSimpleFitter(obj.fSimpleFitter),         // simple fitter for short pads
 299   fSqrtFitter(obj.fSqrtFitter),           // sqrt fitter for medium pads
 300   fLogFitter(obj.fLogFitter),            // log fitter for long pads
 301   fFitter0M(obj.fFitter0M),
 302   fFitter1M(obj.fFitter1M),
 303   fFitter2M(obj.fFitter2M),
 304   fFitter0T(obj.fFitter0T),
 305   fFitter1T(obj.fFitter1T),
 306   fFitter2T(obj.fFitter2T),
 307   //
 308   fDFitter0M(obj.fDFitter0M),
 309   fDFitter1M(obj.fDFitter1M),
 310   fDFitter2M(obj.fDFitter2M),
 311   fDFitter0T(obj.fDFitter0T),
 312   fDFitter1T(obj.fDFitter1T),
 313   fDFitter2T(obj.fDFitter2T),
 314   fSingleSectorFitter(obj.fSingleSectorFitter),   // just for debugging
 315   //
 316   // Counters
 317   //
 318   fTotalTracks(obj.fTotalTracks),         // just for debugging
 319   fAcceptedTracks(obj.fAcceptedTracks),      // just for debugging
 320   fDebugCalPadRaw(obj.fDebugCalPadRaw),      // just for debugging
 321   fDebugCalPadCorr(obj.fDebugCalPadCorr),     // just for debugging
 322   fPrevIter(obj.fPrevIter)        // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!)
 323
 324 {
 325    //
 326    // Copy constructor.
 327    //
 328 }
 329
 330 AliTPCcalibTracksGain& AliTPCcalibTracksGain::operator=(const AliTPCcalibTracksGain& rhs) {
 331    //
 332    // Assignment operator.
 333    //
 334
 335    if (this != &rhs) {
 336       TNamed::operator=(rhs);
 337       fDebugCalPadRaw = new AliTPCCalPad(*(rhs.fDebugCalPadRaw));
 338       fDebugCalPadCorr = new AliTPCCalPad(*(rhs.fDebugCalPadCorr));
 339       fSimpleFitter = new AliTPCFitPad(*(rhs.fSimpleFitter));
 340       fSqrtFitter = new AliTPCFitPad(*(rhs.fSqrtFitter));
 341       fLogFitter = new AliTPCFitPad(*(rhs.fLogFitter));
 342       fSingleSectorFitter = new AliTPCFitPad(*(rhs.fSingleSectorFitter));
 343       fPrevIter = new AliTPCcalibTracksGain(*(rhs.fPrevIter));
 344       fCuts = new AliTPCcalibTracksCuts(*(rhs.fCuts));
 345    }
 346    return *this;
 347 }
 348
 349 AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts, TNamed* /*debugStreamPrefix*/, AliTPCcalibTracksGain* prevIter) :
 350   AliTPCcalibBase(),
 351   fCuts(0),            // cuts that are used for sieving the tracks used for calibration
 352   fArrayQM(0),                // Qmax normalized
 353   fArrayQT(0),                // Qtot normalized
 354   //
 355   // Simple Histograms
 356   //
 357   fProfileArrayQM(0),         // Qmax normalized  versus local X
 358   fProfileArrayQT(0),         // Qtot normalized  versus local X
 359   fProfileArrayQM2D(0),       // Qmax normalized  versus local X and phi
 360   fProfileArrayQT2D(0),       // Qtot normalized  versus local X and phi
 361   //
 362   // Fitters
 363   //
 364   fSimpleFitter(0),         // simple fitter for short pads
 365   fSqrtFitter(0),           // sqrt fitter for medium pads
 366   fLogFitter(0),            // log fitter for long pads
 367   fFitter0M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 368   fFitter1M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 369   fFitter2M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 370   fFitter0T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 371   fFitter1T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 372   fFitter2T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 373   //
 374   fDFitter0M(0),              // fitting of the atenuation, angular correction
 375   fDFitter1M(0),              // fitting of the atenuation, angular correction
 376   fDFitter2M(0),              // fitting of the atenuation, angular correction
 377   fDFitter0T(0),              // fitting of the atenuation, angular correction
 378   fDFitter1T(0),              // fitting of the atenuation, angular correction
 379   fDFitter2T(0),              // fitting of the atenuation, angular correction
 380   fSingleSectorFitter(0),   // just for debugging
 381   //
 382   // Counters
 383   //
 384   fTotalTracks(0),         // just for debugging
 385   fAcceptedTracks(0),      // just for debugging
 386   fDebugCalPadRaw(0),      // just for debugging
 387   fDebugCalPadCorr(0),     // just for debugging
 388   fPrevIter(0)        // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!)
 389
 390 {
 391    //
 392    // Constructor.
 393    //
 394    G__SetCatchException(0);
 395    this->SetNameTitle(name, title);
 396    fCuts = cuts;
 397    fPrevIter = prevIter;
 398    //
 399    // Fitter initialization
 400    //
 401    fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
 402    fSqrtFitter   = new AliTPCFitPad(8, "hyp7", "");
 403    fLogFitter    = new AliTPCFitPad(8, "hyp7", "");
 404    fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
 405    //
 406    fFitter0M      = new TLinearFitter(45,"hyp44");
 407    fFitter1M      = new TLinearFitter(45,"hyp44");
 408    fFitter2M      = new TLinearFitter(45,"hyp44");
 409    fFitter0T      = new TLinearFitter(45,"hyp44");
 410    fFitter1T      = new TLinearFitter(45,"hyp44");
 411    fFitter2T      = new TLinearFitter(45,"hyp44");
 412    //
 413    fDFitter0M      = new TLinearFitter(10,"hyp9");
 414    fDFitter1M      = new TLinearFitter(10,"hyp9");
 415    fDFitter2M      = new TLinearFitter(10,"hyp9");
 416    fDFitter0T      = new TLinearFitter(10,"hyp9");
 417    fDFitter1T      = new TLinearFitter(10,"hyp9");
 418    fDFitter2T      = new TLinearFitter(10,"hyp9");
 419    //
 420    //
 421    fFitter0M->StoreData(kFALSE);
 422    fFitter1M->StoreData(kFALSE);
 423    fFitter2M->StoreData(kFALSE);
 424    fFitter0T->StoreData(kFALSE);
 425    fFitter1T->StoreData(kFALSE);
 426    fFitter2T->StoreData(kFALSE);
 427    //
 428    fDFitter0M->StoreData(kFALSE);
 429    fDFitter1M->StoreData(kFALSE);
 430    fDFitter2M->StoreData(kFALSE);
 431    fDFitter0T->StoreData(kFALSE);
 432    fDFitter1T->StoreData(kFALSE);
 433    fDFitter2T->StoreData(kFALSE);
 434    //
 435    //
 436    // Add profile histograms -JUST for visualization - Not used for real calibration
 437    //
 438    //
 439    fArrayQM=new TObjArray(73);                // Qmax normalized
 440    fArrayQT=new TObjArray(73);                // Qtot normalized
 441    fProfileArrayQM   = new TObjArray(37);         // Qmax normalized  versus local X
 442    fProfileArrayQT   = new TObjArray(37);         // Qtot normalized  versus local X
 443    fProfileArrayQM2D = new TObjArray(37);         // Qmax normalized  versus local X and phi
 444    fProfileArrayQT2D = new TObjArray(37);         // Qtot normalized  versus local X and phi
 445    char hname[1000];
 446    for (Int_t i=0; i<73; i++){
 447      sprintf(hname,"QM_%d",i);
 448      fArrayQM->AddAt(new TH1F(hname,hname,200,0,1000),i);
 449      sprintf(hname,"QT_%d",i);
 450      fArrayQT->AddAt(new TH1F(hname,hname,200,0,1000),i);
 451    }
 452
 453    for (Int_t i=0; i<37;i++){
 454      sprintf(hname,"QMvsx_%d",i);
 455      fProfileArrayQM->AddAt(new TProfile(hname,hname,50,89,250),i);
 456      sprintf(hname,"QTvsx_%d",i);
 457      fProfileArrayQT->AddAt(new TProfile(hname,hname,50,89,250),i);
 458      sprintf(hname,"QM2D_%d",i);
 459      fProfileArrayQM2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
 460      sprintf(hname,"QT2D_%d",i);
 461      fProfileArrayQT2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
 462    }
 463    //
 464    // just for debugging -counters
 465    //
 466    fTotalTracks     = 0;
 467    fAcceptedTracks  = 0;
 468    fDebugCalPadRaw  = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
 469    fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
 470    // this will be gone for the a new ROOT version > v5-17-05
 471    for (UInt_t i = 0; i < 36; i++) {
 472       fNShortClusters[i]  = 0;
 473       fNMediumClusters[i] = 0;
 474       fNLongClusters[i]   = 0;
 475    }
 476  }
 477
 478 AliTPCcalibTracksGain::~AliTPCcalibTracksGain() {
 479    //
 480    // Destructor.
 481    //
 482
 483    Info("Destructor","");
 484    if (fSimpleFitter) delete fSimpleFitter;
 485    if (fSqrtFitter) delete fSqrtFitter;
 486    if (fLogFitter) delete fLogFitter;
 487    if (fSingleSectorFitter) delete fSingleSectorFitter;
 488
 489    if (fDebugCalPadRaw) delete fDebugCalPadRaw;
 490    if (fDebugCalPadCorr) delete fDebugCalPadCorr;
 491 }
 492
 493 void AliTPCcalibTracksGain::Terminate(){
 494    //
 495    // Evaluate fitters and close the debug stream.
 496    // Also move or copy the debug stream, if a debugStreamPrefix is provided.
 497    //
 498
 499    Evaluate();
 500    AliTPCcalibBase::Terminate();
 501 }
 502
 503
 504
 505 void AliTPCcalibTracksGain::Process(AliTPCseed* seed) {
 506    //
 507    // Main method to be called when a new seed is supposed to be processed
 508    // and be used for gain calibration. Its quality is checked before it
 509    // is added.
 510    //
 511
 512
 513    fTotalTracks++;
 514    if (!fCuts->AcceptTrack(seed)) return;
 515    //
 516    // reinint on proof
 517    //   if (gProof){
 518      static Bool_t doinit= kTRUE;
 519      if (doinit){
 520        fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
 521        fSqrtFitter   = new AliTPCFitPad(8, "hyp7", "");
 522        fLogFitter    = new AliTPCFitPad(8, "hyp7", "");
 523        fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
 524        //
 525        fFitter0M      = new TLinearFitter(45,"hyp44");
 526        fFitter1M      = new TLinearFitter(45,"hyp44");
 527        fFitter2M      = new TLinearFitter(45,"hyp44");
 528        fFitter0T      = new TLinearFitter(45,"hyp44");
 529        fFitter1T      = new TLinearFitter(45,"hyp44");
 530        fFitter2T      = new TLinearFitter(45,"hyp44");
 531        //
 532        fDFitter0M      = new TLinearFitter(10,"hyp9");
 533        fDFitter1M      = new TLinearFitter(10,"hyp9");
 534        fDFitter2M      = new TLinearFitter(10,"hyp9");
 535        fDFitter0T      = new TLinearFitter(10,"hyp9");
 536        fDFitter1T      = new TLinearFitter(10,"hyp9");
 537        fDFitter2T      = new TLinearFitter(10,"hyp9");
 538        doinit=kFALSE;
 539      }
 540      //}
 541
 542    fAcceptedTracks++;
 543    AddTrack(seed);
 544 }
 545
 546 Long64_t AliTPCcalibTracksGain::Merge(TCollection *list) {
 547    //
 548    // Merge() merges the results of all AliTPCcalibTracksGain objects contained in
 549    // list, thus allowing a distributed computation of several files, e.g. on PROOF.
 550    // The merged results are merged with the data members of the AliTPCcalibTracksGain
 551    // object used for calling the Merge method.
 552    // The return value is 0 /*the total number of tracks used for calibration*/ if the merge
 553    // is successful, otherwise it is -1.
 554    //
 555
 556    if (!list || list->IsEmpty()) return -1;
 557
 558    if (!fSimpleFitter) fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
 559    if (!fSqrtFitter)   fSqrtFitter   = new AliTPCFitPad(8, "hyp7", "");
 560    if (!fLogFitter)    fLogFitter    = new AliTPCFitPad(8, "hyp7", "");
 561    if (!fSingleSectorFitter) fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
 562
 563
 564    // just for debugging
 565    if (!fDebugCalPadRaw) fDebugCalPadRaw  = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
 566     if (!fDebugCalPadCorr) fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
 567
 568    TIterator* iter = list->MakeIterator();
 569    AliTPCcalibTracksGain* cal = 0;
 570
 571    while ((cal = (AliTPCcalibTracksGain*)iter->Next())) {
 572       if (!cal->InheritsFrom(AliTPCcalibTracksGain::Class())) {
 573          Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 574          return -1;
 575       }
 576
 577       Add(cal);
 578    }
 579    return 0;
 580 }
 581
 582 void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
 583    //
 584    // Adds another AliTPCcalibTracksGain object to this object.
 585    //
 586
 587    fSimpleFitter->Add(cal->fSimpleFitter);
 588    fSqrtFitter->Add(cal->fSqrtFitter);
 589    fLogFitter->Add(cal->fLogFitter);
 590    fSingleSectorFitter->Add(cal->fSingleSectorFitter);
 591    //
 592    //
 593    //
 594    fFitter0M->Add(cal->fFitter0M);
 595    fFitter1M->Add(cal->fFitter1M);
 596    fFitter2M->Add(cal->fFitter2M);
 597    fFitter0T->Add(cal->fFitter0T);
 598    fFitter1T->Add(cal->fFitter1T);
 599    fFitter2T->Add(cal->fFitter2T);
 600    //
 601    fDFitter0M->Add(cal->fDFitter0M);
 602    fDFitter1M->Add(cal->fDFitter1M);
 603    fDFitter2M->Add(cal->fDFitter2M);
 604    fDFitter0T->Add(cal->fDFitter0T);
 605    fDFitter1T->Add(cal->fDFitter1T);
 606    fDFitter2T->Add(cal->fDFitter2T);
 607    //
 608    //
 609    // histograms
 610    //
 611    for (Int_t i=0; i<73; i++){
 612      TH1F *his,*hism;
 613      his = (TH1F*)fArrayQM->At(i);
 614      hism =  (TH1F*)cal->fArrayQM->At(i);
 615      if (his && hism) his->Add(hism);
 616      his =  (TH1F*)fArrayQT->At(i);
 617      hism =  (TH1F*)cal->fArrayQT->At(i);
 618      if (his && hism) his->Add(hism);
 619    }
 620    //
 621    //
 622    for (Int_t i=0; i<37; i++){
 623      TProfile *his,*hism;
 624      his = (TProfile*)fProfileArrayQM->At(i);
 625      hism = (TProfile*)cal->fProfileArrayQM->At(i);
 626      if (his && hism) his->Add(hism);
 627      his = (TProfile*)fProfileArrayQT->At(i);
 628      hism = (TProfile*)cal->fProfileArrayQT->At(i);
 629      if (his && hism) his->Add(hism);
 630    }
 631    //
 632    //
 633    for (Int_t i=0; i<37; i++){
 634      TProfile2D *his,*hism;
 635      his = (TProfile2D*)fProfileArrayQM2D->At(i);
 636      hism = (TProfile2D*)cal->fProfileArrayQM2D->At(i);
 637      if (his && hism) his->Add(hism);
 638      his = (TProfile2D*)fProfileArrayQT2D->At(i);
 639      hism = (TProfile2D*)cal->fProfileArrayQT2D->At(i);
 640      if (his && hism) his->Add(hism);
 641    }
 642    //
 643    // this will be gone for the a new ROOT version > v5-17-05
 644    for (UInt_t iSegment = 0; iSegment < 36; iSegment++) {
 645       fNShortClusters[iSegment] += cal->fNShortClusters[iSegment];
 646       fNMediumClusters[iSegment] += cal->fNMediumClusters[iSegment];
 647       fNLongClusters[iSegment] += cal->fNLongClusters[iSegment];
 648    }
 649
 650    // just for debugging, remove me
 651    fTotalTracks += cal->fTotalTracks;
 652    fAcceptedTracks += cal->fAcceptedTracks;
 653    fDebugCalPadRaw->Add(cal->fDebugCalPadRaw);
 654    fDebugCalPadCorr->Add(cal->fDebugCalPadCorr);
 655
 656 }
 657
 658 void AliTPCcalibTracksGain::AddTrack(AliTPCseed* seed) {
 659    //
 660    // The clusters making up the track (seed) are added to various fit functions.
 661    // See AddCluster(...) for more detail.
 662    //
 663
 664    DumpTrack(seed);
 665    //
 666    // simple histograming part
 667    for (Int_t i=0; i<159; i++){
 668      AliTPCclusterMI* cluster = seed->GetClusterPointer(i);
 669      if (cluster) AddCluster(cluster);
 670    }
 671 }
 672
 673 void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster){
 674   //
 675   // Adding cluster information to the simple histograms
 676   // No correction, fittings are applied
 677   //
 678   Float_t kThreshold=5;
 679   if (cluster->GetX()<=0) return;
 680   if (cluster->GetQ()<=kThreshold) return;
 681   //
 682
 683
 684   Int_t sector = cluster->GetDetector();
 685   TH1F  * his;
 686   his = GetQT(sector);
 687   if (his) his->Fill(cluster->GetQ());
 688   his = GetQT(-1);
 689   if (his) his->Fill(cluster->GetQ());
 690   his = GetQM(sector);
 691   if (his) his->Fill(cluster->GetMax());
 692   his = GetQM(-1);
 693   if (his) his->Fill(cluster->GetMax());
 694   //
 695   sector = sector%36;
 696   TProfile *prof;
 697   prof = GetProfileQT(sector);
 698   if (prof) prof->Fill(cluster->GetX(),cluster->GetQ());
 699   prof = GetProfileQT(-1);
 700   if (prof) prof->Fill(cluster->GetX(),cluster->GetQ());
 701   prof = GetProfileQM(sector);
 702   if (prof) prof->Fill(cluster->GetX(),cluster->GetMax());
 703   prof = GetProfileQM(-1);
 704   if (prof) prof->Fill(cluster->GetX(),cluster->GetMax());
 705   //
 706   Float_t phi = cluster->GetY()/cluster->GetX();
 707   TProfile2D *prof2;
 708   prof2 = GetProfileQT2D(sector);
 709   if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetQ());
 710   prof2 = GetProfileQT2D(-1);
 711   if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetQ());
 712   prof2 = GetProfileQM2D(sector);
 713   if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetMax());
 714   prof2 = GetProfileQM2D(-1);
 715   if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetMax());
 716
 717   //
 718 }
 719
 720
 721
 722 void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t /*momenta*/, Float_t/* mdedx*/, Int_t padType,
 723                                        Float_t xcenter, TVectorD& dedxQ, TVectorD& /*dedxM*/, Float_t /*fraction*/, Float_t fraction2, Float_t dedge,
 724                                        TVectorD& /*parY*/, TVectorD& /*parZ*/, TVectorD& meanPos) {
 725    //
 726    // Adds cluster to the appropriate fitter for later analysis.
 727    // The charge used for the fit is the maximum charge for this specific cluster or the
 728    // accumulated charge per cluster, depending on the value of fgkUseTotalCharge.
 729    // Depending on the pad size where the cluster is registered, the value will be put in
 730    // the appropriate fitter. Furthermore, for each pad size three different types of fitters
 731    // are used. The fit functions are the same for all fitters (parabolic functions), but the value
 732    // added to each fitter is different. The simple fitter gets the charge plugged in as is, the sqrt fitter
 733    // gets the square root of the charge, and the log fitter gets fgkM*(1+q/fgkM), where q is the original charge
 734    // and fgkM==25.
 735    //
 736
 737    if (!cluster) {
 738       Error("AddCluster", "Cluster not valid.");
 739       return;
 740    }
 741
 742    if (dedge < 3.) return;
 743    if (fraction2 > 0.9) return;
 744
 745    //Int_t padType = GetPadType(cluster->GetX());
 746    Double_t xx[7];
 747    //Double_t centerPad[2] = {0};
 748    //AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
 749    //xx[0] = cluster->GetX() - centerPad[0];
 750    //xx[1] = cluster->GetY() - centerPad[1];
 751    xx[0] = cluster->GetX() - xcenter;
 752    xx[1] = cluster->GetY();
 753    xx[2] = xx[0] * xx[0];
 754    xx[3] = xx[1] * xx[1];
 755    xx[4] = xx[0] * xx[1];
 756    xx[5] = TMath::Abs(cluster->GetZ()) - TMath::Abs(meanPos[4]);
 757    xx[6] = xx[5] * xx[5];
 758    //
 759    // Update profile histograms
 760    //
 761
 762    //
 763    // Update fitters
 764    //
 765    Int_t segment = cluster->GetDetector() % 36;
 766    Double_t q = fgkUseTotalCharge ? ((Double_t)(cluster->GetQ())) : ((Double_t)(cluster->GetMax()));  // note: no normalization to pad size!
 767
 768    // just for debugging
 769    Int_t row = 0;
 770    Int_t pad = 0;
 771    GetRowPad(cluster->GetX(), cluster->GetY(), row, pad);
 772    fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
 773
 774    // correct charge by normalising to mean charge per track
 775    q /= dedxQ[2];
 776
 777    // just for debugging
 778    fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
 779
 780    Double_t sqrtQ = TMath::Sqrt(q);
 781    Double_t logQ = fgkM * TMath::Log(1 + q / fgkM);
 782    TLinearFitter * fitter =0;
 783    //
 784    fitter = fSimpleFitter->GetFitter(segment, padType);
 785    fitter->AddPoint(xx, q);
 786    //
 787    fitter = fSqrtFitter->GetFitter(segment, padType);
 788    fitter->AddPoint(xx, sqrtQ);
 789    //
 790    fitter = fLogFitter->GetFitter(segment, padType);
 791    fitter->AddPoint(xx, logQ);
 792    //
 793    fitter=fSingleSectorFitter->GetFitter(0, padType);
 794    fitter->AddPoint(xx, q);
 795
 796    // this will be gone for the a new ROOT version > v5-17-05
 797    if (padType == kShortPads)
 798       fNShortClusters[segment]++;
 799    if (padType == kMediumPads)
 800       fNMediumClusters[segment]++;
 801    if (padType == kLongPads)
 802       fNLongClusters[segment]++;
 803 }
 804
 805 void AliTPCcalibTracksGain::Evaluate(Bool_t robust, Double_t frac) {
 806    //
 807    // Evaluates all fitters contained in this object.
 808    // If the robust option is set to kTRUE a robust fit is performed with frac as
 809    // the minimal fraction of good points (see TLinearFitter::EvalRobust for details).
 810    // Beware: Robust fitting is much slower!
 811    //
 812
 813    fSimpleFitter->Evaluate(robust, frac);
 814    fSqrtFitter->Evaluate(robust, frac);
 815    fLogFitter->Evaluate(robust, frac);
 816    fSingleSectorFitter->Evaluate(robust, frac);
 817    fFitter0M->Eval();
 818    fFitter1M->Eval();
 819    fFitter2M->Eval();
 820    fFitter0T->Eval();
 821    fFitter1T->Eval();
 822    fFitter2T->Eval();
 823    //
 824    fDFitter0M->Eval();
 825    fDFitter1M->Eval();
 826    fDFitter2M->Eval();
 827    fDFitter0T->Eval();
 828    fDFitter1T->Eval();
 829    fDFitter2T->Eval();
 830 }
 831
 832 AliTPCCalPad* AliTPCcalibTracksGain::CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 833   //
 834   // Creates the calibration object AliTPCcalPad using fitted parameterization
 835   //
 836    TObjArray tpc(72);
 837    for (UInt_t iSector = 0; iSector < 72; iSector++)
 838       tpc.Add(CreateFitCalROC(iSector, fitType, undoTransformation, normalizeToPadSize));
 839    return new AliTPCCalPad(&tpc);
 840 }
 841
 842 AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 843   //
 844   // Create the AliTPCCalROC with the values per pad
 845   // sector  - sector of interest
 846   // fitType -
 847   //
 848
 849    TVectorD par(8);
 850    if (sector < 36) {
 851       GetParameters(sector % 36, 0, fitType, par);
 852       return CreateFitCalROC(sector, 0, par, fitType, undoTransformation, normalizeToPadSize);
 853    }
 854    else {
 855       GetParameters(sector % 36, 1, fitType, par);
 856       AliTPCCalROC* roc1 = CreateFitCalROC(sector, 1, par, fitType, undoTransformation, normalizeToPadSize);
 857       GetParameters(sector % 36, 2, fitType, par);
 858       AliTPCCalROC* roc2 = CreateFitCalROC(sector, 2, par, fitType, undoTransformation, normalizeToPadSize);
 859       AliTPCCalROC* roc3 = CreateCombinedCalROC(roc1, roc2);
 860       delete roc1;
 861       delete roc2;
 862       return roc3;
 863    }
 864 }
 865
 866 AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t padType, TVectorD &fitParam, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 867    //
 868    // This function is essentially a copy of AliTPCCalROC::CreateGlobalFitCalROC(...), with the
 869    // modifications, that the center of the region of same pad size is used as the origin
 870    // of the fit function instead of the center of the ROC.
 871    // The possibility of a linear fit is removed as well because it is not needed.
 872    // Only values for pads with the given pad size are calculated, the rest is 0.
 873    // Set undoTransformation for undoing the transformation that was applied to the
 874    // charge values before they were put into the fitter (thus allowing comparison to the original
 875    // charge values). For fitType use 0 for the simple fitter, 1 for the sqrt fitter, 2 for the log fitter.
 876    // If normalizeToPadSize is true, the values are normalized to the pad size.
 877    // Please be aware, that you even need to specify the fitType if you want to normalize to the pad size without
 878    // undoing the transformation (because normalizing involves undoing the trafo first, then normalizing, then
 879    // applying the trafo again).
 880    // Please note: The normalization to the pad size is a simple linear scaling with the pad length, which
 881    //              actually doesn't describe reality!
 882    //
 883
 884    Float_t dlx, dly;
 885    Double_t centerPad[2] = {0};
 886    Float_t localXY[3] = {0};
 887    AliTPCROC* tpcROC = AliTPCROC::Instance();
 888    if ((padType == 0 && sector >= tpcROC->GetNInnerSector()) || (padType > 0 && sector < tpcROC->GetNInnerSector()) || sector >= tpcROC->GetNSector())
 889       return 0;
 890    AliTPCCalROC* lROCfitted = new AliTPCCalROC(sector);
 891    //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
 892    UInt_t startRow = 0;
 893    UInt_t endRow = 0;
 894    switch (padType) {
 895       case kShortPads:
 896          startRow = 0;
 897          endRow = lROCfitted->GetNrows();
 898          break;
 899       case kMediumPads:
 900          startRow = 0;
 901          endRow = 64;
 902          break;
 903       case kLongPads:
 904          startRow = 64;
 905          endRow = lROCfitted->GetNrows();
 906          break;
 907    }
 908
 909    AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
 910    Double_t value = 0;
 911    for (UInt_t irow = startRow; irow < endRow; irow++) {
 912       for (UInt_t ipad = 0; ipad < lROCfitted->GetNPads(irow); ipad++) {
 913          tpcROC->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 914          dlx = localXY[0] - centerPad[0];
 915          dly = localXY[1] - centerPad[1];
 916          value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
 917
 918          // Let q' = value be the transformed value without any pad size corrections,
 919          // let T be the transformation and let l be the pad size
 920          //    1) don't undo transformation, don't normalize: return q'
 921          //    2) undo transformation,       don't normalize: return T^{-1} q'
 922          //    3) undo transformation,       normalize:       return (T^{-1} q') / l
 923          //    4) don't undo transformation, normalize:       return T((T^{-1} q') / l)
 924          if (!undoTransformation && !normalizeToPadSize) {/* value remains unchanged */}  // (1)
 925          else {                                                                           // (2), (3), (4)
 926             //calculate T^{-1}
 927             switch (fitType) {
 928                case  0: /* value remains unchanged */ break;
 929                case  1: value = value * value; break;
 930                case  2: value = (TMath::Exp(value / fgkM) - 1) * fgkM; break;
 931                default: Error("CreateFitCalROC", "Wrong fit type."); break;
 932             }
 933             if (normalizeToPadSize) value /= GetPadLength(localXY[0]);                    // (3)
 934          }
 935          if (!undoTransformation && normalizeToPadSize) {                                 // (4)
 936             // calculate T
 937             switch (fitType) {
 938                case  0: /* value remains unchanged */ break;
 939                case  1: value = TMath::Sqrt(value); break;
 940                case  2: value = fgkM * TMath::Log(1 + value / fgkM); break;
 941                default: Error("CreateFitCalROC", "Wrong fit type."); break;
 942             }
 943          }
 944          lROCfitted->SetValue(irow, ipad, value);
 945       }
 946    }
 947    return lROCfitted;
 948 }
 949
 950 AliTPCCalROC* AliTPCcalibTracksGain::CreateCombinedCalROC(const AliTPCCalROC* roc1, const AliTPCCalROC* roc2) {
 951    //
 952    // Combines the medium pad size values of roc1 with the long pad size values of roc2 into a new
 953    // AliTPCCalROC. Returns a null pointer if any one of the ROCs is an IROC; issues a warning message
 954    // if the sectors of roc1 and roc2 don't match, but still continue and use the sector of roc1 as the
 955    // sector of the new ROC.
 956    //
 957
 958    if (!roc1 || !roc2) return 0;
 959    if ((Int_t)(roc1->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
 960    if ((Int_t)(roc2->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
 961    if (roc1->GetSector() != roc2->GetSector()) Warning("CreateCombinedCalROC", "Sector number mismatch.");
 962    AliTPCCalROC* roc = new AliTPCCalROC(roc1->GetSector());
 963
 964    for (UInt_t iRow = 0; iRow < 64; iRow++) {
 965       for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
 966          roc->SetValue(iRow, iPad, roc1->GetValue(iRow, iPad));
 967    }
 968    for (UInt_t iRow = 64; iRow < roc->GetNrows(); iRow++) {
 969       for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
 970          roc->SetValue(iRow, iPad, roc2->GetValue(iRow, iPad));
 971    }
 972    return roc;
 973 }
 974
 975 Bool_t AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) {
 976    //
 977    // Puts the fit parameters for the specified segment (IROC & OROC), padType and fitType
 978    // into the fitParam TVectorD (which should contain 8 elements).
 979    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 980    // Note: The fitter has to be evaluated first!
 981    //
 982   TLinearFitter * fitter = GetFitter(segment, padType, fitType);
 983   if (fitter){
 984     fitter->Eval();
 985     fitter->GetParameters(fitParam);
 986     return kTRUE;
 987   }else{
 988     Error("AliTPCcalibTracksGain::GetParameters",
 989           Form("Fitter%d_%d_%d not availble", segment, padType, fitType));
 990     return kFALSE;
 991   }
 992   return kFALSE;
 993 }
 994
 995 void AliTPCcalibTracksGain::GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError) {
 996    //
 997    // Puts the fit parameter errors for the specified segment (IROC & OROC), padType and fitType
 998    // into the fitError TVectorD (which should contain 8 elements).
 999    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1000    // Note: The fitter has to be evaluated first!
1001    //
1002
1003    GetFitter(segment, padType, fitType)->GetErrors(fitError);
1004    fitError *= TMath::Sqrt(GetRedChi2(segment, padType, fitType));
1005 }
1006
1007 Double_t AliTPCcalibTracksGain::GetRedChi2(UInt_t segment, UInt_t padType, UInt_t fitType) {
1008    //
1009    // Returns the reduced chi^2 value for the specified segment, padType and fitType.
1010    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1011    // Note: The fitter has to be evaluated first!
1012    //
1013
1014    // this will be gone for the a new ROOT version > v5-17-05
1015    Int_t lNClusters = 0;
1016    switch (padType) {
1017       case kShortPads:
1018          lNClusters = fNShortClusters[segment];
1019          break;
1020       case kMediumPads:
1021          lNClusters = fNMediumClusters[segment];
1022          break;
1023       case kLongPads:
1024          lNClusters = fNLongClusters[segment];
1025          break;
1026    }
1027    return GetFitter(segment, padType, fitType)->GetChisquare()/(lNClusters - 8);
1028 }
1029
1030 void AliTPCcalibTracksGain::GetCovarianceMatrix(UInt_t segment, UInt_t padType, UInt_t fitType, TMatrixD& covMatrix) {
1031    //
1032    // Returns the covariance matrix for the specified segment, padType, fitType.
1033    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1034    //
1035
1036    GetFitter(segment, padType, fitType)->GetCovarianceMatrix(covMatrix);
1037 }
1038
1039 TLinearFitter* AliTPCcalibTracksGain::GetFitter(UInt_t segment, UInt_t padType, UInt_t fitType) {
1040    //
1041    // Returns the TLinearFitter object for the specified segment, padType, fitType.
1042    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1043    //
1044
1045    switch (fitType) {
1046       case kSimpleFitter:
1047          return fSimpleFitter->GetFitter(segment, padType);
1048       case kSqrtFitter:
1049          return fSqrtFitter->GetFitter(segment, padType);
1050       case kLogFitter:
1051          return fLogFitter->GetFitter(segment, padType);
1052       case 3:
1053          return fSingleSectorFitter->GetFitter(0, padType);
1054    }
1055    return 0;
1056 }
1057
1058 Double_t AliTPCcalibTracksGain::GetPadLength(Double_t lx) {
1059    //
1060    // The function returns 0.75 for an IROC, 1. for an OROC at medium pad size position,
1061    // 1.5 for an OROC at long pad size position, -1 if out of bounds.
1062    //
1063
1064    Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
1065    Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
1066    Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
1067    Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
1068    Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
1069    Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
1070
1071    // if IROC
1072    if (lx >= irocLow && lx <= irocUp) return 0.75;
1073    // if OROC medium pads
1074    if (lx >= orocLow1 && lx <= orocUp1) return 1.;
1075    // if OROC long pads
1076    if (lx >= orocLow2 && lx <= orocUp2) return 1.5;
1077    // if out of bounds
1078    return -1;
1079 }
1080
1081 Int_t AliTPCcalibTracksGain::GetPadType(Double_t lx) {
1082    //
1083    // The function returns 0 for an IROC, 1 for an OROC at medium pad size position,
1084    // 2 for an OROC at long pad size position, -1 if out of bounds.
1085    //
1086
1087    if (GetPadLength(lx) == 0.75) return 0;
1088    else if (GetPadLength(lx) == 1.) return 1;
1089    else if (GetPadLength(lx) == 1.5) return 2;
1090    return -1;
1091 }
1092
1093 // ONLY FOR DEBUGGING PURPOSES - REMOVE ME WHEN NOT NEEDED ANYMORE
1094 Bool_t AliTPCcalibTracksGain::GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad) {
1095    //
1096    // Calculate the row and pad number when the local coordinates are given.
1097    // Returns kFALSE if the position is out of range, otherwise return kTRUE.
1098    // WARNING: This function is preliminary and probably isn't very accurate!!
1099    //
1100
1101    Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
1102    //Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
1103    Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
1104    //Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
1105    Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
1106    //Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
1107
1108    if (GetPadType(lx) == 0) {
1109       row = (Int_t)((lx - irocLow) / fgTPCparam->GetInnerPadPitchLength());
1110       pad = (Int_t)((ly + fgTPCparam->GetYInner(row)) / fgTPCparam->GetInnerPadPitchWidth());
1111    } else if (GetPadType(lx) == 1) {
1112       row = (Int_t)((lx - orocLow1) / fgTPCparam->GetOuter1PadPitchLength());
1113       pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
1114    } else if (GetPadType(lx) == 2) {
1115       row = fgTPCparam->GetNRowUp1() + (Int_t)((lx - orocLow2) / fgTPCparam->GetOuter2PadPitchLength());
1116       pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
1117    }
1118    else return kFALSE;
1119    return kTRUE;
1120 }
1121
1122 void AliTPCcalibTracksGain::DumpTrack(AliTPCseed* track) {
1123    //
1124    //  Dump track information to the debug stream
1125    //
1126
1127    Int_t rows[200];
1128    Int_t sector[3];
1129    Int_t npoints[3];
1130    TVectorD dedxM[3];
1131    TVectorD dedxQ[3];
1132    TVectorD parY[3];
1133    TVectorD parZ[3];
1134    TVectorD meanPos[3];
1135    //
1136    Int_t count=0;
1137    for (Int_t ipad = 0; ipad < 3; ipad++) {
1138      dedxM[ipad].ResizeTo(5);
1139      dedxQ[ipad].ResizeTo(5);
1140      parY[ipad].ResizeTo(3);
1141      parZ[ipad].ResizeTo(3);
1142      meanPos[ipad].ResizeTo(6);
1143      Bool_t isOK = GetDedx(track, ipad, rows, sector[ipad], npoints[ipad], dedxM[ipad], dedxQ[ipad], parY[ipad], parZ[ipad], meanPos[ipad]);
1144      if (isOK)
1145        AddTracklet(sector[ipad],ipad, dedxQ[ipad], dedxM[ipad], parY[ipad], parZ[ipad], meanPos[ipad] );
1146      if (isOK) count++;
1147    }
1148
1149    TTreeSRedirector * cstream =  GetDebugStreamer();
1150    if (cstream){
1151      (*cstream) << "Track" <<
1152        "Track.=" << track <<        // track information
1153        "\n";
1154      //
1155      //
1156      //
1157      if ( GetStreamLevel()>1 && count>1){
1158        (*cstream) << "TrackG" <<
1159          "Track.=" << track <<        // track information
1160          "ncount="<<count<<
1161          // info for pad type 0
1162          "sector0="<<sector[0]<<
1163          "npoints0="<<npoints[0]<<
1164          "dedxM0.="<<&dedxM[0]<<
1165          "dedxQ0.="<<&dedxQ[0]<<
1166          "parY0.="<<&parY[0]<<
1167          "parZ0.="<<&parZ[0]<<
1168          "meanPos0.="<<&meanPos[0]<<
1169          //
1170          // info for pad type 1
1171          "sector1="<<sector[1]<<
1172          "npoints1="<<npoints[1]<<
1173          "dedxM1.="<<&dedxM[1]<<
1174          "dedxQ1.="<<&dedxQ[1]<<
1175          "parY1.="<<&parY[1]<<
1176          "parZ1.="<<&parZ[1]<<
1177          "meanPos1.="<<&meanPos[1]<<
1178          //
1179          // info for pad type 2
1180          "sector2="<<sector[2]<<
1181          "npoints2="<<npoints[2]<<
1182          "dedxM2.="<<&dedxM[2]<<
1183          "dedxQ2.="<<&dedxQ[2]<<
1184          "parY2.="<<&parY[2]<<
1185          "parZ2.="<<&parZ[2]<<
1186          "meanPos2.="<<&meanPos[2]<<
1187          //
1188          "\n";
1189      }
1190    }
1191 }
1192
1193 Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /*rows*/,
1194                                       Int_t &sector, Int_t& npoints,
1195                                       TVectorD &dedxM, TVectorD &dedxQ,
1196                                       TVectorD &parY, TVectorD &parZ, TVectorD&meanPos)
1197 {
1198   //
1199   // GetDedx for given sector for given track
1200   // padType - type of pads
1201   //
1202
1203   static TLinearFitter fitY(2, "pol1");
1204   static TLinearFitter fitZ(2, "pol1");
1205   fitY.StoreData(kFALSE);
1206   fitZ.StoreData(kFALSE);
1207   //
1208   //
1209    Int_t firstRow = 0, lastRow = 0;
1210    Int_t minRow = 100;
1211    Float_t xcenter = 0;
1212    const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
1213    const Float_t kedgey = 4.;
1214    if (padType == 0) {
1215       firstRow = 0;
1216       lastRow = fgTPCparam->GetNRowLow();
1217       xcenter = 108.47;
1218    }
1219    if (padType == 1) {
1220       firstRow = fgTPCparam->GetNRowLow();
1221       lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1222       xcenter = 166.60;
1223    }
1224    if (padType == 2) {
1225       firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1226       lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
1227       xcenter = 222.6;
1228    }
1229    minRow = (lastRow - firstRow) / 2;
1230    //
1231    //
1232    Int_t nclusters = 0;
1233    Int_t nclustersNE = 0; // number of not edge clusters
1234    Int_t lastSector = -1;
1235    Float_t amplitudeQ[100];
1236    Float_t amplitudeM[100];
1237    Int_t rowIn[100];
1238    Int_t index[100];
1239    //
1240    //
1241    fitY.ClearPoints();
1242    fitZ.ClearPoints();
1243
1244    for (Int_t iCluster = firstRow; iCluster < lastRow; iCluster++) {
1245       AliTPCclusterMI* cluster = track->GetClusterPointer(iCluster);
1246       if (cluster) {
1247          Int_t detector = cluster->GetDetector() ;
1248          if (lastSector == -1) lastSector = detector;
1249          if (lastSector != detector) continue;
1250          amplitudeQ[nclusters] = cluster->GetQ();
1251          amplitudeM[nclusters] = cluster->GetMax();
1252          rowIn[nclusters] = iCluster;
1253          nclusters++;
1254          Double_t dx = cluster->GetX() - xcenter;
1255          Double_t y = cluster->GetY();
1256          Double_t z = cluster->GetZ();
1257          fitY.AddPoint(&dx, y);
1258          fitZ.AddPoint(&dx, z);
1259          meanPos[0] += dx;
1260          meanPos[1] += dx;
1261          meanPos[2] += y;
1262          meanPos[3] += y*y;
1263          meanPos[4] += z;
1264          meanPos[5] += z*z;
1265          if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) nclustersNE++;
1266       }
1267    }
1268
1269    if (nclusters < minRow / 2) return kFALSE;
1270    if (nclustersNE < minRow / 2) return kFALSE;
1271    for (Int_t i = 0; i < 6; i++) meanPos[i] /= Double_t(nclusters);
1272    fitY.Eval();
1273    fitZ.Eval();
1274    fitY.GetParameters(parY);
1275    fitZ.GetParameters(parZ);
1276    //
1277    // calculate truncated mean
1278    //
1279    TMath::Sort(nclusters, amplitudeQ, index, kFALSE);
1280    //
1281    //
1282    //
1283    Float_t ndedx[5];
1284    for (Int_t i = 0; i < 5; i++) {
1285       dedxQ[i] = 0;
1286       dedxM[i] = 0;
1287       ndedx[i] = 0;
1288    }
1289    //
1290    // dedx calculation
1291    //
1292    Int_t inonEdge = 0;
1293    for (Int_t i = 0; i < nclusters; i++) {
1294       Int_t rowSorted = rowIn[index[i]];
1295       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1296
1297       if (TMath::Abs(cluster->GetY()) > cluster->GetX()*ktany - kedgey) continue;  //don't take edge clusters
1298       inonEdge++;
1299       if (inonEdge < nclustersNE * 0.5) {
1300          ndedx[0]++;
1301          dedxQ[0] += amplitudeQ[index[i]];
1302          dedxM[0] += amplitudeM[index[i]];
1303       }
1304       if (inonEdge < nclustersNE * 0.6) {
1305          ndedx[1]++;
1306          dedxQ[1] += amplitudeQ[index[i]];
1307          dedxM[1] += amplitudeM[index[i]];
1308       }
1309       if (inonEdge < nclustersNE * 0.7) {
1310          ndedx[2]++;
1311          dedxQ[2] += amplitudeQ[index[i]];
1312          dedxM[2] += amplitudeM[index[i]];
1313       }
1314       if (inonEdge < nclustersNE * 0.8) {
1315          ndedx[3]++;
1316          dedxQ[3] += amplitudeQ[index[i]];
1317          dedxM[3] += amplitudeM[index[i]];
1318       }
1319       if (inonEdge < nclustersNE * 0.9) {
1320          ndedx[4]++;
1321          dedxQ[4] += amplitudeQ[index[i]];
1322          dedxM[4] += amplitudeM[index[i]];
1323       }
1324    }
1325    for (Int_t i = 0; i < 5; i++) {
1326       dedxQ[i] /= ndedx[i];
1327       dedxM[i] /= ndedx[i];
1328    }
1329    TTreeSRedirector * cstream =  GetDebugStreamer();
1330    inonEdge = 0;
1331    Float_t momenta = track->GetP();
1332    Float_t mdedx = track->GetdEdx();
1333    for (Int_t i = 0; i < nclusters; i++) {
1334       Int_t rowSorted = rowIn[index[i]];
1335       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1336       if (!cluster) {
1337          printf("Problem\n");
1338          continue;
1339       }
1340       if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) inonEdge++;
1341       Float_t dedge = cluster->GetX()*ktany - TMath::Abs(cluster->GetY());
1342       Float_t fraction = Float_t(i) / Float_t(nclusters);
1343       Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
1344
1345       AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
1346
1347       if (cstream) (*cstream) << "dEdx" <<
1348          "Cl.=" << cluster <<           // cluster of interest
1349          "P=" << momenta <<             // track momenta
1350          "dedx=" << mdedx <<            // mean dedx - corrected for angle
1351          "IPad=" << padType <<          // pad type 0..2
1352          "xc=" << xcenter <<            // x center of chamber
1353          "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
1354          "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
1355          "fraction=" << fraction <<     // fraction - order in statistic (0,1)
1356          "fraction2=" << fraction2 <<   // fraction - order in statistic (0,1)
1357          "dedge=" << dedge <<           // distance to the edge
1358          "parY.=" << &parY <<           // line fit
1359          "parZ.=" << &parZ <<           // line fit
1360          "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
1361          "\n";
1362    }
1363
1364    if (cstream) (*cstream) << "dEdxT" <<
1365      "P=" << momenta <<             // track momenta
1366      "npoints="<<inonEdge<<         // number of points
1367      "sector="<<lastSector<<        // sector number
1368      "dedx=" << mdedx <<            // mean dedx - corrected for angle
1369      "IPad=" << padType <<          // pad type 0..2
1370      "xc=" << xcenter <<            // x center of chamber
1371      "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
1372      "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
1373      "parY.=" << &parY <<           // line fit
1374      "parZ.=" << &parZ <<           // line fit
1375      "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
1376      "\n";
1377
1378    sector = lastSector;
1379    npoints = inonEdge;
1380    return kTRUE;
1381 }
1382
1383 void AliTPCcalibTracksGain::AddTracklet(UInt_t sector, UInt_t padType,TVectorD &dedxQ, TVectorD &dedxM,TVectorD& parY, TVectorD& parZ, TVectorD& meanPos){
1384   //
1385   // Add measured point - dedx to the fitter
1386   //
1387   //
1388   //chain->SetAlias("dr","(250-abs(meanPos.fElements[4]))/250");
1389   //chain->SetAlias("tz","(0+abs(parZ.fElements[1]))");
1390   //chain->SetAlias("ty","(0+abs(parY.fElements[1]))");
1391   //chain->SetAlias("corrg","sqrt((1+ty^2)*(1+tz^2))");
1392   //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);
1393
1394   Double_t xxx[100];
1395   //
1396   // z and angular part
1397   //
1398
1399   xxx[0] = (250.-TMath::Abs(meanPos[4]))/250.;
1400   xxx[1] = TMath::Abs(parY[1]);
1401   xxx[2] = TMath::Abs(parZ[1]);
1402   xxx[3] = xxx[0]*xxx[1];
1403   xxx[4] = xxx[0]*xxx[2];
1404   xxx[5] = xxx[1]*xxx[2];
1405   xxx[6] = xxx[0]*xxx[0];
1406   xxx[7] = xxx[1]*xxx[1];
1407   xxx[8] = xxx[2]*xxx[2];
1408   //
1409   // chamber part
1410   //
1411   Int_t tsector = sector%36;
1412   for (Int_t i=0;i<35;i++){
1413     xxx[9+i]=(i==tsector)?1:0;
1414   }
1415   TLinearFitter *fitterM = fFitter0M;
1416   if (padType==1) fitterM=fFitter1M;
1417   if (padType==2) fitterM=fFitter2M;
1418   fitterM->AddPoint(xxx,dedxM[1]);
1419   //
1420   TLinearFitter *fitterT = fFitter0T;
1421   if (padType==1) fitterT = fFitter1T;
1422   if (padType==2) fitterT = fFitter2T;
1423   fitterT->AddPoint(xxx,dedxQ[1]);
1424   //
1425   TLinearFitter *dfitterM = fDFitter0M;
1426   if (padType==1) dfitterM=fDFitter1M;
1427   if (padType==2) dfitterM=fDFitter2M;
1428   dfitterM->AddPoint(xxx,dedxM[1]);
1429   //
1430   TLinearFitter *dfitterT = fDFitter0T;
1431   if (padType==1) dfitterT = fDFitter1T;
1432   if (padType==2) dfitterT = fDFitter2T;
1433   dfitterT->AddPoint(xxx,dedxQ[1]);
1434 }
1435
1436
1437 TGraph *AliTPCcalibTracksGain::CreateAmpGraph(Int_t ipad, Bool_t qmax){
1438   //
1439   // create the amplitude graph
1440   // The normalized amplitudes are extrapolated to the 0 angle (y,z)  and 0 drift length
1441   //
1442
1443   TVectorD vec;
1444   if (qmax){
1445     if (ipad==0) fFitter0M->GetParameters(vec);
1446     if (ipad==1) fFitter1M->GetParameters(vec);
1447     if (ipad==2) fFitter2M->GetParameters(vec);
1448   }else{
1449     if (ipad==0) fFitter0T->GetParameters(vec);
1450     if (ipad==1) fFitter1T->GetParameters(vec);
1451     if (ipad==2) fFitter2T->GetParameters(vec);
1452   }
1453
1454   Float_t amp[36];
1455   Float_t sec[36];
1456   for (Int_t i=0;i<35;i++){
1457     sec[i]=i;
1458     amp[i]=vec[10+i]+vec[0];
1459   }
1460   amp[35]=vec[0];
1461   Float_t mean = TMath::Mean(36,amp);
1462   for (Int_t i=0;i<36;i++){
1463     sec[i]=i;
1464     amp[i]=(amp[i]-mean)/mean;
1465   }
1466   TGraph *gr = new TGraph(36,sec,amp);
1467   return gr;
1468 }
1469
1470
1471 void   AliTPCcalibTracksGain::UpdateClusterParam(AliTPCClusterParam* clparam){
1472   //
1473   //   SetQ normalization parameters
1474   //
1475   //  void SetQnorm(Int_t ipad, Int_t itype,  TVectorD * norm);
1476
1477   TVectorD vec;
1478
1479   //
1480   fDFitter0T->Eval();
1481   fDFitter1T->Eval();
1482   fDFitter2T->Eval();
1483   fDFitter0M->Eval();
1484   fDFitter1M->Eval();
1485   fDFitter2M->Eval();
1486   fDFitter0T->GetParameters(vec);
1487   clparam->SetQnorm(0,0,&vec);
1488   fDFitter1T->GetParameters(vec);
1489   clparam->SetQnorm(1,0,&vec);
1490   fDFitter2T->GetParameters(vec);
1491   clparam->SetQnorm(2,0,&vec);
1492   //
1493   fDFitter0M->GetParameters(vec);
1494   clparam->SetQnorm(0,1,&vec);
1495   fDFitter1M->GetParameters(vec);
1496   clparam->SetQnorm(1,1,&vec);
1497   fDFitter2M->GetParameters(vec);
1498   clparam->SetQnorm(2,1,&vec);
1499   //
1500
1501 }
1502
1503
1504 void   AliTPCcalibTracksGain::Analyze(){
1505
1506  Evaluate();
1507
1508 }
1509
1510