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) :
 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        "run="<<fRun<<              //  run number
1223        "event="<<fEvent<<          //  event number
1224        "time="<<fTime<<            //  time stamp of event
1225        "trigger="<<fTrigger<<      //  trigger
1226        "mag="<<fMagF<<             //  magnetic field
1227        "Track.=" << track <<        // track information
1228        "\n";
1229      //
1230      //
1231      //
1232      if ( GetStreamLevel()>1 && count>1){
1233        (*cstream) << "TrackG" <<
1234          "run="<<fRun<<              //  run number
1235          "event="<<fEvent<<          //  event number
1236          "time="<<fTime<<            //  time stamp of event
1237          "trigger="<<fTrigger<<      //  trigger
1238          "mag="<<fMagF<<             //  magnetic field
1239          "Track.=" << track <<        // track information
1240          "ncount="<<count<<
1241          // info for pad type 0
1242          "sector0="<<sector[0]<<
1243          "npoints0="<<npoints[0]<<
1244          "dedxM0.="<<&dedxM[0]<<
1245          "dedxQ0.="<<&dedxQ[0]<<
1246          "parY0.="<<&parY[0]<<
1247          "parZ0.="<<&parZ[0]<<
1248          "meanPos0.="<<&meanPos[0]<<
1249          //
1250          // info for pad type 1
1251          "sector1="<<sector[1]<<
1252          "npoints1="<<npoints[1]<<
1253          "dedxM1.="<<&dedxM[1]<<
1254          "dedxQ1.="<<&dedxQ[1]<<
1255          "parY1.="<<&parY[1]<<
1256          "parZ1.="<<&parZ[1]<<
1257          "meanPos1.="<<&meanPos[1]<<
1258          //
1259          // info for pad type 2
1260          "sector2="<<sector[2]<<
1261          "npoints2="<<npoints[2]<<
1262          "dedxM2.="<<&dedxM[2]<<
1263          "dedxQ2.="<<&dedxQ[2]<<
1264          "parY2.="<<&parY[2]<<
1265          "parZ2.="<<&parZ[2]<<
1266          "meanPos2.="<<&meanPos[2]<<
1267          //
1268          "\n";
1269      }
1270    }
1271 }
1272
1273 Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /*rows*/,
1274                                       Int_t &sector, Int_t& npoints,
1275                                       TVectorD &dedxM, TVectorD &dedxQ,
1276                                       TVectorD &parY, TVectorD &parZ, TVectorD&meanPos)
1277 {
1278   //
1279   // GetDedx for given sector for given track
1280   // padType - type of pads
1281   //
1282
1283   static TLinearFitter fitY(2, "pol1");
1284   static TLinearFitter fitZ(2, "pol1");
1285   fitY.StoreData(kFALSE);
1286   fitZ.StoreData(kFALSE);
1287   //
1288   //
1289    Int_t firstRow = 0, lastRow = 0;
1290    Int_t minRow = 100;
1291    Float_t xcenter = 0;
1292    const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
1293    const Float_t kedgey = 4.;
1294    if (padType == 0) {
1295       firstRow = 0;
1296       lastRow = fgTPCparam->GetNRowLow();
1297       xcenter = 108.47;
1298    }
1299    if (padType == 1) {
1300       firstRow = fgTPCparam->GetNRowLow();
1301       lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1302       xcenter = 166.60;
1303    }
1304    if (padType == 2) {
1305       firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1306       lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
1307       xcenter = 222.6;
1308    }
1309    minRow = (lastRow - firstRow) / 2;
1310    //
1311    //
1312    Int_t nclusters = 0;
1313    Int_t nclustersNE = 0; // number of not edge clusters
1314    Int_t lastSector = -1;
1315    Float_t amplitudeQ[100];
1316    Float_t amplitudeM[100];
1317    Int_t rowIn[100];
1318    Int_t index[100];
1319    //
1320    //
1321    fitY.ClearPoints();
1322    fitZ.ClearPoints();
1323
1324    for (Int_t iCluster = firstRow; iCluster < lastRow; iCluster++) {
1325       AliTPCclusterMI* cluster = track->GetClusterPointer(iCluster);
1326       if (cluster) {
1327          Int_t detector = cluster->GetDetector() ;
1328          if (lastSector == -1) lastSector = detector;
1329          if (lastSector != detector) continue;
1330          amplitudeQ[nclusters] = GetQNorm(cluster);
1331          amplitudeM[nclusters] = GetMaxNorm(cluster);
1332          rowIn[nclusters] = iCluster;
1333          nclusters++;
1334          Double_t dx = cluster->GetX() - xcenter;
1335          Double_t y = cluster->GetY();
1336          Double_t z = cluster->GetZ();
1337          fitY.AddPoint(&dx, y);
1338          fitZ.AddPoint(&dx, z);
1339          meanPos[0] += dx;
1340          meanPos[1] += dx;
1341          meanPos[2] += y;
1342          meanPos[3] += y*y;
1343          meanPos[4] += z;
1344          meanPos[5] += z*z;
1345          if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) nclustersNE++;
1346       }
1347    }
1348
1349    if (nclusters < minRow / 2) return kFALSE;
1350    if (nclustersNE < minRow / 2) return kFALSE;
1351    for (Int_t i = 0; i < 6; i++) meanPos[i] /= Double_t(nclusters);
1352    fitY.Eval();
1353    fitZ.Eval();
1354    fitY.GetParameters(parY);
1355    fitZ.GetParameters(parZ);
1356    //
1357    // calculate truncated mean
1358    //
1359    TMath::Sort(nclusters, amplitudeQ, index, kFALSE);
1360    //
1361    //
1362    //
1363    Float_t ndedx[5];
1364    for (Int_t i = 0; i < 5; i++) {
1365       dedxQ[i] = 0;
1366       dedxM[i] = 0;
1367       ndedx[i] = 0;
1368    }
1369    //
1370    // dedx calculation
1371    //
1372    Int_t inonEdge = 0;
1373    for (Int_t i = 0; i < nclusters; i++) {
1374       Int_t rowSorted = rowIn[index[i]];
1375       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1376
1377       if (TMath::Abs(cluster->GetY()) > cluster->GetX()*ktany - kedgey) continue;  //don't take edge clusters
1378       inonEdge++;
1379       if (inonEdge < nclustersNE * 0.5) {
1380          ndedx[0]++;
1381          dedxQ[0] += amplitudeQ[index[i]];
1382          dedxM[0] += amplitudeM[index[i]];
1383       }
1384       if (inonEdge < nclustersNE * 0.6) {
1385          ndedx[1]++;
1386          dedxQ[1] += amplitudeQ[index[i]];
1387          dedxM[1] += amplitudeM[index[i]];
1388       }
1389       if (inonEdge < nclustersNE * 0.7) {
1390          ndedx[2]++;
1391          dedxQ[2] += amplitudeQ[index[i]];
1392          dedxM[2] += amplitudeM[index[i]];
1393       }
1394       if (inonEdge < nclustersNE * 0.8) {
1395          ndedx[3]++;
1396          dedxQ[3] += amplitudeQ[index[i]];
1397          dedxM[3] += amplitudeM[index[i]];
1398       }
1399       if (inonEdge < nclustersNE * 0.9) {
1400          ndedx[4]++;
1401          dedxQ[4] += amplitudeQ[index[i]];
1402          dedxM[4] += amplitudeM[index[i]];
1403       }
1404    }
1405    for (Int_t i = 0; i < 5; i++) {
1406       dedxQ[i] /= ndedx[i];
1407       dedxM[i] /= ndedx[i];
1408    }
1409    TTreeSRedirector * cstream =  GetDebugStreamer();
1410    inonEdge = 0;
1411    Float_t momenta = track->GetP();
1412    Float_t mdedx = track->GetdEdx();
1413    for (Int_t i = 0; i < nclusters; i++) {
1414       Int_t rowSorted = rowIn[index[i]];
1415       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1416       if (!cluster) {
1417          printf("Problem\n");
1418          continue;
1419       }
1420       if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) inonEdge++;
1421       Float_t dedge = cluster->GetX()*ktany - TMath::Abs(cluster->GetY());
1422       Float_t fraction = Float_t(i) / Float_t(nclusters);
1423       Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
1424
1425       AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
1426       Float_t gain = GetGain(cluster);
1427       if (cstream) (*cstream) << "dEdx" <<
1428         "run="<<fRun<<              //  run number
1429         "event="<<fEvent<<          //  event number
1430         "time="<<fTime<<            //  time stamp of event
1431         "trigger="<<fTrigger<<      //  trigger
1432         "mag="<<fMagF<<             //  magnetic field
1433
1434                      "Cl.=" << cluster <<           // cluster of interest
1435                      "gain="<<gain<<                // gain at cluster position
1436                      "P=" << momenta <<             // track momenta
1437                      "dedx=" << mdedx <<            // mean dedx - corrected for angle
1438                      "IPad=" << padType <<          // pad type 0..2
1439                      "xc=" << xcenter <<            // x center of chamber
1440                      "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
1441                      "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
1442                      "fraction=" << fraction <<     // fraction - order in statistic (0,1)
1443                      "fraction2=" << fraction2 <<   // fraction - order in statistic (0,1)
1444                      "dedge=" << dedge <<           // distance to the edge
1445                      "parY.=" << &parY <<           // line fit
1446                      "parZ.=" << &parZ <<           // line fit
1447                      "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
1448                      "\n";
1449    }
1450
1451    if (cstream) (*cstream) << "dEdxT" <<
1452      "run="<<fRun<<              //  run number
1453      "event="<<fEvent<<          //  event number
1454      "time="<<fTime<<            //  time stamp of event
1455      "trigger="<<fTrigger<<      //  trigger
1456      "mag="<<fMagF<<             //  magnetic field
1457      "P=" << momenta <<             // track momenta
1458      "npoints="<<inonEdge<<         // number of points
1459      "sector="<<lastSector<<        // sector number
1460      "dedx=" << mdedx <<            // mean dedx - corrected for angle
1461      "IPad=" << padType <<          // pad type 0..2
1462      "xc=" << xcenter <<            // x center of chamber
1463      "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
1464      "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
1465      "parY.=" << &parY <<           // line fit
1466      "parZ.=" << &parZ <<           // line fit
1467      "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
1468      "\n";
1469
1470    sector = lastSector;
1471    npoints = inonEdge;
1472    return kTRUE;
1473 }
1474
1475 void AliTPCcalibTracksGain::AddTracklet(UInt_t sector, UInt_t padType,TVectorD &dedxQ, TVectorD &dedxM,TVectorD& parY, TVectorD& parZ, TVectorD& meanPos){
1476   //
1477   // Add measured point - dedx to the fitter
1478   //
1479   //
1480   //chain->SetAlias("dr","(250-abs(meanPos.fElements[4]))/250");
1481   //chain->SetAlias("tz","(0+abs(parZ.fElements[1]))");
1482   //chain->SetAlias("ty","(0+abs(parY.fElements[1]))");
1483   //chain->SetAlias("corrg","sqrt((1+ty^2)*(1+tz^2))");
1484   //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);
1485
1486   Double_t xxx[100];
1487   //
1488   // z and angular part
1489   //
1490
1491   xxx[0] = (250.-TMath::Abs(meanPos[4]))/250.;
1492   xxx[1] = TMath::Abs(parY[1]);
1493   xxx[2] = TMath::Abs(parZ[1]);
1494   xxx[3] = xxx[0]*xxx[1];
1495   xxx[4] = xxx[0]*xxx[2];
1496   xxx[5] = xxx[1]*xxx[2];
1497   xxx[6] = xxx[0]*xxx[0];
1498   xxx[7] = xxx[1]*xxx[1];
1499   xxx[8] = xxx[2]*xxx[2];
1500   //
1501   // chamber part
1502   //
1503   Int_t tsector = sector%36;
1504   for (Int_t i=0;i<35;i++){
1505     xxx[9+i]=(i==tsector)?1:0;
1506   }
1507   TLinearFitter *fitterM = fFitter0M;
1508   if (padType==1) fitterM=fFitter1M;
1509   if (padType==2) fitterM=fFitter2M;
1510   fitterM->AddPoint(xxx,dedxM[1]);
1511   //
1512   TLinearFitter *fitterT = fFitter0T;
1513   if (padType==1) fitterT = fFitter1T;
1514   if (padType==2) fitterT = fFitter2T;
1515   fitterT->AddPoint(xxx,dedxQ[1]);
1516   //
1517   TLinearFitter *dfitterM = fDFitter0M;
1518   if (padType==1) dfitterM=fDFitter1M;
1519   if (padType==2) dfitterM=fDFitter2M;
1520   dfitterM->AddPoint(xxx,dedxM[1]);
1521   //
1522   TLinearFitter *dfitterT = fDFitter0T;
1523   if (padType==1) dfitterT = fDFitter1T;
1524   if (padType==2) dfitterT = fDFitter2T;
1525   dfitterT->AddPoint(xxx,dedxQ[1]);
1526 }
1527
1528
1529 TGraph *AliTPCcalibTracksGain::CreateAmpGraph(Int_t ipad, Bool_t qmax){
1530   //
1531   // create the amplitude graph
1532   // The normalized amplitudes are extrapolated to the 0 angle (y,z)  and 0 drift length
1533   //
1534
1535   TVectorD vec;
1536   if (qmax){
1537     if (ipad==0) fFitter0M->GetParameters(vec);
1538     if (ipad==1) fFitter1M->GetParameters(vec);
1539     if (ipad==2) fFitter2M->GetParameters(vec);
1540   }else{
1541     if (ipad==0) fFitter0T->GetParameters(vec);
1542     if (ipad==1) fFitter1T->GetParameters(vec);
1543     if (ipad==2) fFitter2T->GetParameters(vec);
1544   }
1545
1546   Float_t amp[36];
1547   Float_t sec[36];
1548   for (Int_t i=0;i<35;i++){
1549     sec[i]=i;
1550     amp[i]=vec[10+i]+vec[0];
1551   }
1552   amp[35]=vec[0];
1553   Float_t mean = TMath::Mean(36,amp);
1554   for (Int_t i=0;i<36;i++){
1555     sec[i]=i;
1556     amp[i]=(amp[i]-mean)/mean;
1557   }
1558   TGraph *gr = new TGraph(36,sec,amp);
1559   return gr;
1560 }
1561
1562
1563 void   AliTPCcalibTracksGain::UpdateClusterParam(AliTPCClusterParam* clparam){
1564   //
1565   //   SetQ normalization parameters
1566   //
1567   //  void SetQnorm(Int_t ipad, Int_t itype,  TVectorD * norm);
1568
1569   TVectorD vec;
1570
1571   //
1572   fDFitter0T->Eval();
1573   fDFitter1T->Eval();
1574   fDFitter2T->Eval();
1575   fDFitter0M->Eval();
1576   fDFitter1M->Eval();
1577   fDFitter2M->Eval();
1578   fDFitter0T->GetParameters(vec);
1579   clparam->SetQnorm(0,0,&vec);
1580   fDFitter1T->GetParameters(vec);
1581   clparam->SetQnorm(1,0,&vec);
1582   fDFitter2T->GetParameters(vec);
1583   clparam->SetQnorm(2,0,&vec);
1584   //
1585   fDFitter0M->GetParameters(vec);
1586   clparam->SetQnorm(0,1,&vec);
1587   fDFitter1M->GetParameters(vec);
1588   clparam->SetQnorm(1,1,&vec);
1589   fDFitter2M->GetParameters(vec);
1590   clparam->SetQnorm(2,1,&vec);
1591   //
1592
1593 }
1594
1595
1596 void   AliTPCcalibTracksGain::Analyze(){
1597
1598  Evaluate();
1599
1600 }
1601
1602
1603
1604 TVectorD *  AliTPCcalibTracksGain::MakeQPosNorm(TTree * chain0, Int_t ipad, Bool_t isMax, Int_t maxPoints, Int_t verbose){
1605   //
1606   // Input parameters
1607   // chain0 - the tree with information  -Debug stream
1608   // ipad   - 0 IROC
1609   //        - 1 OROC medium
1610   //        - 2 OROC LONG
1611   // isMax  - kFALSE - total charge param
1612   //          kTRUE  - Max charge param
1613   //
1614   // maxPoints - number of points for fit
1615   //
1616   // verbose -
1617   //
1618   /* e.g
1619     ipad=0
1620     isMax=kTRUE;
1621     maxPoints=1000000;
1622   */
1623   // Make Q normalization as function of following parameters
1624   // 1 - dp   - relative pad position
1625   // 2 - dt   - relative time position
1626   // 3 - di   - drift length (norm to 1);
1627   // 4 - dq0  - Tot/Max charge
1628   // 5 - dq1  - Max/Tot charge
1629   // 6 - sy   - sigma y - shape
1630   // 7 - sz   - sigma z - shape
1631   //
1632   // Coeficient of Taylor expansion fitted
1633   // Fit parameters returned as TVectorD
1634   // Fit parameters to be used in corresponding correction function
1635   // in AliTPCclusterParam
1636   //
1637   //
1638   TStatToolkit toolkit;
1639   Double_t chi2;
1640   TVectorD fitParam;
1641   TMatrixD covMatrix;
1642   Int_t npoints;
1643   TCut cutA("dedge>3&&fraction2<0.7");
1644   chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1645   chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1646   chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1647   chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1648   chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1649   chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1650   chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1651   //
1652   TString fstring="";
1653   fstring+="dp++";                               //1
1654   fstring+="dt++";                               //2
1655   fstring+="dp*dp++";                             //3
1656   fstring+="dt*dt++";                             //4
1657   fstring+="dt*dt*dt++";                             //5
1658   fstring+="dp*dt++";                            //6
1659   fstring+="dp*dt*dt++";                          //7
1660   fstring+="(dq0)++";                            //8
1661   fstring+="(dq1)++";                            //9
1662   //
1663   //
1664   fstring+="dp*dp*(di)++";                        //10
1665   fstring+="dt*dt*(di)++";                        //11
1666   fstring+="dp*dp*sy++";                          //12
1667   fstring+="dt*sz++";                          //13
1668   fstring+="dt*dt*sz++";                          //14
1669   fstring+="dt*dt*dt*sz++";                          //15
1670   //
1671   fstring+="dp*dp*1*sy*sz++";                     //16
1672   fstring+="dt*sy*sz++";                       //17
1673   fstring+="dt*dt*sy*sz++";                       //18
1674   fstring+="dt*dt*dt*sy*sz++";                       //19
1675   //
1676   fstring+="dp*dp*(dq0)++";                       //20
1677   fstring+="dt*1*(dq0)++";                       //21
1678   fstring+="dt*dt*(dq0)++";                       //22
1679   fstring+="dt*dt*dt*(dq0)++";                       //23
1680   //
1681   fstring+="dp*dp*(dq1)++";                       //24
1682   fstring+="dt*(dq1)++";                       //25
1683   fstring+="dt*dt*(dq1)++";                       //26
1684   fstring+="dt*dt*dt*(dq1)++";                       //27
1685
1686   TString var;
1687   if (isMax)  var = "Cl.fMax/gain/dedxM.fElements[2]";
1688   if (!isMax) var = "Cl.fQ/gain/dedxQ.fElements[2]";
1689   TString cutP="IPad==";
1690   cutP+=ipad;
1691   //
1692   TString *strq0 = toolkit.FitPlane(chain0,var.Data(),fstring.Data(), cutP.Data()+cutA, chi2,npoints,fitParam,covMatrix,-1,0,maxPoints);
1693   //
1694   //
1695   if (verbose){
1696     printf("Chi2/npoints = %f",TMath::Sqrt(chi2/npoints));
1697     printf("\nFit function\n:%s\n",strq0->Data());
1698   }
1699   TVectorD *vec = new TVectorD(fitParam);
1700   return vec;
1701 }
1702
1703 void  AliTPCcalibTracksGain::MakeQPosNormAll(TTree * chain, AliTPCClusterParam * param, Int_t maxPoints, Int_t verbose){
1704   //
1705   // Fill the content of the of the AliTPCclusterParam
1706   // with fitted values of corrections
1707   //
1708   param->fPosQTnorm[0] = MakeQPosNorm(chain,0,kTRUE,100000,kTRUE);
1709   param->fPosQTnorm[1] = MakeQPosNorm(chain,1,kTRUE,100000,kTRUE);
1710   param->fPosQTnorm[2] = MakeQPosNorm(chain,1,kTRUE,100000,kTRUE);
1711   //
1712   param->fPosQMnorm[0] = MakeQPosNorm(chain,0,kFALSE,100000,kTRUE);
1713   param->fPosQMnorm[1] = MakeQPosNorm(chain,1,kFALSE,100000,kTRUE);
1714   param->fPosQMnorm[2] = MakeQPosNorm(chain,2,kFALSE,100000,kTRUE);
1715 }
1716
1717
1718
1719 /*
1720
1721  Position correction fit:
1722  //
1723 TStatToolkit toolkit;
1724 Double_t chi2;
1725 TVectorD fitParam;
1726 TMatrixD covMatrix;
1727 Int_t npoints;
1728 //
1729 TCut cutA("dedge>3&&fraction2<0.7");
1730 chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1731 chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1732 chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1733 chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1734 chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1735 chain0->SetAlias("sy","(0.2/sqrt(0.01^2+Cl.fSigmaY2))");
1736 chain0->SetAlias("sz","(0.2/sqrt(0.01^2+Cl.fSigmaZ2))");
1737
1738 TString fstring="";
1739
1740 fstring+="dp++";                               //1
1741 fstring+="dt++";                               //2
1742 fstring+="dp*dp++";                             //3
1743 fstring+="dt*dt++";                             //4
1744 fstring+="dt*dt*dt++";                             //5
1745 fstring+="dp*dt++";                            //6
1746 fstring+="dp*dt*dt++";                          //7
1747 fstring+="(dq0)++";                            //8
1748 fstring+="(dq1)++";                            //9
1749 //
1750 //
1751 fstring+="dp*dp*(di)++";                        //10
1752 fstring+="dt*dt*(di)++";                        //11
1753 fstring+="dp*dp*sy++";                          //12
1754 fstring+="dt*sz++";                          //13
1755 fstring+="dt*dt*sz++";                          //14
1756 fstring+="dt*dt*dt*sz++";                          //15
1757 //
1758 fstring+="dp*dp*1*sy*sz++";                     //16
1759 fstring+="dt*sy*sz++";                       //17
1760 fstring+="dt*dt*sy*sz++";                       //18
1761 fstring+="dt*dt*dt*sy*sz++";                       //19
1762 //
1763 fstring+="dp*dp*(dq0)++";                       //20
1764 fstring+="dt*1*(dq0)++";                       //21
1765 fstring+="dt*dt*(dq0)++";                       //22
1766 fstring+="dt*dt*dt*(dq0)++";                       //23
1767
1768 fstring+="dp*dp*(dq1)++";                       //24
1769 fstring+="dt*(dq1)++";                       //25
1770 fstring+="dt*dt*(dq1)++";                       //26
1771 fstring+="dt*dt*dt*(dq1)++";                       //27
1772
1773
1774  TString *strq0 = toolkit.FitPlane(chain0,"Cl.fMax/gain/dedxM.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
1775  TString *strqt0 = toolkit.FitPlane(chain0,"Cl.fQ/gain/dedxQ.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
1776
1777  chain0->SetAlias("qcorM0",strq0->Data());
1778  chain0->SetAlias("qcorT0",strqt0->Data());
1779 //chain0->SetAlias("mmqcorM0","min(max(qcorM0,0.75),1.15)");
1780  chain0->Draw("(Cl.fMax/gain/dedxM.fElements[2]):min(max(qcorM0,0.75),1.15)","IPad==0"+cutA,"prof",100000)
1781
1782  fraction05 -
1783      sigma                  0.2419
1784      sigma fit              0.2302
1785      sigma fit with shape   0.2257
1786  fraction 07
1787      qtot sigma                  0.322
1788      qmax sigma                  0.292
1789      qmax sigma fit              0.2702
1790      qmax sigma fit+ratio        0.2638
1791
1792 */
1793