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