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   AliTPCCalROC * roc = fGainMap->GetCalROC(cl->GetDetector());
 595   Int_t irow = cl->GetRow();
 596   if (roc){
 597     if (irow < 63) { // IROC
 598       factor = roc->GetValue(irow, TMath::Nint(cl->GetPad()));
 599     } else {         // OROC
 600       factor = roc->GetValue(irow - 63, TMath::Nint(cl->GetPad()));
 601     }
 602   }
 603   if (factor<0.1) factor=0.1;
 604   return factor;
 605 }
 606
 607
 608 Float_t   AliTPCcalibTracksGain::GetMaxNorm(AliTPCclusterMI * cl){
 609   //
 610   // Get normalized amplituded if the gain map provided
 611   //
 612   return cl->GetMax()/GetGain(cl);
 613 }
 614
 615
 616 Float_t   AliTPCcalibTracksGain::GetQNorm(AliTPCclusterMI * cl){
 617   //
 618   // Get normalized amplituded if the gain map provided
 619   //
 620   return cl->GetQ()/GetGain(cl);
 621 }
 622
 623
 624
 625 void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
 626    //
 627    // Adds another AliTPCcalibTracksGain object to this object.
 628    //
 629
 630    fSimpleFitter->Add(cal->fSimpleFitter);
 631    fSqrtFitter->Add(cal->fSqrtFitter);
 632    fLogFitter->Add(cal->fLogFitter);
 633    fSingleSectorFitter->Add(cal->fSingleSectorFitter);
 634    //
 635    //
 636    //
 637    fFitter0M->Add(cal->fFitter0M);
 638    fFitter1M->Add(cal->fFitter1M);
 639    fFitter2M->Add(cal->fFitter2M);
 640    fFitter0T->Add(cal->fFitter0T);
 641    fFitter1T->Add(cal->fFitter1T);
 642    fFitter2T->Add(cal->fFitter2T);
 643    //
 644    fDFitter0M->Add(cal->fDFitter0M);
 645    fDFitter1M->Add(cal->fDFitter1M);
 646    fDFitter2M->Add(cal->fDFitter2M);
 647    fDFitter0T->Add(cal->fDFitter0T);
 648    fDFitter1T->Add(cal->fDFitter1T);
 649    fDFitter2T->Add(cal->fDFitter2T);
 650    //
 651    //
 652    // histograms
 653    //
 654    for (Int_t i=0; i<73; i++){
 655      TH1F *his,*hism;
 656      his = (TH1F*)fArrayQM->At(i);
 657      hism =  (TH1F*)cal->fArrayQM->At(i);
 658      if (his && hism) his->Add(hism);
 659      his =  (TH1F*)fArrayQT->At(i);
 660      hism =  (TH1F*)cal->fArrayQT->At(i);
 661      if (his && hism) his->Add(hism);
 662    }
 663    //
 664    //
 665    for (Int_t i=0; i<37; i++){
 666      TProfile *his,*hism;
 667      his = (TProfile*)fProfileArrayQM->At(i);
 668      hism = (TProfile*)cal->fProfileArrayQM->At(i);
 669      if (his && hism) his->Add(hism);
 670      his = (TProfile*)fProfileArrayQT->At(i);
 671      hism = (TProfile*)cal->fProfileArrayQT->At(i);
 672      if (his && hism) his->Add(hism);
 673    }
 674    //
 675    //
 676    for (Int_t i=0; i<37; i++){
 677      TProfile2D *his,*hism;
 678      his = (TProfile2D*)fProfileArrayQM2D->At(i);
 679      hism = (TProfile2D*)cal->fProfileArrayQM2D->At(i);
 680      if (his && hism) his->Add(hism);
 681      his = (TProfile2D*)fProfileArrayQT2D->At(i);
 682      hism = (TProfile2D*)cal->fProfileArrayQT2D->At(i);
 683      if (his && hism) his->Add(hism);
 684    }
 685    //
 686    // this will be gone for the a new ROOT version > v5-17-05
 687    for (UInt_t iSegment = 0; iSegment < 36; iSegment++) {
 688       fNShortClusters[iSegment] += cal->fNShortClusters[iSegment];
 689       fNMediumClusters[iSegment] += cal->fNMediumClusters[iSegment];
 690       fNLongClusters[iSegment] += cal->fNLongClusters[iSegment];
 691    }
 692
 693    // just for debugging, remove me
 694    fTotalTracks += cal->fTotalTracks;
 695    fAcceptedTracks += cal->fAcceptedTracks;
 696    fDebugCalPadRaw->Add(cal->fDebugCalPadRaw);
 697    fDebugCalPadCorr->Add(cal->fDebugCalPadCorr);
 698
 699 }
 700
 701 void AliTPCcalibTracksGain::AddTrack(AliTPCseed* seed) {
 702    //
 703    // The clusters making up the track (seed) are added to various fit functions.
 704    // See AddCluster(...) for more detail.
 705    //
 706
 707    DumpTrack(seed);
 708    //
 709    // simple histograming part
 710    for (Int_t i=0; i<159; i++){
 711      AliTPCclusterMI* cluster = seed->GetClusterPointer(i);
 712      if (cluster) AddCluster(cluster);
 713    }
 714 }
 715
 716 void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster){
 717   //
 718   // Adding cluster information to the simple histograms
 719   // No correction, fittings are applied
 720
 721
 722   Float_t kThreshold=5;
 723   if (cluster->GetX()<=0) return;
 724   if (cluster->GetQ()<=kThreshold) return;
 725   //
 726
 727   Int_t sector = cluster->GetDetector();
 728   TH1F  * his;
 729   his = GetQT(sector);
 730   if (his) his->Fill(GetQNorm(cluster));
 731   his = GetQT(-1);
 732   if (his) his->Fill(GetQNorm(cluster));
 733   his = GetQM(sector);
 734   if (his) his->Fill(GetMaxNorm(cluster));
 735   his = GetQM(-1);
 736   if (his) his->Fill(GetMaxNorm(cluster));
 737   //
 738   sector = sector%36;
 739   TProfile *prof;
 740   prof = GetProfileQT(sector);
 741   if (prof) prof->Fill(cluster->GetX(),GetQNorm(cluster));
 742   prof = GetProfileQT(-1);
 743   if (prof) prof->Fill(cluster->GetX(),GetQNorm(cluster));
 744   prof = GetProfileQM(sector);
 745   if (prof) prof->Fill(cluster->GetX(),GetMaxNorm(cluster));
 746   prof = GetProfileQM(-1);
 747   if (prof) prof->Fill(cluster->GetX(),GetMaxNorm(cluster));
 748   //
 749   Float_t phi = cluster->GetY()/cluster->GetX();
 750   TProfile2D *prof2;
 751   prof2 = GetProfileQT2D(sector);
 752   if (prof2) prof2->Fill(cluster->GetX(),phi,GetQNorm(cluster));
 753   prof2 = GetProfileQT2D(-1);
 754   if (prof2) prof2->Fill(cluster->GetX(),phi,GetQNorm(cluster));
 755   prof2 = GetProfileQM2D(sector);
 756   if (prof2) prof2->Fill(cluster->GetX(),phi,GetMaxNorm(cluster));
 757   prof2 = GetProfileQM2D(-1);
 758   if (prof2) prof2->Fill(cluster->GetX(),phi,GetMaxNorm(cluster));
 759
 760   //
 761 }
 762
 763
 764
 765 void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t /*momenta*/, Float_t/* mdedx*/, Int_t padType,
 766                                        Float_t xcenter, TVectorD& dedxQ, TVectorD& /*dedxM*/, Float_t /*fraction*/, Float_t fraction2, Float_t dedge,
 767                                        TVectorD& /*parY*/, TVectorD& /*parZ*/, TVectorD& meanPos) {
 768    //
 769    // Adds cluster to the appropriate fitter for later analysis.
 770    // The charge used for the fit is the maximum charge for this specific cluster or the
 771    // accumulated charge per cluster, depending on the value of fgkUseTotalCharge.
 772    // Depending on the pad size where the cluster is registered, the value will be put in
 773    // the appropriate fitter. Furthermore, for each pad size three different types of fitters
 774    // are used. The fit functions are the same for all fitters (parabolic functions), but the value
 775    // added to each fitter is different. The simple fitter gets the charge plugged in as is, the sqrt fitter
 776    // gets the square root of the charge, and the log fitter gets fgkM*(1+q/fgkM), where q is the original charge
 777    // and fgkM==25.
 778    //
 779   Float_t kedge     = 3;
 780   Float_t kfraction = 0.7;
 781   Int_t   kinorm    = 2;
 782
 783
 784   // Where to put selection on threshold?
 785   // Defined by the Q/dEdxT variable - see debug streamer:
 786   //
 787   // Debug stream variables:  (Where tu cut ?)
 788   // chain0->Draw("Cl.fQ/dedxQ.fElements[1]>>his(100,0,3)","fraction2<0.6&&dedge>3","",1000000);
 789   //         mean 1  sigma 0.25
 790   // chain0->Draw("Cl.fMax/dedxM.fElements[1]>>his(100,0,3)","fraction2<0.6&&dedge>3","",1000000)
 791   //         mean 1  sigma 0.25
 792   // chain0->Draw("Cl.fQ/dedxQ.fElements[2]>>his(100,0,3)","fraction2<0.7&&dedge>3","",1000000)
 793   //         mean 1 sigma 0.29
 794   // chain0->Draw("Cl.fMax/dedxM.fElements[2]>>his(100,0,3)","fraction2<0.7&&dedge>3","",1000000)
 795   //         mean 1 sigma 0.27
 796   // chain0->Draw("Cl.fQ/dedxQ.fElements[3]>>his(100,0,3)","fraction2<0.8&&dedge>3","",1000000)
 797   //         mean 1 sigma 0.32
 798   //
 799   // chain0->Draw("Cl.fQ/dedxQ.fElements[4]>>his(100,0,3)","fraction2<0.9&&dedge>3","",1000000)
 800   //         mean 1 sigma 0.4
 801
 802   // Fraction choosen 0.7
 803
 804    if (!cluster) {
 805       Error("AddCluster", "Cluster not valid.");
 806       return;
 807    }
 808
 809    if (dedge < kedge) return;
 810    if (fraction2 > kfraction) return;
 811
 812    //Int_t padType = GetPadType(cluster->GetX());
 813    Double_t xx[7];
 814    //Double_t centerPad[2] = {0};
 815    //AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
 816    //xx[0] = cluster->GetX() - centerPad[0];
 817    //xx[1] = cluster->GetY() - centerPad[1];
 818    xx[0] = cluster->GetX() - xcenter;
 819    xx[1] = cluster->GetY();
 820    xx[2] = xx[0] * xx[0];
 821    xx[3] = xx[1] * xx[1];
 822    xx[4] = xx[0] * xx[1];
 823    xx[5] = TMath::Abs(cluster->GetZ()) - TMath::Abs(meanPos[4]);
 824    xx[6] = xx[5] * xx[5];
 825    //
 826    // Update profile histograms
 827    //
 828
 829    //
 830    // Update fitters
 831    //
 832    Int_t segment = cluster->GetDetector() % 36;
 833    Double_t q = fgkUseTotalCharge ? ((Double_t)(GetQNorm(cluster))) : ((Double_t)(GetMaxNorm(cluster)));  // note: no normalization to pad size!
 834
 835    // just for debugging
 836    Int_t row = 0;
 837    Int_t pad = 0;
 838    GetRowPad(cluster->GetX(), cluster->GetY(), row, pad);
 839    fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
 840
 841    // correct charge by normalising to mean charge per track
 842    q /= dedxQ[kinorm];
 843
 844    // just for debugging
 845    fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
 846
 847    Double_t sqrtQ = TMath::Sqrt(q);
 848    Double_t logQ = fgkM * TMath::Log(1 + q / fgkM);
 849    TLinearFitter * fitter =0;
 850    //
 851    fitter = fSimpleFitter->GetFitter(segment, padType);
 852    fitter->AddPoint(xx, q);
 853    //
 854    fitter = fSqrtFitter->GetFitter(segment, padType);
 855    fitter->AddPoint(xx, sqrtQ);
 856    //
 857    fitter = fLogFitter->GetFitter(segment, padType);
 858    fitter->AddPoint(xx, logQ);
 859    //
 860    fitter=fSingleSectorFitter->GetFitter(0, padType);
 861    fitter->AddPoint(xx, q);
 862
 863    // this will be gone for the a new ROOT version > v5-17-05
 864    if (padType == kShortPads)
 865       fNShortClusters[segment]++;
 866    if (padType == kMediumPads)
 867       fNMediumClusters[segment]++;
 868    if (padType == kLongPads)
 869       fNLongClusters[segment]++;
 870 }
 871
 872 void AliTPCcalibTracksGain::Evaluate(Bool_t robust, Double_t frac) {
 873    //
 874    // Evaluates all fitters contained in this object.
 875    // If the robust option is set to kTRUE a robust fit is performed with frac as
 876    // the minimal fraction of good points (see TLinearFitter::EvalRobust for details).
 877    // Beware: Robust fitting is much slower!
 878    //
 879
 880    fSimpleFitter->Evaluate(robust, frac);
 881    fSqrtFitter->Evaluate(robust, frac);
 882    fLogFitter->Evaluate(robust, frac);
 883    fSingleSectorFitter->Evaluate(robust, frac);
 884    fFitter0M->Eval();
 885    fFitter1M->Eval();
 886    fFitter2M->Eval();
 887    fFitter0T->Eval();
 888    fFitter1T->Eval();
 889    fFitter2T->Eval();
 890    //
 891    fDFitter0M->Eval();
 892    fDFitter1M->Eval();
 893    fDFitter2M->Eval();
 894    fDFitter0T->Eval();
 895    fDFitter1T->Eval();
 896    fDFitter2T->Eval();
 897 }
 898
 899 AliTPCCalPad* AliTPCcalibTracksGain::CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 900   //
 901   // Creates the calibration object AliTPCcalPad using fitted parameterization
 902   //
 903    TObjArray tpc(72);
 904    for (UInt_t iSector = 0; iSector < 72; iSector++)
 905       tpc.Add(CreateFitCalROC(iSector, fitType, undoTransformation, normalizeToPadSize));
 906    return new AliTPCCalPad(&tpc);
 907 }
 908
 909 AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 910   //
 911   // Create the AliTPCCalROC with the values per pad
 912   // sector  - sector of interest
 913   // fitType -
 914   //
 915
 916    TVectorD par(8);
 917    if (sector < 36) {
 918       GetParameters(sector % 36, 0, fitType, par);
 919       return CreateFitCalROC(sector, 0, par, fitType, undoTransformation, normalizeToPadSize);
 920    }
 921    else {
 922       GetParameters(sector % 36, 1, fitType, par);
 923       AliTPCCalROC* roc1 = CreateFitCalROC(sector, 1, par, fitType, undoTransformation, normalizeToPadSize);
 924       GetParameters(sector % 36, 2, fitType, par);
 925       AliTPCCalROC* roc2 = CreateFitCalROC(sector, 2, par, fitType, undoTransformation, normalizeToPadSize);
 926       AliTPCCalROC* roc3 = CreateCombinedCalROC(roc1, roc2);
 927       delete roc1;
 928       delete roc2;
 929       return roc3;
 930    }
 931 }
 932
 933 AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t padType, TVectorD &fitParam, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 934    //
 935    // This function is essentially a copy of AliTPCCalROC::CreateGlobalFitCalROC(...), with the
 936    // modifications, that the center of the region of same pad size is used as the origin
 937    // of the fit function instead of the center of the ROC.
 938    // The possibility of a linear fit is removed as well because it is not needed.
 939    // Only values for pads with the given pad size are calculated, the rest is 0.
 940    // Set undoTransformation for undoing the transformation that was applied to the
 941    // charge values before they were put into the fitter (thus allowing comparison to the original
 942    // charge values). For fitType use 0 for the simple fitter, 1 for the sqrt fitter, 2 for the log fitter.
 943    // If normalizeToPadSize is true, the values are normalized to the pad size.
 944    // Please be aware, that you even need to specify the fitType if you want to normalize to the pad size without
 945    // undoing the transformation (because normalizing involves undoing the trafo first, then normalizing, then
 946    // applying the trafo again).
 947    // Please note: The normalization to the pad size is a simple linear scaling with the pad length, which
 948    //              actually doesn't describe reality!
 949    //
 950
 951    Float_t dlx, dly;
 952    Double_t centerPad[2] = {0};
 953    Float_t localXY[3] = {0};
 954    AliTPCROC* tpcROC = AliTPCROC::Instance();
 955    if ((padType == 0 && sector >= tpcROC->GetNInnerSector()) || (padType > 0 && sector < tpcROC->GetNInnerSector()) || sector >= tpcROC->GetNSector())
 956       return 0;
 957    AliTPCCalROC* lROCfitted = new AliTPCCalROC(sector);
 958    //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
 959    UInt_t startRow = 0;
 960    UInt_t endRow = 0;
 961    switch (padType) {
 962       case kShortPads:
 963          startRow = 0;
 964          endRow = lROCfitted->GetNrows();
 965          break;
 966       case kMediumPads:
 967          startRow = 0;
 968          endRow = 64;
 969          break;
 970       case kLongPads:
 971          startRow = 64;
 972          endRow = lROCfitted->GetNrows();
 973          break;
 974    }
 975
 976    AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
 977    Double_t value = 0;
 978    for (UInt_t irow = startRow; irow < endRow; irow++) {
 979       for (UInt_t ipad = 0; ipad < lROCfitted->GetNPads(irow); ipad++) {
 980          tpcROC->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 981          dlx = localXY[0] - centerPad[0];
 982          dly = localXY[1] - centerPad[1];
 983          value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
 984
 985          // Let q' = value be the transformed value without any pad size corrections,
 986          // let T be the transformation and let l be the pad size
 987          //    1) don't undo transformation, don't normalize: return q'
 988          //    2) undo transformation,       don't normalize: return T^{-1} q'
 989          //    3) undo transformation,       normalize:       return (T^{-1} q') / l
 990          //    4) don't undo transformation, normalize:       return T((T^{-1} q') / l)
 991          if (!undoTransformation && !normalizeToPadSize) {/* value remains unchanged */}  // (1)
 992          else {                                                                           // (2), (3), (4)
 993             //calculate T^{-1}
 994             switch (fitType) {
 995                case  0: /* value remains unchanged */ break;
 996                case  1: value = value * value; break;
 997                case  2: value = (TMath::Exp(value / fgkM) - 1) * fgkM; break;
 998                default: Error("CreateFitCalROC", "Wrong fit type."); break;
 999             }
1000             if (normalizeToPadSize) value /= GetPadLength(localXY[0]);                    // (3)
1001          }
1002          if (!undoTransformation && normalizeToPadSize) {                                 // (4)
1003             // calculate T
1004             switch (fitType) {
1005                case  0: /* value remains unchanged */ break;
1006                case  1: value = TMath::Sqrt(value); break;
1007                case  2: value = fgkM * TMath::Log(1 + value / fgkM); break;
1008                default: Error("CreateFitCalROC", "Wrong fit type."); break;
1009             }
1010          }
1011          lROCfitted->SetValue(irow, ipad, value);
1012       }
1013    }
1014    return lROCfitted;
1015 }
1016
1017 AliTPCCalROC* AliTPCcalibTracksGain::CreateCombinedCalROC(const AliTPCCalROC* roc1, const AliTPCCalROC* roc2) {
1018    //
1019    // Combines the medium pad size values of roc1 with the long pad size values of roc2 into a new
1020    // AliTPCCalROC. Returns a null pointer if any one of the ROCs is an IROC; issues a warning message
1021    // if the sectors of roc1 and roc2 don't match, but still continue and use the sector of roc1 as the
1022    // sector of the new ROC.
1023    //
1024
1025    if (!roc1 || !roc2) return 0;
1026    if ((Int_t)(roc1->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
1027    if ((Int_t)(roc2->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
1028    if (roc1->GetSector() != roc2->GetSector()) Warning("CreateCombinedCalROC", "Sector number mismatch.");
1029    AliTPCCalROC* roc = new AliTPCCalROC(roc1->GetSector());
1030
1031    for (UInt_t iRow = 0; iRow < 64; iRow++) {
1032       for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
1033          roc->SetValue(iRow, iPad, roc1->GetValue(iRow, iPad));
1034    }
1035    for (UInt_t iRow = 64; iRow < roc->GetNrows(); iRow++) {
1036       for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
1037          roc->SetValue(iRow, iPad, roc2->GetValue(iRow, iPad));
1038    }
1039    return roc;
1040 }
1041
1042 Bool_t AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) {
1043    //
1044    // Puts the fit parameters for the specified segment (IROC & OROC), padType and fitType
1045    // into the fitParam TVectorD (which should contain 8 elements).
1046    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1047    // Note: The fitter has to be evaluated first!
1048    //
1049   TLinearFitter * fitter = GetFitter(segment, padType, fitType);
1050   if (fitter){
1051     fitter->Eval();
1052     fitter->GetParameters(fitParam);
1053     return kTRUE;
1054   }else{
1055     Error("AliTPCcalibTracksGain::GetParameters",
1056           Form("Fitter%d_%d_%d not availble", segment, padType, fitType));
1057     return kFALSE;
1058   }
1059   return kFALSE;
1060 }
1061
1062 void AliTPCcalibTracksGain::GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError) {
1063    //
1064    // Puts the fit parameter errors for the specified segment (IROC & OROC), padType and fitType
1065    // into the fitError TVectorD (which should contain 8 elements).
1066    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1067    // Note: The fitter has to be evaluated first!
1068    //
1069
1070    GetFitter(segment, padType, fitType)->GetErrors(fitError);
1071    fitError *= TMath::Sqrt(GetRedChi2(segment, padType, fitType));
1072 }
1073
1074 Double_t AliTPCcalibTracksGain::GetRedChi2(UInt_t segment, UInt_t padType, UInt_t fitType) {
1075    //
1076    // Returns the reduced chi^2 value for the specified segment, padType and fitType.
1077    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1078    // Note: The fitter has to be evaluated first!
1079    //
1080
1081    // this will be gone for the a new ROOT version > v5-17-05
1082    Int_t lNClusters = 0;
1083    switch (padType) {
1084       case kShortPads:
1085          lNClusters = fNShortClusters[segment];
1086          break;
1087       case kMediumPads:
1088          lNClusters = fNMediumClusters[segment];
1089          break;
1090       case kLongPads:
1091          lNClusters = fNLongClusters[segment];
1092          break;
1093    }
1094    return GetFitter(segment, padType, fitType)->GetChisquare()/(lNClusters - 8);
1095 }
1096
1097 void AliTPCcalibTracksGain::GetCovarianceMatrix(UInt_t segment, UInt_t padType, UInt_t fitType, TMatrixD& covMatrix) {
1098    //
1099    // Returns the covariance matrix for the specified segment, padType, fitType.
1100    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1101    //
1102
1103    GetFitter(segment, padType, fitType)->GetCovarianceMatrix(covMatrix);
1104 }
1105
1106 TLinearFitter* AliTPCcalibTracksGain::GetFitter(UInt_t segment, UInt_t padType, UInt_t fitType) {
1107    //
1108    // Returns the TLinearFitter object for the specified segment, padType, fitType.
1109    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
1110    //
1111
1112    switch (fitType) {
1113       case kSimpleFitter:
1114          return fSimpleFitter->GetFitter(segment, padType);
1115       case kSqrtFitter:
1116          return fSqrtFitter->GetFitter(segment, padType);
1117       case kLogFitter:
1118          return fLogFitter->GetFitter(segment, padType);
1119       case 3:
1120          return fSingleSectorFitter->GetFitter(0, padType);
1121    }
1122    return 0;
1123 }
1124
1125 Double_t AliTPCcalibTracksGain::GetPadLength(Double_t lx) {
1126    //
1127    // The function returns 0.75 for an IROC, 1. for an OROC at medium pad size position,
1128    // 1.5 for an OROC at long pad size position, -1 if out of bounds.
1129    //
1130
1131    Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
1132    Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
1133    Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
1134    Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
1135    Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
1136    Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
1137
1138    // if IROC
1139    if (lx >= irocLow && lx <= irocUp) return 0.75;
1140    // if OROC medium pads
1141    if (lx >= orocLow1 && lx <= orocUp1) return 1.;
1142    // if OROC long pads
1143    if (lx >= orocLow2 && lx <= orocUp2) return 1.5;
1144    // if out of bounds
1145    return -1;
1146 }
1147
1148 Int_t AliTPCcalibTracksGain::GetPadType(Double_t lx) {
1149    //
1150    // The function returns 0 for an IROC, 1 for an OROC at medium pad size position,
1151    // 2 for an OROC at long pad size position, -1 if out of bounds.
1152    //
1153
1154    if (GetPadLength(lx) == 0.75) return 0;
1155    else if (GetPadLength(lx) == 1.) return 1;
1156    else if (GetPadLength(lx) == 1.5) return 2;
1157    return -1;
1158 }
1159
1160 // ONLY FOR DEBUGGING PURPOSES - REMOVE ME WHEN NOT NEEDED ANYMORE
1161 Bool_t AliTPCcalibTracksGain::GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad) {
1162    //
1163    // Calculate the row and pad number when the local coordinates are given.
1164    // Returns kFALSE if the position is out of range, otherwise return kTRUE.
1165    // WARNING: This function is preliminary and probably isn't very accurate!!
1166    //
1167
1168    Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
1169    //Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
1170    Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
1171    //Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
1172    Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
1173    //Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
1174
1175    if (GetPadType(lx) == 0) {
1176       row = (Int_t)((lx - irocLow) / fgTPCparam->GetInnerPadPitchLength());
1177       pad = (Int_t)((ly + fgTPCparam->GetYInner(row)) / fgTPCparam->GetInnerPadPitchWidth());
1178    } else if (GetPadType(lx) == 1) {
1179       row = (Int_t)((lx - orocLow1) / fgTPCparam->GetOuter1PadPitchLength());
1180       pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
1181    } else if (GetPadType(lx) == 2) {
1182       row = fgTPCparam->GetNRowUp1() + (Int_t)((lx - orocLow2) / fgTPCparam->GetOuter2PadPitchLength());
1183       pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
1184    }
1185    else return kFALSE;
1186    return kTRUE;
1187 }
1188
1189 void AliTPCcalibTracksGain::DumpTrack(AliTPCseed* track) {
1190    //
1191    //  Dump track information to the debug stream
1192    //
1193
1194    Int_t rows[200];
1195    Int_t sector[3];
1196    Int_t npoints[3];
1197    TVectorD dedxM[3];
1198    TVectorD dedxQ[3];
1199    TVectorD parY[3];
1200    TVectorD parZ[3];
1201    TVectorD meanPos[3];
1202    //
1203    Int_t count=0;
1204    for (Int_t ipad = 0; ipad < 3; ipad++) {
1205      dedxM[ipad].ResizeTo(5);
1206      dedxQ[ipad].ResizeTo(5);
1207      parY[ipad].ResizeTo(3);
1208      parZ[ipad].ResizeTo(3);
1209      meanPos[ipad].ResizeTo(6);
1210      Bool_t isOK = GetDedx(track, ipad, rows, sector[ipad], npoints[ipad], dedxM[ipad], dedxQ[ipad], parY[ipad], parZ[ipad], meanPos[ipad]);
1211      if (isOK)
1212        AddTracklet(sector[ipad],ipad, dedxQ[ipad], dedxM[ipad], parY[ipad], parZ[ipad], meanPos[ipad] );
1213      if (isOK) count++;
1214    }
1215
1216    TTreeSRedirector * cstream =  GetDebugStreamer();
1217    if (cstream){
1218      (*cstream) << "Track" <<
1219        "Track.=" << track <<        // track information
1220        "\n";
1221      //
1222      //
1223      //
1224      if ( GetStreamLevel()>1 && count>1){
1225        (*cstream) << "TrackG" <<
1226          "Track.=" << track <<        // track information
1227          "ncount="<<count<<
1228          // info for pad type 0
1229          "sector0="<<sector[0]<<
1230          "npoints0="<<npoints[0]<<
1231          "dedxM0.="<<&dedxM[0]<<
1232          "dedxQ0.="<<&dedxQ[0]<<
1233          "parY0.="<<&parY[0]<<
1234          "parZ0.="<<&parZ[0]<<
1235          "meanPos0.="<<&meanPos[0]<<
1236          //
1237          // info for pad type 1
1238          "sector1="<<sector[1]<<
1239          "npoints1="<<npoints[1]<<
1240          "dedxM1.="<<&dedxM[1]<<
1241          "dedxQ1.="<<&dedxQ[1]<<
1242          "parY1.="<<&parY[1]<<
1243          "parZ1.="<<&parZ[1]<<
1244          "meanPos1.="<<&meanPos[1]<<
1245          //
1246          // info for pad type 2
1247          "sector2="<<sector[2]<<
1248          "npoints2="<<npoints[2]<<
1249          "dedxM2.="<<&dedxM[2]<<
1250          "dedxQ2.="<<&dedxQ[2]<<
1251          "parY2.="<<&parY[2]<<
1252          "parZ2.="<<&parZ[2]<<
1253          "meanPos2.="<<&meanPos[2]<<
1254          //
1255          "\n";
1256      }
1257    }
1258 }
1259
1260 Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /*rows*/,
1261                                       Int_t &sector, Int_t& npoints,
1262                                       TVectorD &dedxM, TVectorD &dedxQ,
1263                                       TVectorD &parY, TVectorD &parZ, TVectorD&meanPos)
1264 {
1265   //
1266   // GetDedx for given sector for given track
1267   // padType - type of pads
1268   //
1269
1270   static TLinearFitter fitY(2, "pol1");
1271   static TLinearFitter fitZ(2, "pol1");
1272   fitY.StoreData(kFALSE);
1273   fitZ.StoreData(kFALSE);
1274   //
1275   //
1276    Int_t firstRow = 0, lastRow = 0;
1277    Int_t minRow = 100;
1278    Float_t xcenter = 0;
1279    const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
1280    const Float_t kedgey = 4.;
1281    if (padType == 0) {
1282       firstRow = 0;
1283       lastRow = fgTPCparam->GetNRowLow();
1284       xcenter = 108.47;
1285    }
1286    if (padType == 1) {
1287       firstRow = fgTPCparam->GetNRowLow();
1288       lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1289       xcenter = 166.60;
1290    }
1291    if (padType == 2) {
1292       firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1293       lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
1294       xcenter = 222.6;
1295    }
1296    minRow = (lastRow - firstRow) / 2;
1297    //
1298    //
1299    Int_t nclusters = 0;
1300    Int_t nclustersNE = 0; // number of not edge clusters
1301    Int_t lastSector = -1;
1302    Float_t amplitudeQ[100];
1303    Float_t amplitudeM[100];
1304    Int_t rowIn[100];
1305    Int_t index[100];
1306    //
1307    //
1308    fitY.ClearPoints();
1309    fitZ.ClearPoints();
1310
1311    for (Int_t iCluster = firstRow; iCluster < lastRow; iCluster++) {
1312       AliTPCclusterMI* cluster = track->GetClusterPointer(iCluster);
1313       if (cluster) {
1314          Int_t detector = cluster->GetDetector() ;
1315          if (lastSector == -1) lastSector = detector;
1316          if (lastSector != detector) continue;
1317          amplitudeQ[nclusters] = GetQNorm(cluster);
1318          amplitudeM[nclusters] = GetMaxNorm(cluster);
1319          rowIn[nclusters] = iCluster;
1320          nclusters++;
1321          Double_t dx = cluster->GetX() - xcenter;
1322          Double_t y = cluster->GetY();
1323          Double_t z = cluster->GetZ();
1324          fitY.AddPoint(&dx, y);
1325          fitZ.AddPoint(&dx, z);
1326          meanPos[0] += dx;
1327          meanPos[1] += dx;
1328          meanPos[2] += y;
1329          meanPos[3] += y*y;
1330          meanPos[4] += z;
1331          meanPos[5] += z*z;
1332          if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) nclustersNE++;
1333       }
1334    }
1335
1336    if (nclusters < minRow / 2) return kFALSE;
1337    if (nclustersNE < minRow / 2) return kFALSE;
1338    for (Int_t i = 0; i < 6; i++) meanPos[i] /= Double_t(nclusters);
1339    fitY.Eval();
1340    fitZ.Eval();
1341    fitY.GetParameters(parY);
1342    fitZ.GetParameters(parZ);
1343    //
1344    // calculate truncated mean
1345    //
1346    TMath::Sort(nclusters, amplitudeQ, index, kFALSE);
1347    //
1348    //
1349    //
1350    Float_t ndedx[5];
1351    for (Int_t i = 0; i < 5; i++) {
1352       dedxQ[i] = 0;
1353       dedxM[i] = 0;
1354       ndedx[i] = 0;
1355    }
1356    //
1357    // dedx calculation
1358    //
1359    Int_t inonEdge = 0;
1360    for (Int_t i = 0; i < nclusters; i++) {
1361       Int_t rowSorted = rowIn[index[i]];
1362       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1363
1364       if (TMath::Abs(cluster->GetY()) > cluster->GetX()*ktany - kedgey) continue;  //don't take edge clusters
1365       inonEdge++;
1366       if (inonEdge < nclustersNE * 0.5) {
1367          ndedx[0]++;
1368          dedxQ[0] += amplitudeQ[index[i]];
1369          dedxM[0] += amplitudeM[index[i]];
1370       }
1371       if (inonEdge < nclustersNE * 0.6) {
1372          ndedx[1]++;
1373          dedxQ[1] += amplitudeQ[index[i]];
1374          dedxM[1] += amplitudeM[index[i]];
1375       }
1376       if (inonEdge < nclustersNE * 0.7) {
1377          ndedx[2]++;
1378          dedxQ[2] += amplitudeQ[index[i]];
1379          dedxM[2] += amplitudeM[index[i]];
1380       }
1381       if (inonEdge < nclustersNE * 0.8) {
1382          ndedx[3]++;
1383          dedxQ[3] += amplitudeQ[index[i]];
1384          dedxM[3] += amplitudeM[index[i]];
1385       }
1386       if (inonEdge < nclustersNE * 0.9) {
1387          ndedx[4]++;
1388          dedxQ[4] += amplitudeQ[index[i]];
1389          dedxM[4] += amplitudeM[index[i]];
1390       }
1391    }
1392    for (Int_t i = 0; i < 5; i++) {
1393       dedxQ[i] /= ndedx[i];
1394       dedxM[i] /= ndedx[i];
1395    }
1396    TTreeSRedirector * cstream =  GetDebugStreamer();
1397    inonEdge = 0;
1398    Float_t momenta = track->GetP();
1399    Float_t mdedx = track->GetdEdx();
1400    for (Int_t i = 0; i < nclusters; i++) {
1401       Int_t rowSorted = rowIn[index[i]];
1402       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1403       if (!cluster) {
1404          printf("Problem\n");
1405          continue;
1406       }
1407       if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) inonEdge++;
1408       Float_t dedge = cluster->GetX()*ktany - TMath::Abs(cluster->GetY());
1409       Float_t fraction = Float_t(i) / Float_t(nclusters);
1410       Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
1411
1412       AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
1413       Float_t gain = GetGain(cluster);
1414       if (cstream) (*cstream) << "dEdx" <<
1415                      "Cl.=" << cluster <<           // cluster of interest
1416                      "gain="<<gain<<                // gain at cluster position
1417                      "P=" << momenta <<             // track momenta
1418                      "dedx=" << mdedx <<            // mean dedx - corrected for angle
1419                      "IPad=" << padType <<          // pad type 0..2
1420                      "xc=" << xcenter <<            // x center of chamber
1421                      "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
1422                      "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
1423                      "fraction=" << fraction <<     // fraction - order in statistic (0,1)
1424                      "fraction2=" << fraction2 <<   // fraction - order in statistic (0,1)
1425                      "dedge=" << dedge <<           // distance to the edge
1426                      "parY.=" << &parY <<           // line fit
1427                      "parZ.=" << &parZ <<           // line fit
1428                      "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
1429                      "\n";
1430    }
1431
1432    if (cstream) (*cstream) << "dEdxT" <<
1433      "P=" << momenta <<             // track momenta
1434      "npoints="<<inonEdge<<         // number of points
1435      "sector="<<lastSector<<        // sector number
1436      "dedx=" << mdedx <<            // mean dedx - corrected for angle
1437      "IPad=" << padType <<          // pad type 0..2
1438      "xc=" << xcenter <<            // x center of chamber
1439      "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
1440      "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
1441      "parY.=" << &parY <<           // line fit
1442      "parZ.=" << &parZ <<           // line fit
1443      "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
1444      "\n";
1445
1446    sector = lastSector;
1447    npoints = inonEdge;
1448    return kTRUE;
1449 }
1450
1451 void AliTPCcalibTracksGain::AddTracklet(UInt_t sector, UInt_t padType,TVectorD &dedxQ, TVectorD &dedxM,TVectorD& parY, TVectorD& parZ, TVectorD& meanPos){
1452   //
1453   // Add measured point - dedx to the fitter
1454   //
1455   //
1456   //chain->SetAlias("dr","(250-abs(meanPos.fElements[4]))/250");
1457   //chain->SetAlias("tz","(0+abs(parZ.fElements[1]))");
1458   //chain->SetAlias("ty","(0+abs(parY.fElements[1]))");
1459   //chain->SetAlias("corrg","sqrt((1+ty^2)*(1+tz^2))");
1460   //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);
1461
1462   Double_t xxx[100];
1463   //
1464   // z and angular part
1465   //
1466
1467   xxx[0] = (250.-TMath::Abs(meanPos[4]))/250.;
1468   xxx[1] = TMath::Abs(parY[1]);
1469   xxx[2] = TMath::Abs(parZ[1]);
1470   xxx[3] = xxx[0]*xxx[1];
1471   xxx[4] = xxx[0]*xxx[2];
1472   xxx[5] = xxx[1]*xxx[2];
1473   xxx[6] = xxx[0]*xxx[0];
1474   xxx[7] = xxx[1]*xxx[1];
1475   xxx[8] = xxx[2]*xxx[2];
1476   //
1477   // chamber part
1478   //
1479   Int_t tsector = sector%36;
1480   for (Int_t i=0;i<35;i++){
1481     xxx[9+i]=(i==tsector)?1:0;
1482   }
1483   TLinearFitter *fitterM = fFitter0M;
1484   if (padType==1) fitterM=fFitter1M;
1485   if (padType==2) fitterM=fFitter2M;
1486   fitterM->AddPoint(xxx,dedxM[1]);
1487   //
1488   TLinearFitter *fitterT = fFitter0T;
1489   if (padType==1) fitterT = fFitter1T;
1490   if (padType==2) fitterT = fFitter2T;
1491   fitterT->AddPoint(xxx,dedxQ[1]);
1492   //
1493   TLinearFitter *dfitterM = fDFitter0M;
1494   if (padType==1) dfitterM=fDFitter1M;
1495   if (padType==2) dfitterM=fDFitter2M;
1496   dfitterM->AddPoint(xxx,dedxM[1]);
1497   //
1498   TLinearFitter *dfitterT = fDFitter0T;
1499   if (padType==1) dfitterT = fDFitter1T;
1500   if (padType==2) dfitterT = fDFitter2T;
1501   dfitterT->AddPoint(xxx,dedxQ[1]);
1502 }
1503
1504
1505 TGraph *AliTPCcalibTracksGain::CreateAmpGraph(Int_t ipad, Bool_t qmax){
1506   //
1507   // create the amplitude graph
1508   // The normalized amplitudes are extrapolated to the 0 angle (y,z)  and 0 drift length
1509   //
1510
1511   TVectorD vec;
1512   if (qmax){
1513     if (ipad==0) fFitter0M->GetParameters(vec);
1514     if (ipad==1) fFitter1M->GetParameters(vec);
1515     if (ipad==2) fFitter2M->GetParameters(vec);
1516   }else{
1517     if (ipad==0) fFitter0T->GetParameters(vec);
1518     if (ipad==1) fFitter1T->GetParameters(vec);
1519     if (ipad==2) fFitter2T->GetParameters(vec);
1520   }
1521
1522   Float_t amp[36];
1523   Float_t sec[36];
1524   for (Int_t i=0;i<35;i++){
1525     sec[i]=i;
1526     amp[i]=vec[10+i]+vec[0];
1527   }
1528   amp[35]=vec[0];
1529   Float_t mean = TMath::Mean(36,amp);
1530   for (Int_t i=0;i<36;i++){
1531     sec[i]=i;
1532     amp[i]=(amp[i]-mean)/mean;
1533   }
1534   TGraph *gr = new TGraph(36,sec,amp);
1535   return gr;
1536 }
1537
1538
1539 void   AliTPCcalibTracksGain::UpdateClusterParam(AliTPCClusterParam* clparam){
1540   //
1541   //   SetQ normalization parameters
1542   //
1543   //  void SetQnorm(Int_t ipad, Int_t itype,  TVectorD * norm);
1544
1545   TVectorD vec;
1546
1547   //
1548   fDFitter0T->Eval();
1549   fDFitter1T->Eval();
1550   fDFitter2T->Eval();
1551   fDFitter0M->Eval();
1552   fDFitter1M->Eval();
1553   fDFitter2M->Eval();
1554   fDFitter0T->GetParameters(vec);
1555   clparam->SetQnorm(0,0,&vec);
1556   fDFitter1T->GetParameters(vec);
1557   clparam->SetQnorm(1,0,&vec);
1558   fDFitter2T->GetParameters(vec);
1559   clparam->SetQnorm(2,0,&vec);
1560   //
1561   fDFitter0M->GetParameters(vec);
1562   clparam->SetQnorm(0,1,&vec);
1563   fDFitter1M->GetParameters(vec);
1564   clparam->SetQnorm(1,1,&vec);
1565   fDFitter2M->GetParameters(vec);
1566   clparam->SetQnorm(2,1,&vec);
1567   //
1568
1569 }
1570
1571
1572 void   AliTPCcalibTracksGain::Analyze(){
1573
1574  Evaluate();
1575
1576 }
1577
1578
1579
1580 TVectorD *  AliTPCcalibTracksGain::MakeQPosNorm(TTree * chain0, Int_t ipad, Bool_t isMax, Int_t maxPoints, Int_t verbose){
1581   //
1582   // Input parameters
1583   // chain0 - the tree with information  -Debug stream
1584   // ipad   - 0 IROC
1585   //        - 1 OROC medium
1586   //        - 2 OROC LONG
1587   // isMax  - kFALSE - total charge param
1588   //          kTRUE  - Max charge param
1589   //
1590   // maxPoints - number of points for fit
1591   //
1592   // verbose -
1593   //
1594   /* e.g
1595     ipad=0
1596     isMax=kTRUE;
1597     maxPoints=1000000;
1598   */
1599   // Make Q normalization as function of following parameters
1600   // 1 - dp   - relative pad position
1601   // 2 - dt   - relative time position
1602   // 3 - di   - drift length (norm to 1);
1603   // 4 - dq0  - Tot/Max charge
1604   // 5 - dq1  - Max/Tot charge
1605   // 6 - sy   - sigma y - shape
1606   // 7 - sz   - sigma z - shape
1607   //
1608   // Coeficient of Taylor expansion fitted
1609   // Fit parameters returned as TVectorD
1610   // Fit parameters to be used in corresponding correction function
1611   // in AliTPCclusterParam
1612   //
1613   //
1614   TStatToolkit toolkit;
1615   Double_t chi2;
1616   TVectorD fitParam;
1617   TMatrixD covMatrix;
1618   Int_t npoints;
1619   TCut cutA("dedge>3&&fraction2<0.7");
1620   chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1621   chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1622   chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1623   chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1624   chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1625   chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1626   chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1627   //
1628   TString fstring="";
1629   fstring+="dp++";                               //1
1630   fstring+="dt++";                               //2
1631   fstring+="dp*dp++";                             //3
1632   fstring+="dt*dt++";                             //4
1633   fstring+="dt*dt*dt++";                             //5
1634   fstring+="dp*dt++";                            //6
1635   fstring+="dp*dt*dt++";                          //7
1636   fstring+="(dq0)++";                            //8
1637   fstring+="(dq1)++";                            //9
1638   //
1639   //
1640   fstring+="dp*dp*(di)++";                        //10
1641   fstring+="dt*dt*(di)++";                        //11
1642   fstring+="dp*dp*sy++";                          //12
1643   fstring+="dt*sz++";                          //13
1644   fstring+="dt*dt*sz++";                          //14
1645   fstring+="dt*dt*dt*sz++";                          //15
1646   //
1647   fstring+="dp*dp*1*sy*sz++";                     //16
1648   fstring+="dt*sy*sz++";                       //17
1649   fstring+="dt*dt*sy*sz++";                       //18
1650   fstring+="dt*dt*dt*sy*sz++";                       //19
1651   //
1652   fstring+="dp*dp*(dq0)++";                       //20
1653   fstring+="dt*1*(dq0)++";                       //21
1654   fstring+="dt*dt*(dq0)++";                       //22
1655   fstring+="dt*dt*dt*(dq0)++";                       //23
1656   //
1657   fstring+="dp*dp*(dq1)++";                       //24
1658   fstring+="dt*(dq1)++";                       //25
1659   fstring+="dt*dt*(dq1)++";                       //26
1660   fstring+="dt*dt*dt*(dq1)++";                       //27
1661
1662   TString var;
1663   if (isMax)  var = "Cl.fMax/gain/dedxM.fElements[2]";
1664   if (!isMax) var = "Cl.fQ/gain/dedxQ.fElements[2]";
1665   TString cutP="IPad==";
1666   cutP+=ipad;
1667   //
1668   TString *strq0 = toolkit.FitPlane(chain0,var.Data(),fstring.Data(), cutP.Data()+cutA, chi2,npoints,fitParam,covMatrix,-1,0,maxPoints);
1669   //
1670   //
1671   if (verbose){
1672     printf("Chi2/npoints = %f",TMath::Sqrt(chi2/npoints));
1673     printf("\nFit function\n:%s\n",strq0->Data());
1674   }
1675   TVectorD *vec = new TVectorD(fitParam);
1676   return vec;
1677 }
1678
1679 void  AliTPCcalibTracksGain::MakeQPosNormAll(TTree * chain, AliTPCClusterParam * param, Int_t maxPoints, Int_t verbose){
1680   //
1681   // Fill the content of the of the AliTPCclusterParam
1682   // with fitted values of corrections
1683   //
1684
1685 }
1686
1687
1688
1689 /*
1690
1691  Position correction fit:
1692  //
1693 TStatToolkit toolkit;
1694 Double_t chi2;
1695 TVectorD fitParam;
1696 TMatrixD covMatrix;
1697 Int_t npoints;
1698 //
1699 TCut cutA("dedge>3&&fraction2<0.7");
1700 chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1701 chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1702 chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1703 chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1704 chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1705 chain0->SetAlias("sy","(0.2/sqrt(0.01^2+Cl.fSigmaY2))");
1706 chain0->SetAlias("sz","(0.2/sqrt(0.01^2+Cl.fSigmaZ2))");
1707
1708 TString fstring="";
1709
1710 fstring+="dp++";                               //1
1711 fstring+="dt++";                               //2
1712 fstring+="dp*dp++";                             //3
1713 fstring+="dt*dt++";                             //4
1714 fstring+="dt*dt*dt++";                             //5
1715 fstring+="dp*dt++";                            //6
1716 fstring+="dp*dt*dt++";                          //7
1717 fstring+="(dq0)++";                            //8
1718 fstring+="(dq1)++";                            //9
1719 //
1720 //
1721 fstring+="dp*dp*(di)++";                        //10
1722 fstring+="dt*dt*(di)++";                        //11
1723 fstring+="dp*dp*sy++";                          //12
1724 fstring+="dt*sz++";                          //13
1725 fstring+="dt*dt*sz++";                          //14
1726 fstring+="dt*dt*dt*sz++";                          //15
1727 //
1728 fstring+="dp*dp*1*sy*sz++";                     //16
1729 fstring+="dt*sy*sz++";                       //17
1730 fstring+="dt*dt*sy*sz++";                       //18
1731 fstring+="dt*dt*dt*sy*sz++";                       //19
1732 //
1733 fstring+="dp*dp*(dq0)++";                       //20
1734 fstring+="dt*1*(dq0)++";                       //21
1735 fstring+="dt*dt*(dq0)++";                       //22
1736 fstring+="dt*dt*dt*(dq0)++";                       //23
1737
1738 fstring+="dp*dp*(dq1)++";                       //24
1739 fstring+="dt*(dq1)++";                       //25
1740 fstring+="dt*dt*(dq1)++";                       //26
1741 fstring+="dt*dt*dt*(dq1)++";                       //27
1742
1743
1744  TString *strq0 = toolkit.FitPlane(chain0,"Cl.fMax/gain/dedxM.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
1745  TString *strqt0 = toolkit.FitPlane(chain0,"Cl.fQ/gain/dedxQ.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
1746
1747  chain0->SetAlias("qcorM0",strq0->Data());
1748  chain0->SetAlias("qcorT0",strqt0->Data());
1749 //chain0->SetAlias("mmqcorM0","min(max(qcorM0,0.75),1.15)");
1750  chain0->Draw("(Cl.fMax/gain/dedxM.fElements[2]):min(max(qcorM0,0.75),1.15)","IPad==0"+cutA,"prof",100000)
1751
1752  fraction05 -
1753      sigma                  0.2419
1754      sigma fit              0.2302
1755      sigma fit with shape   0.2257
1756  fraction 07
1757      qtot sigma                  0.322
1758      qmax sigma                  0.292
1759      qmax sigma fit              0.2702
1760      qmax sigma fit+ratio        0.2638
1761
1762 */
1763