TPC/Calib/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 #include "AliTPCcalibTracksGain.h"
 170
 171
 172 #include <TPDGCode.h>
 173 #include <TStyle.h>
 174 #include "TSystem.h"
 175 #include "TMatrixD.h"
 176 #include "TTreeStream.h"
 177 #include "TF1.h"
 178 #include "AliTPCParamSR.h"
 179 #include "AliTPCClusterParam.h"
 180 #include "AliTrackPointArray.h"
 181 #include <TH1.h>
 182 #include <TH3F.h>
 183 #include <TLinearFitter.h>
 184 #include <TTreeStream.h>
 185 #include <TFile.h>
 186 #include <TCollection.h>
 187 #include <TIterator.h>
 188 #include <TProfile.h>
 189 #include <TProfile2D.h>
 190 #include <TProof.h>
 191 #include <TStatToolkit.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 #include "AliTPCcalibDB.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("libPWGPP.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   //
 248   // Fitters
 249   //
 250   fSimpleFitter(0),         // simple fitter for short pads
 251   fSqrtFitter(0),           // sqrt fitter for medium pads
 252   fLogFitter(0),            // log fitter for long pads
 253
 254   fFitter0M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 255   fFitter1M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 256   fFitter2M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 257   fFitter0T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 258   fFitter1T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 259   fFitter2T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 260   //
 261   fDFitter0M(0),              // fitting of the atenuation, angular correction
 262   fDFitter1M(0),              // fitting of the atenuation, angular correction
 263   fDFitter2M(0),              // fitting of the atenuation, angular correction
 264   fDFitter0T(0),              // fitting of the atenuation, angular correction
 265   fDFitter1T(0),              // fitting of the atenuation, angular correction
 266   fDFitter2T(0),              // fitting of the atenuation, angular correction
 267   //
 268   fSingleSectorFitter(0),   // just for debugging
 269   //
 270   // Counters
 271   //
 272   fTotalTracks(0),         // just for debugging
 273   fAcceptedTracks(0)      // just for debugging
 274
 275 {
 276    //
 277    // Default constructor.
 278    //
 279 }
 280
 281 AliTPCcalibTracksGain::AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj) :
 282   AliTPCcalibBase(obj),
 283   fCuts(obj.fCuts),            // cuts that are used for sieving the tracks used for calibration
 284   //
 285   // Fitters
 286   //
 287   fSimpleFitter(obj.fSimpleFitter),         // simple fitter for short pads
 288   fSqrtFitter(obj.fSqrtFitter),           // sqrt fitter for medium pads
 289   fLogFitter(obj.fLogFitter),            // log fitter for long pads
 290   fFitter0M(obj.fFitter0M),
 291   fFitter1M(obj.fFitter1M),
 292   fFitter2M(obj.fFitter2M),
 293   fFitter0T(obj.fFitter0T),
 294   fFitter1T(obj.fFitter1T),
 295   fFitter2T(obj.fFitter2T),
 296   //
 297   fDFitter0M(obj.fDFitter0M),
 298   fDFitter1M(obj.fDFitter1M),
 299   fDFitter2M(obj.fDFitter2M),
 300   fDFitter0T(obj.fDFitter0T),
 301   fDFitter1T(obj.fDFitter1T),
 302   fDFitter2T(obj.fDFitter2T),
 303   fSingleSectorFitter(obj.fSingleSectorFitter),   // just for debugging
 304   //
 305   // Counters
 306   //
 307   fTotalTracks(obj.fTotalTracks),         // just for debugging
 308   fAcceptedTracks(obj.fAcceptedTracks)      // just for debugging
 309
 310 {
 311    //
 312    // Copy constructor.
 313    //
 314 }
 315
 316 AliTPCcalibTracksGain& AliTPCcalibTracksGain::operator=(const AliTPCcalibTracksGain& rhs) {
 317    //
 318    // Assignment operator.
 319    //
 320
 321    if (this != &rhs) {
 322       TNamed::operator=(rhs);
 323       fSimpleFitter = new AliTPCFitPad(*(rhs.fSimpleFitter));
 324       fSqrtFitter = new AliTPCFitPad(*(rhs.fSqrtFitter));
 325       fLogFitter = new AliTPCFitPad(*(rhs.fLogFitter));
 326       fSingleSectorFitter = new AliTPCFitPad(*(rhs.fSingleSectorFitter));
 327       fCuts = new AliTPCcalibTracksCuts(*(rhs.fCuts));
 328    }
 329    return *this;
 330 }
 331
 332 AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts) :
 333   AliTPCcalibBase(),
 334   fCuts(0),            // cuts that are used for sieving the tracks used for calibration
 335   //
 336   // Fitters
 337   //
 338   fSimpleFitter(0),         // simple fitter for short pads
 339   fSqrtFitter(0),           // sqrt fitter for medium pads
 340   fLogFitter(0),            // log fitter for long pads
 341   fFitter0M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 342   fFitter1M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 343   fFitter2M(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 344   fFitter0T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 345   fFitter1T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 346   fFitter2T(0),              // fitting of the atenuation, angular correction, and mean chamber gain
 347   //
 348   fDFitter0M(0),              // fitting of the atenuation, angular correction
 349   fDFitter1M(0),              // fitting of the atenuation, angular correction
 350   fDFitter2M(0),              // fitting of the atenuation, angular correction
 351   fDFitter0T(0),              // fitting of the atenuation, angular correction
 352   fDFitter1T(0),              // fitting of the atenuation, angular correction
 353   fDFitter2T(0),              // fitting of the atenuation, angular correction
 354   fSingleSectorFitter(0),   // just for debugging
 355   //
 356   // Counters
 357   //
 358   fTotalTracks(0),         // just for debugging
 359   fAcceptedTracks(0)      // just for debugging
 360
 361 {
 362    //
 363    // Constructor.
 364    //
 365    this->SetNameTitle(name, title);
 366    fCuts = cuts;
 367    //
 368    // Fitter initialization
 369    //
 370    fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
 371    fSqrtFitter   = new AliTPCFitPad(8, "hyp7", "");
 372    fLogFitter    = new AliTPCFitPad(8, "hyp7", "");
 373    fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
 374    //
 375    fFitter0M      = new TLinearFitter(45,"hyp44");
 376    fFitter1M      = new TLinearFitter(45,"hyp44");
 377    fFitter2M      = new TLinearFitter(45,"hyp44");
 378    fFitter0T      = new TLinearFitter(45,"hyp44");
 379    fFitter1T      = new TLinearFitter(45,"hyp44");
 380    fFitter2T      = new TLinearFitter(45,"hyp44");
 381    //
 382    fDFitter0M      = new TLinearFitter(10,"hyp9");
 383    fDFitter1M      = new TLinearFitter(10,"hyp9");
 384    fDFitter2M      = new TLinearFitter(10,"hyp9");
 385    fDFitter0T      = new TLinearFitter(10,"hyp9");
 386    fDFitter1T      = new TLinearFitter(10,"hyp9");
 387    fDFitter2T      = new TLinearFitter(10,"hyp9");
 388    //
 389    //
 390    fFitter0M->StoreData(kFALSE);
 391    fFitter1M->StoreData(kFALSE);
 392    fFitter2M->StoreData(kFALSE);
 393    fFitter0T->StoreData(kFALSE);
 394    fFitter1T->StoreData(kFALSE);
 395    fFitter2T->StoreData(kFALSE);
 396    //
 397    fDFitter0M->StoreData(kFALSE);
 398    fDFitter1M->StoreData(kFALSE);
 399    fDFitter2M->StoreData(kFALSE);
 400    fDFitter0T->StoreData(kFALSE);
 401    fDFitter1T->StoreData(kFALSE);
 402    fDFitter2T->StoreData(kFALSE);
 403    //
 404    //
 405    // just for debugging -counters
 406    //
 407    fTotalTracks     = 0;
 408    fAcceptedTracks  = 0;
 409 }
 410
 411 AliTPCcalibTracksGain::~AliTPCcalibTracksGain() {
 412    //
 413    // Destructor.
 414    //
 415
 416    Info("Destructor",":");
 417    if (fSimpleFitter) delete fSimpleFitter;
 418    if (fSqrtFitter) delete fSqrtFitter;
 419    if (fLogFitter) delete fLogFitter;
 420    if (fSingleSectorFitter) delete fSingleSectorFitter;
 421
 422 }
 423
 424
 425
 426
 427 void AliTPCcalibTracksGain::Terminate(){
 428    //
 429    // Evaluate fitters and close the debug stream.
 430    // Also move or copy the debug stream, if a debugStreamPrefix is provided.
 431    //
 432
 433    Evaluate();
 434    AliTPCcalibBase::Terminate();
 435 }
 436
 437
 438
 439 void AliTPCcalibTracksGain::Process(AliTPCseed* seed) {
 440    //
 441    // Main method to be called when a new seed is supposed to be processed
 442    // and be used for gain calibration. Its quality is checked before it
 443    // is added.
 444    //
 445
 446
 447    fTotalTracks++;
 448    if (!fCuts->AcceptTrack(seed)) return;
 449    //
 450    // reinint on proof
 451    //   if (gProof){
 452      static Bool_t doinit= kTRUE;
 453      if (doinit){
 454        fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
 455        fSqrtFitter   = new AliTPCFitPad(8, "hyp7", "");
 456        fLogFitter    = new AliTPCFitPad(8, "hyp7", "");
 457        fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
 458        //
 459        fFitter0M      = new TLinearFitter(45,"hyp44");
 460        fFitter1M      = new TLinearFitter(45,"hyp44");
 461        fFitter2M      = new TLinearFitter(45,"hyp44");
 462        fFitter0T      = new TLinearFitter(45,"hyp44");
 463        fFitter1T      = new TLinearFitter(45,"hyp44");
 464        fFitter2T      = new TLinearFitter(45,"hyp44");
 465        //
 466        fDFitter0M      = new TLinearFitter(10,"hyp9");
 467        fDFitter1M      = new TLinearFitter(10,"hyp9");
 468        fDFitter2M      = new TLinearFitter(10,"hyp9");
 469        fDFitter0T      = new TLinearFitter(10,"hyp9");
 470        fDFitter1T      = new TLinearFitter(10,"hyp9");
 471        fDFitter2T      = new TLinearFitter(10,"hyp9");
 472        doinit=kFALSE;
 473      }
 474      //}
 475
 476    fAcceptedTracks++;
 477    AddTrack(seed);
 478 }
 479
 480 Long64_t AliTPCcalibTracksGain::Merge(TCollection *list) {
 481    //
 482    // Merge() merges the results of all AliTPCcalibTracksGain objects contained in
 483    // list, thus allowing a distributed computation of several files, e.g. on PROOF.
 484    // The merged results are merged with the data members of the AliTPCcalibTracksGain
 485    // object used for calling the Merge method.
 486    // The return value is 0 /*the total number of tracks used for calibration*/ if the merge
 487    // is successful, otherwise it is -1.
 488    //
 489
 490    if (!list || list->IsEmpty()) return -1;
 491
 492    if (!fSimpleFitter) fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
 493    if (!fSqrtFitter)   fSqrtFitter   = new AliTPCFitPad(8, "hyp7", "");
 494    if (!fLogFitter)    fLogFitter    = new AliTPCFitPad(8, "hyp7", "");
 495    if (!fSingleSectorFitter) fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
 496
 497
 498
 499    TIterator* iter = list->MakeIterator();
 500    AliTPCcalibTracksGain* cal = 0;
 501
 502    while ((cal = (AliTPCcalibTracksGain*)iter->Next())) {
 503       if (!cal->InheritsFrom(AliTPCcalibTracksGain::Class())) {
 504         //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 505          return -1;
 506       }
 507
 508       Add(cal);
 509    }
 510    return 0;
 511 }
 512
 513 Float_t  AliTPCcalibTracksGain::GetGain(AliTPCclusterMI *cluster){
 514   //
 515   // get gain
 516   //
 517   Float_t gainPad= 1;
 518   AliTPCCalPad * gainMap =  AliTPCcalibDB::Instance()->GetDedxGainFactor();
 519   if (gainMap) {
 520     AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector());
 521     gainPad  = roc->GetValue(cluster->GetRow(), TMath::Nint(cluster->GetPad()));
 522   }
 523   return gainPad;
 524 }
 525
 526 Float_t   AliTPCcalibTracksGain::GetMaxNorm(AliTPCclusterMI * cluster, Float_t ky, Float_t kz){
 527   //
 528   // Get normalized amplituded if the gain map provided
 529   //
 530   AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
 531   Float_t maxNorm =
 532     parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),
 533                           cluster->GetTimeBin(),ky,kz,0.5,0.5,1.6);
 534
 535   return GetGain(cluster)*maxNorm;
 536 }
 537
 538
 539 Float_t   AliTPCcalibTracksGain::GetQNorm(AliTPCclusterMI * cluster,  Float_t ky, Float_t kz){
 540   //
 541   // Get normalized amplituded if the gain map provided
 542   //
 543   AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
 544   Float_t totNorm = parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),                        cluster->GetTimeBin(),ky,kz,0.5,0.5,cluster->GetQ(),2.5,1.6);
 545   return GetGain(cluster)*totNorm;
 546 }
 547
 548
 549
 550 void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
 551    //
 552    // Adds another AliTPCcalibTracksGain object to this object.
 553    //
 554
 555   fSimpleFitter->Add(cal->fSimpleFitter);
 556   fSqrtFitter->Add(cal->fSqrtFitter);
 557   fLogFitter->Add(cal->fLogFitter);
 558   fSingleSectorFitter->Add(cal->fSingleSectorFitter);
 559   //
 560   //
 561   //
 562   if (cal->fFitter0M->GetNpoints()>0) fFitter0M->Add(cal->fFitter0M);
 563   if (cal->fFitter1M->GetNpoints()>0) fFitter1M->Add(cal->fFitter1M);
 564   if (cal->fFitter2M->GetNpoints()>0) fFitter2M->Add(cal->fFitter2M);
 565   if (cal->fFitter0T->GetNpoints()>0) fFitter0T->Add(cal->fFitter0T);
 566   if (cal->fFitter1T->GetNpoints()>0) fFitter1T->Add(cal->fFitter1T);
 567   if (cal->fFitter2T->GetNpoints()>0) fFitter2T->Add(cal->fFitter2T);
 568   //
 569   if (cal->fDFitter0M->GetNpoints()>0) fDFitter0M->Add(cal->fDFitter0M);
 570   if (cal->fDFitter1M->GetNpoints()>0) fDFitter1M->Add(cal->fDFitter1M);
 571   if (cal->fDFitter2M->GetNpoints()>0) fDFitter2M->Add(cal->fDFitter2M);
 572   if (cal->fDFitter0T->GetNpoints()>0) fDFitter0T->Add(cal->fDFitter0T);
 573   if (cal->fDFitter1T->GetNpoints()>0) fDFitter1T->Add(cal->fDFitter1T);
 574   if (cal->fDFitter2T->GetNpoints()>0) fDFitter2T->Add(cal->fDFitter2T);
 575    //
 576
 577    // just for debugging, remove me
 578    fTotalTracks += cal->fTotalTracks;
 579    fAcceptedTracks += cal->fAcceptedTracks;
 580
 581 }
 582
 583 void AliTPCcalibTracksGain::AddTrack(AliTPCseed* seed) {
 584    //
 585    // The clusters making up the track (seed) are added to various fit functions.
 586    // See AddCluster(...) for more detail.
 587    //
 588
 589    DumpTrack(seed);
 590 }
 591
 592
 593
 594
 595 void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t /*momenta*/, Float_t/* mdedx*/, Int_t padType,
 596                                        Float_t xcenter, TVectorD& dedxQ, TVectorD& /*dedxM*/, Float_t /*fraction*/, Float_t fraction2, Float_t dedge,
 597                                        TVectorD& parY, TVectorD& parZ, TVectorD& meanPos) {
 598    //
 599    // Adds cluster to the appropriate fitter for later analysis.
 600    // The charge used for the fit is the maximum charge for this specific cluster or the
 601    // accumulated charge per cluster, depending on the value of fgkUseTotalCharge.
 602    // Depending on the pad size where the cluster is registered, the value will be put in
 603    // the appropriate fitter. Furthermore, for each pad size three different types of fitters
 604    // are used. The fit functions are the same for all fitters (parabolic functions), but the value
 605    // added to each fitter is different. The simple fitter gets the charge plugged in as is, the sqrt fitter
 606    // gets the square root of the charge, and the log fitter gets fgkM*(1+q/fgkM), where q is the original charge
 607    // and fgkM==25.
 608    //
 609   Float_t kedge     = 3;
 610   Float_t kfraction = 0.7;
 611   Int_t   kinorm    = 2;
 612
 613
 614   // Where to put selection on threshold?
 615   // Defined by the Q/dEdxT variable - see debug streamer:
 616   //
 617   // Debug stream variables:  (Where tu cut ?)
 618   // chain0->Draw("Cl.fQ/dedxQ.fElements[1]>>his(100,0,3)","fraction2<0.6&&dedge>3","",1000000);
 619   //         mean 1  sigma 0.25
 620   // chain0->Draw("Cl.fMax/dedxM.fElements[1]>>his(100,0,3)","fraction2<0.6&&dedge>3","",1000000)
 621   //         mean 1  sigma 0.25
 622   // chain0->Draw("Cl.fQ/dedxQ.fElements[2]>>his(100,0,3)","fraction2<0.7&&dedge>3","",1000000)
 623   //         mean 1 sigma 0.29
 624   // chain0->Draw("Cl.fMax/dedxM.fElements[2]>>his(100,0,3)","fraction2<0.7&&dedge>3","",1000000)
 625   //         mean 1 sigma 0.27
 626   // chain0->Draw("Cl.fQ/dedxQ.fElements[3]>>his(100,0,3)","fraction2<0.8&&dedge>3","",1000000)
 627   //         mean 1 sigma 0.32
 628   //
 629   // chain0->Draw("Cl.fQ/dedxQ.fElements[4]>>his(100,0,3)","fraction2<0.9&&dedge>3","",1000000)
 630   //         mean 1 sigma 0.4
 631
 632   // Fraction choosen 0.7
 633
 634    if (!cluster) {
 635       Error("AddCluster", "Cluster not valid.");
 636       return;
 637    }
 638
 639    if (dedge < kedge) return;
 640    if (fraction2 > kfraction) return;
 641
 642    //Int_t padType = GetPadType(cluster->GetX());
 643    Double_t xx[7];
 644    //Double_t centerPad[2] = {0};
 645    //AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
 646    //xx[0] = cluster->GetX() - centerPad[0];
 647    //xx[1] = cluster->GetY() - centerPad[1];
 648    xx[0] = cluster->GetX() - xcenter;
 649    xx[1] = cluster->GetY();
 650    xx[2] = xx[0] * xx[0];
 651    xx[3] = xx[1] * xx[1];
 652    xx[4] = xx[0] * xx[1];
 653    xx[5] = TMath::Abs(cluster->GetZ()) - TMath::Abs(meanPos[4]);
 654    xx[6] = xx[5] * xx[5];
 655
 656    //
 657    // Update fitters
 658    //
 659    Int_t segment = cluster->GetDetector() % 36;
 660    Double_t q = fgkUseTotalCharge ?
 661      ((Double_t)(cluster->GetQ()/GetQNorm(cluster,parY[1], parZ[1]))) : ((Double_t)(cluster->GetMax()/GetMaxNorm(cluster,parY[1], parZ[1])));
 662
 663    // correct charge by normalising to mean charge per track
 664    q /= dedxQ[kinorm];
 665
 666    // just for debugging
 667
 668    Double_t sqrtQ = TMath::Sqrt(q);
 669    Double_t logQ = fgkM * TMath::Log(1 + q / fgkM);
 670    TLinearFitter * fitter =0;
 671    //
 672    fitter = fSimpleFitter->GetFitter(segment, padType);
 673    fitter->AddPoint(xx, q);
 674    //
 675    fitter = fSqrtFitter->GetFitter(segment, padType);
 676    fitter->AddPoint(xx, sqrtQ);
 677    //
 678    fitter = fLogFitter->GetFitter(segment, padType);
 679    fitter->AddPoint(xx, logQ);
 680    //
 681    fitter=fSingleSectorFitter->GetFitter(0, padType);
 682    fitter->AddPoint(xx, q);
 683
 684 }
 685
 686 void AliTPCcalibTracksGain::Evaluate(Bool_t robust, Double_t frac) {
 687    //
 688    // Evaluates all fitters contained in this object.
 689    // If the robust option is set to kTRUE a robust fit is performed with frac as
 690    // the minimal fraction of good points (see TLinearFitter::EvalRobust for details).
 691    // Beware: Robust fitting is much slower!
 692    //
 693
 694    fSimpleFitter->Evaluate(robust, frac);
 695    fSqrtFitter->Evaluate(robust, frac);
 696    fLogFitter->Evaluate(robust, frac);
 697    fSingleSectorFitter->Evaluate(robust, frac);
 698    fFitter0M->Eval();
 699    fFitter1M->Eval();
 700    fFitter2M->Eval();
 701    fFitter0T->Eval();
 702    fFitter1T->Eval();
 703    fFitter2T->Eval();
 704    //
 705    fDFitter0M->Eval();
 706    fDFitter1M->Eval();
 707    fDFitter2M->Eval();
 708    fDFitter0T->Eval();
 709    fDFitter1T->Eval();
 710    fDFitter2T->Eval();
 711 }
 712
 713 AliTPCCalPad* AliTPCcalibTracksGain::CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 714   //
 715   // Creates the calibration object AliTPCcalPad using fitted parameterization
 716   //
 717    TObjArray tpc(72);
 718    for (UInt_t iSector = 0; iSector < 72; iSector++)
 719       tpc.Add(CreateFitCalROC(iSector, fitType, undoTransformation, normalizeToPadSize));
 720    return new AliTPCCalPad(&tpc);
 721 }
 722
 723 AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 724   //
 725   // Create the AliTPCCalROC with the values per pad
 726   // sector  - sector of interest
 727   // fitType -
 728   //
 729
 730    TVectorD par(8);
 731    if (sector < 36) {
 732       GetParameters(sector % 36, 0, fitType, par);
 733       return CreateFitCalROC(sector, 0, par, fitType, undoTransformation, normalizeToPadSize);
 734    }
 735    else {
 736       GetParameters(sector % 36, 1, fitType, par);
 737       AliTPCCalROC* roc1 = CreateFitCalROC(sector, 1, par, fitType, undoTransformation, normalizeToPadSize);
 738       GetParameters(sector % 36, 2, fitType, par);
 739       AliTPCCalROC* roc2 = CreateFitCalROC(sector, 2, par, fitType, undoTransformation, normalizeToPadSize);
 740       AliTPCCalROC* roc3 = CreateCombinedCalROC(roc1, roc2);
 741       delete roc1;
 742       delete roc2;
 743       return roc3;
 744    }
 745 }
 746
 747 AliTPCCalROC* AliTPCcalibTracksGain::CreateFitCalROC(UInt_t sector, UInt_t padType, TVectorD &fitParam, UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
 748    //
 749    // This function is essentially a copy of AliTPCCalROC::CreateGlobalFitCalROC(...), with the
 750    // modifications, that the center of the region of same pad size is used as the origin
 751    // of the fit function instead of the center of the ROC.
 752    // The possibility of a linear fit is removed as well because it is not needed.
 753    // Only values for pads with the given pad size are calculated, the rest is 0.
 754    // Set undoTransformation for undoing the transformation that was applied to the
 755    // charge values before they were put into the fitter (thus allowing comparison to the original
 756    // charge values). For fitType use 0 for the simple fitter, 1 for the sqrt fitter, 2 for the log fitter.
 757    // If normalizeToPadSize is true, the values are normalized to the pad size.
 758    // Please be aware, that you even need to specify the fitType if you want to normalize to the pad size without
 759    // undoing the transformation (because normalizing involves undoing the trafo first, then normalizing, then
 760    // applying the trafo again).
 761    // Please note: The normalization to the pad size is a simple linear scaling with the pad length, which
 762    //              actually doesn't describe reality!
 763    //
 764
 765    Float_t dlx, dly;
 766    Double_t centerPad[2] = {0};
 767    Float_t localXY[3] = {0};
 768    AliTPCROC* tpcROC = AliTPCROC::Instance();
 769    if ((padType == 0 && sector >= tpcROC->GetNInnerSector()) || (padType > 0 && sector < tpcROC->GetNInnerSector()) || sector >= tpcROC->GetNSector())
 770       return 0;
 771    AliTPCCalROC* lROCfitted = new AliTPCCalROC(sector);
 772    //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
 773    UInt_t startRow = 0;
 774    UInt_t endRow = 0;
 775    switch (padType) {
 776       case kShortPads:
 777          startRow = 0;
 778          endRow = lROCfitted->GetNrows();
 779          break;
 780       case kMediumPads:
 781          startRow = 0;
 782          endRow = 64;
 783          break;
 784       case kLongPads:
 785          startRow = 64;
 786          endRow = lROCfitted->GetNrows();
 787          break;
 788    }
 789
 790    AliTPCFitPad::GetPadRegionCenterLocal(padType, centerPad);
 791    Double_t value = 0;
 792    for (UInt_t irow = startRow; irow < endRow; irow++) {
 793       for (UInt_t ipad = 0; ipad < lROCfitted->GetNPads(irow); ipad++) {
 794          tpcROC->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 795          dlx = localXY[0] - centerPad[0];
 796          dly = localXY[1] - centerPad[1];
 797          value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
 798
 799          // Let q' = value be the transformed value without any pad size corrections,
 800          // let T be the transformation and let l be the pad size
 801          //    1) don't undo transformation, don't normalize: return q'
 802          //    2) undo transformation,       don't normalize: return T^{-1} q'
 803          //    3) undo transformation,       normalize:       return (T^{-1} q') / l
 804          //    4) don't undo transformation, normalize:       return T((T^{-1} q') / l)
 805          if (!undoTransformation && !normalizeToPadSize) {/* value remains unchanged */}  // (1)
 806          else {                                                                           // (2), (3), (4)
 807             //calculate T^{-1}
 808             switch (fitType) {
 809                case  0: /* value remains unchanged */ break;
 810                case  1: value = value * value; break;
 811                case  2: value = (TMath::Exp(value / fgkM) - 1) * fgkM; break;
 812                default: Error("CreateFitCalROC", "Wrong fit type."); break;
 813             }
 814             if (normalizeToPadSize) value /= GetPadLength(localXY[0]);                    // (3)
 815          }
 816          if (!undoTransformation && normalizeToPadSize) {                                 // (4)
 817             // calculate T
 818             switch (fitType) {
 819                case  0: /* value remains unchanged */ break;
 820                case  1: value = TMath::Sqrt(value); break;
 821                case  2: value = fgkM * TMath::Log(1 + value / fgkM); break;
 822                default: Error("CreateFitCalROC", "Wrong fit type."); break;
 823             }
 824          }
 825          lROCfitted->SetValue(irow, ipad, value);
 826       }
 827    }
 828    return lROCfitted;
 829 }
 830
 831 AliTPCCalROC* AliTPCcalibTracksGain::CreateCombinedCalROC(const AliTPCCalROC* roc1, const AliTPCCalROC* roc2) {
 832    //
 833    // Combines the medium pad size values of roc1 with the long pad size values of roc2 into a new
 834    // AliTPCCalROC. Returns a null pointer if any one of the ROCs is an IROC; issues a warning message
 835    // if the sectors of roc1 and roc2 don't match, but still continue and use the sector of roc1 as the
 836    // sector of the new ROC.
 837    //
 838
 839    if (!roc1 || !roc2) return 0;
 840    if ((Int_t)(roc1->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
 841    if ((Int_t)(roc2->GetSector()) < fgTPCparam->GetNInnerSector()) return 0;
 842    if (roc1->GetSector() != roc2->GetSector()) Warning("CreateCombinedCalROC", "Sector number mismatch.");
 843    AliTPCCalROC* roc = new AliTPCCalROC(roc1->GetSector());
 844
 845    for (UInt_t iRow = 0; iRow < 64; iRow++) {
 846       for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
 847          roc->SetValue(iRow, iPad, roc1->GetValue(iRow, iPad));
 848    }
 849    for (UInt_t iRow = 64; iRow < roc->GetNrows(); iRow++) {
 850       for (UInt_t iPad = 0; iPad < roc->GetNPads(iRow); iPad++)
 851          roc->SetValue(iRow, iPad, roc2->GetValue(iRow, iPad));
 852    }
 853    return roc;
 854 }
 855
 856 Bool_t AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) {
 857    //
 858    // Puts the fit parameters for the specified segment (IROC & OROC), padType and fitType
 859    // into the fitParam TVectorD (which should contain 8 elements).
 860    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 861    // Note: The fitter has to be evaluated first!
 862    //
 863   TLinearFitter * fitter = GetFitter(segment, padType, fitType);
 864   if (fitter){
 865     fitter->Eval();
 866     fitter->GetParameters(fitParam);
 867     return kTRUE;
 868   }else{
 869     Error("AliTPCcalibTracksGain::GetParameters","Fitter%d_%d_%d not available", segment, padType, fitType);
 870     return kFALSE;
 871   }
 872   return kFALSE;
 873 }
 874
 875 void AliTPCcalibTracksGain::GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError) {
 876    //
 877    // Puts the fit parameter errors for the specified segment (IROC & OROC), padType and fitType
 878    // into the fitError TVectorD (which should contain 8 elements).
 879    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 880    // Note: The fitter has to be evaluated first!
 881    //
 882
 883    GetFitter(segment, padType, fitType)->GetErrors(fitError);
 884    //fitError *= TMath::Sqrt(GetRedChi2(segment, padType, fitType));
 885 }
 886
 887
 888 void AliTPCcalibTracksGain::GetCovarianceMatrix(UInt_t segment, UInt_t padType, UInt_t fitType, TMatrixD& covMatrix) {
 889    //
 890    // Returns the covariance matrix for the specified segment, padType, fitType.
 891    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 892    //
 893
 894    GetFitter(segment, padType, fitType)->GetCovarianceMatrix(covMatrix);
 895 }
 896
 897 TLinearFitter* AliTPCcalibTracksGain::GetFitter(UInt_t segment, UInt_t padType, UInt_t fitType) {
 898    //
 899    // Returns the TLinearFitter object for the specified segment, padType, fitType.
 900    // padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
 901    //
 902
 903    switch (fitType) {
 904       case kSimpleFitter:
 905          return fSimpleFitter->GetFitter(segment, padType);
 906       case kSqrtFitter:
 907          return fSqrtFitter->GetFitter(segment, padType);
 908       case kLogFitter:
 909          return fLogFitter->GetFitter(segment, padType);
 910       case 3:
 911          return fSingleSectorFitter->GetFitter(0, padType);
 912    }
 913    return 0;
 914 }
 915
 916 Double_t AliTPCcalibTracksGain::GetPadLength(Double_t lx) {
 917    //
 918    // The function returns 0.75 for an IROC, 1. for an OROC at medium pad size position,
 919    // 1.5 for an OROC at long pad size position, -1 if out of bounds.
 920    //
 921
 922    Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
 923    Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
 924    Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
 925    Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
 926    Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
 927    Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
 928
 929    // if IROC
 930    if (lx >= irocLow && lx <= irocUp) return 0.75;
 931    // if OROC medium pads
 932    if (lx >= orocLow1 && lx <= orocUp1) return 1.;
 933    // if OROC long pads
 934    if (lx >= orocLow2 && lx <= orocUp2) return 1.5;
 935    // if out of bounds
 936    return -1;
 937 }
 938
 939 Int_t AliTPCcalibTracksGain::GetPadType(Double_t lx) {
 940    //
 941    // The function returns 0 for an IROC, 1 for an OROC at medium pad size position,
 942    // 2 for an OROC at long pad size position, -1 if out of bounds.
 943    //
 944
 945    if (GetPadLength(lx) == 0.75) return 0;
 946    else if (GetPadLength(lx) == 1.) return 1;
 947    else if (GetPadLength(lx) == 1.5) return 2;
 948    return -1;
 949 }
 950
 951 void AliTPCcalibTracksGain::DumpTrack(AliTPCseed* track) {
 952    //
 953    //  Dump track information to the debug stream
 954    //
 955
 956    Int_t rows[200];
 957    Int_t sector[3];
 958    Int_t npoints[3];
 959    TVectorD dedxM[3];
 960    TVectorD dedxQ[3];
 961    TVectorD parY[3];
 962    TVectorD parZ[3];
 963    TVectorD meanPos[3];
 964    //
 965    Int_t count=0;
 966    for (Int_t ipad = 0; ipad < 3; ipad++) {
 967      dedxM[ipad].ResizeTo(5);
 968      dedxQ[ipad].ResizeTo(5);
 969      parY[ipad].ResizeTo(3);
 970      parZ[ipad].ResizeTo(3);
 971      meanPos[ipad].ResizeTo(6);
 972      Bool_t isOK = GetDedx(track, ipad, rows, sector[ipad], npoints[ipad], dedxM[ipad], dedxQ[ipad], parY[ipad], parZ[ipad], meanPos[ipad]);
 973      if (isOK)
 974        AddTracklet(sector[ipad],ipad, dedxQ[ipad], dedxM[ipad], parY[ipad], parZ[ipad], meanPos[ipad] );
 975      if (isOK) count++;
 976    }
 977
 978    TTreeSRedirector * cstream =  GetDebugStreamer();
 979    if (cstream){
 980      (*cstream) << "Track" <<
 981        "run="<<fRun<<              //  run number
 982        "event="<<fEvent<<          //  event number
 983        "time="<<fTime<<            //  time stamp of event
 984        "trigger="<<fTrigger<<      //  trigger
 985        "mag="<<fMagF<<             //  magnetic field
 986        "Track.=" << track <<        // track information
 987        "\n";
 988      //
 989      //
 990      //
 991      if ( GetStreamLevel()>1 && count>1){
 992        (*cstream) << "TrackG" <<
 993          "run="<<fRun<<              //  run number
 994          "event="<<fEvent<<          //  event number
 995          "time="<<fTime<<            //  time stamp of event
 996          "trigger="<<fTrigger<<      //  trigger
 997          "mag="<<fMagF<<             //  magnetic field
 998          "Track.=" << track <<        // track information
 999          "ncount="<<count<<
1000          // info for pad type 0
1001          "sector0="<<sector[0]<<
1002          "npoints0="<<npoints[0]<<
1003          "dedxM0.="<<&dedxM[0]<<
1004          "dedxQ0.="<<&dedxQ[0]<<
1005          "parY0.="<<&parY[0]<<
1006          "parZ0.="<<&parZ[0]<<
1007          "meanPos0.="<<&meanPos[0]<<
1008          //
1009          // info for pad type 1
1010          "sector1="<<sector[1]<<
1011          "npoints1="<<npoints[1]<<
1012          "dedxM1.="<<&dedxM[1]<<
1013          "dedxQ1.="<<&dedxQ[1]<<
1014          "parY1.="<<&parY[1]<<
1015          "parZ1.="<<&parZ[1]<<
1016          "meanPos1.="<<&meanPos[1]<<
1017          //
1018          // info for pad type 2
1019          "sector2="<<sector[2]<<
1020          "npoints2="<<npoints[2]<<
1021          "dedxM2.="<<&dedxM[2]<<
1022          "dedxQ2.="<<&dedxQ[2]<<
1023          "parY2.="<<&parY[2]<<
1024          "parZ2.="<<&parZ[2]<<
1025          "meanPos2.="<<&meanPos[2]<<
1026          //
1027          "\n";
1028      }
1029    }
1030 }
1031
1032 Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /*rows*/,
1033                                       Int_t &sector, Int_t& npoints,
1034                                       TVectorD &dedxM, TVectorD &dedxQ,
1035                                       TVectorD &parY, TVectorD &parZ, TVectorD&meanPos)
1036 {
1037   //
1038   // GetDedx for given sector for given track
1039   // padType - type of pads
1040   //
1041
1042   static TLinearFitter fitY(2, "pol1");
1043   static TLinearFitter fitZ(2, "pol1");
1044   fitY.StoreData(kFALSE);
1045   fitZ.StoreData(kFALSE);
1046   //
1047   //
1048   Int_t firstRow = 0, lastRow = 0;
1049   Int_t minRow = 100;
1050   Float_t xcenter = 0;
1051   const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
1052   const Float_t kedgey = 4.;
1053   if (padType == 0) {
1054     firstRow = 0;
1055     lastRow = fgTPCparam->GetNRowLow();
1056     xcenter = 108.47;
1057   }
1058   if (padType == 1) {
1059     firstRow = fgTPCparam->GetNRowLow();
1060     lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1061     xcenter = 166.60;
1062    }
1063   if (padType == 2) {
1064     firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
1065     lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
1066     xcenter = 222.6;
1067   }
1068   minRow = (lastRow - firstRow) / 2;
1069   //
1070   //
1071   Int_t nclusters = 0;
1072   Int_t nclustersNE = 0; // number of not edge clusters
1073   Int_t lastSector = -1;
1074   Float_t amplitudeQ[100];
1075   Float_t amplitudeM[100];
1076    Int_t rowIn[100];
1077    Int_t index[100];
1078    //
1079    //
1080    fitY.ClearPoints();
1081    fitZ.ClearPoints();
1082
1083    for (Int_t iCluster = firstRow; iCluster < lastRow; iCluster++) {
1084       AliTPCclusterMI* cluster = track->GetClusterPointer(iCluster);
1085       if (cluster) {
1086          Int_t detector = cluster->GetDetector() ;
1087          if (lastSector == -1) lastSector = detector;
1088          if (lastSector != detector) continue;
1089          amplitudeQ[nclusters] = cluster->GetQ()/GetQNorm(cluster,parY[1], parZ[1]);
1090          amplitudeM[nclusters] = cluster->GetMax()/GetMaxNorm(cluster,parY[1], parZ[1]);
1091          rowIn[nclusters] = iCluster;
1092          nclusters++;
1093          Double_t dx = cluster->GetX() - xcenter;
1094          Double_t y = cluster->GetY();
1095          Double_t z = cluster->GetZ();
1096          fitY.AddPoint(&dx, y);
1097          fitZ.AddPoint(&dx, z);
1098          meanPos[0] += dx;
1099          meanPos[1] += dx;
1100          meanPos[2] += y;
1101          meanPos[3] += y*y;
1102          meanPos[4] += z;
1103          meanPos[5] += z*z;
1104          if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) nclustersNE++;
1105       }
1106    }
1107
1108    if (nclusters < minRow / 2) return kFALSE;
1109    if (nclustersNE < minRow / 2) return kFALSE;
1110    for (Int_t i = 0; i < 6; i++) meanPos[i] /= Double_t(nclusters);
1111    fitY.Eval();
1112    fitZ.Eval();
1113    fitY.GetParameters(parY);
1114    fitZ.GetParameters(parZ);
1115    //
1116    // calculate truncated mean
1117    //
1118    TMath::Sort(nclusters, amplitudeQ, index, kFALSE);
1119    //
1120    //
1121    //
1122    Float_t ndedx[5];
1123    for (Int_t i = 0; i < 5; i++) {
1124       dedxQ[i] = 0;
1125       dedxM[i] = 0;
1126       ndedx[i] = 0;
1127    }
1128    //
1129    // dedx calculation
1130    //
1131    Int_t inonEdge = 0;
1132    for (Int_t i = 0; i < nclusters; i++) {
1133       Int_t rowSorted = rowIn[index[i]];
1134       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1135
1136       if (TMath::Abs(cluster->GetY()) > cluster->GetX()*ktany - kedgey) continue;  //don't take edge clusters
1137       inonEdge++;
1138       if (inonEdge < nclustersNE * 0.5) {
1139          ndedx[0]++;
1140          dedxQ[0] += amplitudeQ[index[i]];
1141          dedxM[0] += amplitudeM[index[i]];
1142       }
1143       if (inonEdge < nclustersNE * 0.6) {
1144          ndedx[1]++;
1145          dedxQ[1] += amplitudeQ[index[i]];
1146          dedxM[1] += amplitudeM[index[i]];
1147       }
1148       if (inonEdge < nclustersNE * 0.7) {
1149          ndedx[2]++;
1150          dedxQ[2] += amplitudeQ[index[i]];
1151          dedxM[2] += amplitudeM[index[i]];
1152       }
1153       if (inonEdge < nclustersNE * 0.8) {
1154          ndedx[3]++;
1155          dedxQ[3] += amplitudeQ[index[i]];
1156          dedxM[3] += amplitudeM[index[i]];
1157       }
1158       if (inonEdge < nclustersNE * 0.9) {
1159          ndedx[4]++;
1160          dedxQ[4] += amplitudeQ[index[i]];
1161          dedxM[4] += amplitudeM[index[i]];
1162       }
1163    }
1164    for (Int_t i = 0; i < 5; i++) {
1165       dedxQ[i] /= ndedx[i];
1166       dedxM[i] /= ndedx[i];
1167    }
1168    TTreeSRedirector * cstream =  GetDebugStreamer();
1169    inonEdge = 0;
1170    Float_t momenta = track->GetP();
1171    Float_t mdedx = track->GetdEdx();
1172    for (Int_t i = 0; i < nclusters; i++) {
1173       Int_t rowSorted = rowIn[index[i]];
1174       AliTPCclusterMI* cluster = track->GetClusterPointer(rowSorted);
1175       if (!cluster) {
1176          printf("Problem\n");
1177          continue;
1178       }
1179       if (TMath::Abs(cluster->GetY()) < cluster->GetX()*ktany - kedgey) inonEdge++;
1180       Float_t dedge = cluster->GetX()*ktany - TMath::Abs(cluster->GetY());
1181       Float_t fraction = Float_t(i) / Float_t(nclusters);
1182       Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
1183
1184       AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
1185
1186       Double_t qNorm = GetQNorm(cluster,parY[1], parZ[1]);
1187       Double_t mNorm = GetMaxNorm(cluster,parY[1], parZ[1]);
1188
1189
1190       Float_t gain = GetGain(cluster);
1191       if (cstream) (*cstream) << "dEdx" <<
1192         "run="<<fRun<<              //  run number
1193         "event="<<fEvent<<          //  event number
1194         "time="<<fTime<<            //  time stamp of event
1195         "trigger="<<fTrigger<<      //  trigger
1196         "mag="<<fMagF<<             //  magnetic field
1197
1198         "Cl.=" << cluster <<           // cluster of interest
1199         "gain="<<gain<<                // gain at cluster position
1200         "mNorm="<<mNorm<<              // Q max normalization
1201         "qNorm="<<qNorm<<              // Q tot normalization
1202         "P=" << momenta <<             // track momenta
1203         "dedx=" << mdedx <<            // mean dedx - corrected for angle
1204         "IPad=" << padType <<          // pad type 0..2
1205         "xc=" << xcenter <<            // x center of chamber
1206         "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
1207         "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
1208         "fraction=" << fraction <<     // fraction - order in statistic (0,1)
1209         "fraction2=" << fraction2 <<   // fraction - order in statistic (0,1)
1210         "dedge=" << dedge <<           // distance to the edge
1211         "parY.=" << &parY <<           // line fit
1212         "parZ.=" << &parZ <<           // line fit
1213         "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
1214         "\n";
1215    }
1216
1217    if (cstream) (*cstream) << "dEdxT" <<
1218      "run="<<fRun<<              //  run number
1219      "event="<<fEvent<<          //  event number
1220      "time="<<fTime<<            //  time stamp of event
1221      "trigger="<<fTrigger<<      //  trigger
1222      "mag="<<fMagF<<             //  magnetic field
1223      "P=" << momenta <<             // track momenta
1224      "npoints="<<inonEdge<<         // number of points
1225      "sector="<<lastSector<<        // sector number
1226      "dedx=" << mdedx <<            // mean dedx - corrected for angle
1227      "IPad=" << padType <<          // pad type 0..2
1228      "xc=" << xcenter <<            // x center of chamber
1229      "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
1230      "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
1231      "parY.=" << &parY <<           // line fit
1232      "parZ.=" << &parZ <<           // line fit
1233      "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
1234      "\n";
1235
1236    sector = lastSector;
1237    npoints = inonEdge;
1238    return kTRUE;
1239 }
1240
1241 void AliTPCcalibTracksGain::AddTracklet(UInt_t sector, UInt_t padType,TVectorD &dedxQ, TVectorD &dedxM,TVectorD& parY, TVectorD& parZ, TVectorD& meanPos){
1242   //
1243   // Add measured point - dedx to the fitter
1244   //
1245   //
1246   //chain->SetAlias("dr","(250-abs(meanPos.fElements[4]))/250");
1247   //chain->SetAlias("tz","(0+abs(parZ.fElements[1]))");
1248   //chain->SetAlias("ty","(0+abs(parY.fElements[1]))");
1249   //chain->SetAlias("corrg","sqrt((1+ty^2)*(1+tz^2))");
1250   //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);
1251
1252   Double_t xxx[100];
1253   //
1254   // z and angular part
1255   //
1256
1257   xxx[0] = (250.-TMath::Abs(meanPos[4]))/250.;
1258   xxx[1] = TMath::Abs(parY[1]);
1259   xxx[2] = TMath::Abs(parZ[1]);
1260   xxx[3] = xxx[0]*xxx[1];
1261   xxx[4] = xxx[0]*xxx[2];
1262   xxx[5] = xxx[1]*xxx[2];
1263   xxx[6] = xxx[0]*xxx[0];
1264   xxx[7] = xxx[1]*xxx[1];
1265   xxx[8] = xxx[2]*xxx[2];
1266   //
1267   // chamber part
1268   //
1269   Int_t tsector = sector%36;
1270   for (Int_t i=0;i<35;i++){
1271     xxx[9+i]=(i==tsector)?1:0;
1272   }
1273   TLinearFitter *fitterM = fFitter0M;
1274   if (padType==1) fitterM=fFitter1M;
1275   if (padType==2) fitterM=fFitter2M;
1276   fitterM->AddPoint(xxx,dedxM[1]);
1277   //
1278   TLinearFitter *fitterT = fFitter0T;
1279   if (padType==1) fitterT = fFitter1T;
1280   if (padType==2) fitterT = fFitter2T;
1281   fitterT->AddPoint(xxx,dedxQ[1]);
1282   //
1283   TLinearFitter *dfitterM = fDFitter0M;
1284   if (padType==1) dfitterM=fDFitter1M;
1285   if (padType==2) dfitterM=fDFitter2M;
1286   dfitterM->AddPoint(xxx,dedxM[1]);
1287   //
1288   TLinearFitter *dfitterT = fDFitter0T;
1289   if (padType==1) dfitterT = fDFitter1T;
1290   if (padType==2) dfitterT = fDFitter2T;
1291   dfitterT->AddPoint(xxx,dedxQ[1]);
1292 }
1293
1294
1295 TGraph *AliTPCcalibTracksGain::CreateAmpGraph(Int_t ipad, Bool_t qmax){
1296   //
1297   // create the amplitude graph
1298   // The normalized amplitudes are extrapolated to the 0 angle (y,z)  and 0 drift length
1299   //
1300
1301   TVectorD vec;
1302   if (qmax){
1303     if (ipad==0) fFitter0M->GetParameters(vec);
1304     if (ipad==1) fFitter1M->GetParameters(vec);
1305     if (ipad==2) fFitter2M->GetParameters(vec);
1306   }else{
1307     if (ipad==0) fFitter0T->GetParameters(vec);
1308     if (ipad==1) fFitter1T->GetParameters(vec);
1309     if (ipad==2) fFitter2T->GetParameters(vec);
1310   }
1311
1312   Float_t amp[36];
1313   Float_t sec[36];
1314   for (Int_t i=0;i<35;i++){
1315     sec[i]=i;
1316     amp[i]=vec[10+i]+vec[0];
1317   }
1318   amp[35]=vec[0];
1319   Float_t mean = TMath::Mean(36,amp);
1320   for (Int_t i=0;i<36;i++){
1321     sec[i]=i;
1322     amp[i]=(amp[i]-mean)/mean;
1323   }
1324   TGraph *gr = new TGraph(36,sec,amp);
1325   return gr;
1326 }
1327
1328
1329 void   AliTPCcalibTracksGain::UpdateClusterParam(AliTPCClusterParam* clparam){
1330   //
1331   //   SetQ normalization parameters
1332   //
1333   //  void SetQnorm(Int_t ipad, Int_t itype,  TVectorD * norm);
1334
1335   TVectorD vec;
1336
1337   //
1338   fDFitter0T->Eval();
1339   fDFitter1T->Eval();
1340   fDFitter2T->Eval();
1341   fDFitter0M->Eval();
1342   fDFitter1M->Eval();
1343   fDFitter2M->Eval();
1344   fDFitter0T->GetParameters(vec);
1345   clparam->SetQnorm(0,0,&vec);
1346   fDFitter1T->GetParameters(vec);
1347   clparam->SetQnorm(1,0,&vec);
1348   fDFitter2T->GetParameters(vec);
1349   clparam->SetQnorm(2,0,&vec);
1350   //
1351   fDFitter0M->GetParameters(vec);
1352   clparam->SetQnorm(0,1,&vec);
1353   fDFitter1M->GetParameters(vec);
1354   clparam->SetQnorm(1,1,&vec);
1355   fDFitter2M->GetParameters(vec);
1356   clparam->SetQnorm(2,1,&vec);
1357   //
1358
1359 }
1360
1361
1362 void   AliTPCcalibTracksGain::Analyze(){
1363
1364  Evaluate();
1365
1366 }
1367
1368