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