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