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