TPC/AliTPCcalibTracks.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 //     Class for calibration of the cluster properties:                       //
  20 //     Cluster position resolution (RMS)  and short term distortions (within pad or within time bin)
  21 //
  22 //  Algorithm:
  23 //    1. Creation of the residual/properties  histograms
  24 //    2. Filling of the histograms.
  25 //       2.a Traklet fitting
  26 //       2.b Resudual filling
  27 //    3. Postprocessing: Creation of the tree with the mean values and RMS at different bins
  28 //    4.               : Paramaterization fitting.
  29 //                       In old AliRoot the RMS paramterization was used to characterize cluster resolution
  30 //                       Mean value /bias was neglected
  31 //    5. To be implemented
  32 //          a.) lookup table for the distortion parmaterization - THn
  33 //          b.) Usage in the reconstruction
  34 //
  35 //
  36 //   1. Creation of the histograms:   MakeHistos()
  37 //
  38 //         6 n dimensional histograms  are filled during the calibration:
  39 //         cluster  performance bins
  40 //         0 - variable of interest
  41 //         1 - pad type   - 0- IROC Short 1- OCROC medium 2 - OROC long pads
  42 //         2 - drift length - drift length -0-1
  43 //         3 - Qmax         - Qmax  - 2- 400
  44 //         4 - cog          - COG position - 0-1
  45 //         5 - tan(phi)     - local  angle
  46 //        Histograms:
  47 //        THnSparse  *fHisDeltaY;    // THnSparse - delta Y between the cluster and tracklet interpolation (+-X(5?) rows)
  48 //        THnSparse  *fHisDeltaZ;    // THnSparse - delta Z
  49 //        THnSparse  *fHisRMSY;      // THnSparse - rms Y
  50 //        THnSparse  *fHisRMSZ;      // THnSparse - rms Z
  51 //        THnSparse  *fHisQmax;      // THnSparse - qmax
  52 //        THnSparse  *fHisQtot;      // THnSparse - qtot
  53 //
  54 //
  55 //
  56
  57
  58                      //
  59 //     The parameter 'clusterParam', a AliTPCClusterParam object             //
  60 //      (needed for TPC cluster error and shape parameterization)            //
  61 //
  62 //     Normally you get this object out of the file 'TPCClusterParam.root'   //
  63 //     In the parameter 'cuts' the cuts are specified, that decide           //
  64 //     weather a track will be accepted for calibration or not.              //
  65 //                                                                           //
  66 //
  67 //     The data flow:
  68 //
  69 /*
  70    Raw Data -> Local Reconstruction -> Tracking ->     Calibration -> RefData (component itself)
  71                Offline/HLT             Offline/HLT                    OCDB entries (AliTPCClusterParam)
  72 */
  73
  74 /*
  75   Example usage - creation of the residual trees:
  76   TFile f("CalibObjects.root");
  77   AliTPCcalibTracks *calibTracks = ( AliTPCcalibTracks *)TPCCalib->FindObject("calibTracks");
  78   TH2 * his2 = calibTracks->fHisDeltaY->Projection(0,4);
  79   his2->SetName("his2");
  80   his2->FitSlicesY();
  81
  82
  83   TTreeSRedirector *pcstream = new TTreeSRedirector("clusterLookup.root");
  84   AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaY, pcstream, 0);
  85   AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaZ, pcstream, 1);
  86   delete pcstream;
  87   TFile fl("clusterLookup.root");
  88   TCut cutNI="cogA==0&&angleA==0&&qmaxA==0&&zA==0&&ipadA==0"
  89
  90 */
  91
  92                                                                //
  93 ///////////////////////////////////////////////////////////////////////////////
  94
  95 //
  96 // ROOT includes
  97 //
  98 #include <iostream>
  99 #include <fstream>
 100 using namespace std;
 101
 102 #include <TROOT.h>
 103 #include <TChain.h>
 104 #include <TFile.h>
 105 #include <TH3F.h>
 106 #include <TProfile.h>
 107
 108 //
 109 //#include <TPDGCode.h>
 110 #include <TStyle.h>
 111 #include "TLinearFitter.h"
 112 //#include "TMatrixD.h"
 113 #include "TTreeStream.h"
 114 #include "TF1.h"
 115 #include <TCanvas.h>
 116 #include <TGraph2DErrors.h>
 117 #include "TPostScript.h"
 118 #include "TCint.h"
 119
 120 #include <TH2D.h>
 121 #include <TF2.h>
 122 #include <TSystem.h>
 123 #include <TCollection.h>
 124 #include <iostream>
 125 #include <TLinearFitter.h>
 126 #include <TString.h>
 127
 128 //
 129 // AliROOT includes
 130 //
 131 #include "AliMagF.h"
 132 #include "AliTracker.h"
 133 #include "AliESD.h"
 134 #include "AliESDtrack.h"
 135 #include "AliESDfriend.h"
 136 #include "AliESDfriendTrack.h"
 137 #include "AliTPCseed.h"
 138 #include "AliTPCclusterMI.h"
 139 #include "AliTPCROC.h"
 140
 141 #include "AliTPCParamSR.h"
 142 #include "AliTrackPointArray.h"
 143 #include "AliTPCcalibTracks.h"
 144 #include "AliTPCClusterParam.h"
 145 #include "AliTPCcalibTracksCuts.h"
 146 #include "AliTPCCalPad.h"
 147 #include "AliTPCCalROC.h"
 148 #include "TText.h"
 149 #include "TPaveText.h"
 150 #include "TSystem.h"
 151 #include "TStatToolkit.h"
 152 #include "TCut.h"
 153 #include "THnSparse.h"
 154 #include "AliRieman.h"
 155
 156
 157
 158 ClassImp(AliTPCcalibTracks)
 159
 160
 161 AliTPCcalibTracks::AliTPCcalibTracks():
 162   AliTPCcalibBase(),
 163   fClusterParam(0),
 164   fROC(0),
 165   fHisDeltaY(0),    // THnSparse - delta Y
 166   fHisDeltaZ(0),    // THnSparse - delta Z
 167   fHisRMSY(0),      // THnSparse - rms Y
 168   fHisRMSZ(0),      // THnSparse - rms Z
 169   fHisQmax(0),      // THnSparse - qmax
 170   fHisQtot(0),      // THnSparse - qtot
 171   fArrayQDY(0),
 172   fArrayQDZ(0),
 173   fArrayQRMSY(0),
 174   fArrayQRMSZ(0),
 175   fResolY(0),
 176   fResolZ(0),
 177   fRMSY(0),
 178   fRMSZ(0),
 179   fCuts(0),
 180   fRejectedTracksHisto(0),
 181   fClusterCutHisto(0),
 182   fCalPadClusterPerPad(0),
 183   fCalPadClusterPerPadRaw(0)
 184 {
 185    //
 186    // AliTPCcalibTracks default constructor
 187    //
 188   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;
 189 }
 190
 191
 192
 193 AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
 194   AliTPCcalibBase(calibTracks),
 195   fClusterParam(0),
 196   fROC(0),
 197   fHisDeltaY(0),    // THnSparse - delta Y
 198   fHisDeltaZ(0),    // THnSparse - delta Z
 199   fHisRMSY(0),      // THnSparse - rms Y
 200   fHisRMSZ(0),      // THnSparse - rms Z
 201   fHisQmax(0),      // THnSparse - qmax
 202   fHisQtot(0),      // THnSparse - qtot
 203   fArrayQDY(0),
 204   fArrayQDZ(0),
 205   fArrayQRMSY(0),
 206   fArrayQRMSZ(0),
 207   fResolY(0),
 208   fResolZ(0),
 209   fRMSY(0),
 210   fRMSZ(0),
 211   fCuts(0),
 212   fRejectedTracksHisto(0),
 213   fClusterCutHisto(0),
 214   fCalPadClusterPerPad(0),
 215   fCalPadClusterPerPadRaw(0)
 216 {
 217    //
 218    // AliTPCcalibTracks copy constructor
 219    //
 220   if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
 221
 222    Bool_t dirStatus = TH1::AddDirectoryStatus();
 223    TH1::AddDirectory(kFALSE);
 224
 225    Int_t length = -1;
 226
 227    (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
 228    fArrayQDY= new TObjArray(length);
 229    fArrayQDZ= new TObjArray(length);
 230    fArrayQRMSY= new TObjArray(length);
 231    fArrayQRMSZ= new TObjArray(length);
 232    for (Int_t i = 0; i < length; i++) {
 233       fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
 234       fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
 235       fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
 236       fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
 237    }
 238
 239    (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
 240    fResolY = new TObjArray(length);
 241    fResolZ = new TObjArray(length);
 242    fRMSY = new TObjArray(length);
 243    fRMSZ = new TObjArray(length);
 244    for (Int_t i = 0; i < length; i++) {
 245       fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
 246       fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
 247       fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
 248       fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
 249    }
 250
 251
 252    fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
 253    fRejectedTracksHisto    = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
 254    fCalPadClusterPerPad    = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
 255    fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
 256
 257    fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(),
 258       calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
 259    SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
 260    TH1::AddDirectory(dirStatus); // set status back to original status
 261 //    cout << "+++++ end of copy constructor +++++" << endl;   // TO BE REMOVED
 262 }
 263
 264
 265 AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
 266   //
 267   // assgnment operator
 268   //
 269   if (this != &calibTracks) {
 270     new (this) AliTPCcalibTracks(calibTracks);
 271   }
 272   return *this;
 273
 274 }
 275
 276
 277 AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam,  AliTPCcalibTracksCuts* cuts, Int_t logLevel) :
 278   AliTPCcalibBase(),
 279   fClusterParam(0),
 280   fROC(0),
 281   fHisDeltaY(0),    // THnSparse - delta Y
 282   fHisDeltaZ(0),    // THnSparse - delta Z
 283   fHisRMSY(0),      // THnSparse - rms Y
 284   fHisRMSZ(0),      // THnSparse - rms Z
 285   fHisQmax(0),      // THnSparse - qmax
 286   fHisQtot(0),      // THnSparse - qtot
 287   fArrayQDY(0),
 288   fArrayQDZ(0),
 289   fArrayQRMSY(0),
 290   fArrayQRMSZ(0),
 291   fResolY(0),
 292   fResolZ(0),
 293   fRMSY(0),
 294   fRMSZ(0),
 295   fCuts(0),
 296   fRejectedTracksHisto(0),
 297   fClusterCutHisto(0),
 298   fCalPadClusterPerPad(0),
 299   fCalPadClusterPerPadRaw(0)
 300  {
 301    //
 302    // AliTPCcalibTracks constructor
 303    // specify 'name' and 'title' of your object
 304    // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
 305    // In the parameter 'cuts' the cuts are specified, that decide
 306    // weather a track will be accepted for calibration or not.
 307    //
 308    // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
 309    //
 310    // All histograms are instatiated in this constructor.
 311    //
 312    this->SetName(name);
 313    this->SetTitle(title);
 314
 315    if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
 316    G__SetCatchException(0);
 317
 318    fClusterParam = clusterParam;
 319    if (fClusterParam){
 320      fClusterParam->SetInstance(fClusterParam);
 321    }
 322    else {
 323      Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
 324    }
 325    fCuts = cuts;
 326    SetDebugLevel(logLevel);
 327    MakeHistos();
 328
 329    TH1::AddDirectory(kFALSE);
 330
 331    fRejectedTracksHisto    = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
 332    fCalPadClusterPerPad    = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
 333    fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
 334    fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
 335
 336
 337    TH1::AddDirectory(kFALSE);
 338
 339
 340    fResolY = new TObjArray(3);
 341    fResolZ = new TObjArray(3);
 342    fRMSY   = new TObjArray(3);
 343    fRMSZ   = new TObjArray(3);
 344    TH3F * his3D;
 345    //
 346    his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
 347    fResolY->AddAt(his3D,0);
 348    his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
 349    fResolY->AddAt(his3D,1);
 350    his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
 351    fResolY->AddAt(his3D,2);
 352    //
 353    his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
 354    fResolZ->AddAt(his3D,0);
 355    his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
 356    fResolZ->AddAt(his3D,1);
 357    his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
 358    fResolZ->AddAt(his3D,2);
 359    //
 360    his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
 361    fRMSY->AddAt(his3D,0);
 362    his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
 363    fRMSY->AddAt(his3D,1);
 364    his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
 365    fRMSY->AddAt(his3D,2);
 366    //
 367    his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
 368    fRMSZ->AddAt(his3D,0);
 369    his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
 370    fRMSZ->AddAt(his3D,1);
 371    his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
 372    fRMSZ->AddAt(his3D,2);
 373    //
 374
 375    TH1::AddDirectory(kFALSE);
 376
 377    fArrayQDY = new TObjArray(300);
 378    fArrayQDZ = new TObjArray(300);
 379    fArrayQRMSY = new TObjArray(300);
 380    fArrayQRMSZ = new TObjArray(300);
 381    for (Int_t iq = 0; iq <= 10; iq++){
 382       for (Int_t ipad = 0; ipad < 3; ipad++){
 383          Int_t   bin   = GetBin(iq, ipad);
 384          Float_t qmean = GetQ(bin);
 385          char hname[200];
 386          snprintf(hname,100,"ResolY Pad%d Qmiddle%f",ipad, qmean);
 387          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
 388          fArrayQDY->AddAt(his3D, bin);
 389          snprintf(hname,100,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
 390          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
 391          fArrayQDZ->AddAt(his3D, bin);
 392          snprintf(hname,100,"RMSY Pad%d Qmiddle%f",ipad, qmean);
 393          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
 394          fArrayQRMSY->AddAt(his3D, bin);
 395          snprintf(hname,100,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
 396          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
 397          fArrayQRMSZ->AddAt(his3D, bin);
 398       }
 399    }
 400
 401
 402
 403    if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl;
 404    cout << "end of main constructor" << endl; // TO BE REMOVED
 405 }
 406
 407
 408 AliTPCcalibTracks::~AliTPCcalibTracks() {
 409    //
 410    // AliTPCcalibTracks destructor
 411    //
 412
 413   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
 414    Int_t length = 0;
 415
 416
 417    if (fResolY) {
 418      length = fResolY->GetEntriesFast();
 419      for (Int_t i = 0; i < length; i++){
 420        delete fResolY->At(i);
 421        delete fResolZ->At(i);
 422        delete fRMSY->At(i);
 423        delete fRMSZ->At(i);
 424      }
 425      delete fResolY;
 426      delete fResolZ;
 427      delete fRMSY;
 428      delete fRMSZ;
 429    }
 430
 431    if (fArrayQDY) {
 432      length = fArrayQDY->GetEntriesFast();
 433      for (Int_t i = 0; i < length; i++){
 434        delete fArrayQDY->At(i);
 435        delete fArrayQDZ->At(i);
 436        delete fArrayQRMSY->At(i);
 437        delete fArrayQRMSZ->At(i);
 438      }
 439    }
 440
 441
 442
 443    delete fArrayQDY;
 444    delete fArrayQDZ;
 445    delete fArrayQRMSY;
 446    delete fArrayQRMSZ;
 447
 448   delete fRejectedTracksHisto;
 449   delete fClusterCutHisto;
 450   if (fCalPadClusterPerPad)    delete fCalPadClusterPerPad;
 451   if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
 452   delete fHisDeltaY;    // THnSparse - delta Y
 453   delete fHisDeltaZ;    // THnSparse - delta Z
 454   delete fHisRMSY;      // THnSparse - rms Y
 455   delete fHisRMSZ;      // THnSparse - rms Z
 456   delete fHisQmax;      // THnSparse - qmax
 457   delete fHisQtot;      // THnSparse - qtot
 458
 459 }
 460
 461
 462
 463 void AliTPCcalibTracks::Process(AliTPCseed *track){
 464    //
 465    // To be called in the selector
 466    // first AcceptTrack is evaluated, then calls all the following analyse functions:
 467    // FillResolutionHistoLocal(track)
 468
 469    //
 470   if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
 471    Int_t accpetStatus = AcceptTrack(track);
 472    if (accpetStatus == 0) {
 473       FillResolutionHistoLocal(track);
 474    }
 475    else fRejectedTracksHisto->Fill(accpetStatus);
 476 }
 477
 478
 479
 480 Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
 481   //
 482   // calculate bins for given q and pad type
 483   // used in TObjArray
 484   //
 485   Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );
 486   res *= 3;
 487   res += pad;
 488   return res;
 489 }
 490
 491
 492 Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
 493   //
 494   // calculate bins for given iq and pad type
 495   // used in TObjArray
 496   //
 497   return iq * 3 + pad;;
 498 }
 499
 500
 501 Float_t AliTPCcalibTracks::GetQ(Int_t bin){
 502    //
 503    // returns to bin belonging charge
 504    // (bin / 3 + 3)^2
 505    //
 506    Int_t bin0 = bin / 3;
 507    bin0 += 3;
 508    return bin0 * bin0;
 509 }
 510
 511
 512 Float_t AliTPCcalibTracks::GetPad(Int_t bin){
 513    //
 514    // returns to bin belonging pad
 515    // bin % 3
 516    //
 517    return bin % 3;
 518 }
 519
 520
 521
 522 Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
 523   //
 524   // Function, that decides wheather a given track is accepted for
 525   // the analysis or not.
 526   // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
 527   // Returns 0 if a track is accepted or an integer different from 0
 528   // to indicate the failed cut
 529   //
 530   const Int_t   kMinClusters  = fCuts->GetMinClusters();
 531   const Float_t kMinRatio     = fCuts->GetMinRatio();
 532   const Float_t kMax1pt       = fCuts->GetMax1pt();
 533   const Float_t kEdgeYXCutNoise    = fCuts->GetEdgeYXCutNoise();
 534   const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
 535
 536   //
 537   // edge induced noise tracks - NEXT RELEASE will be removed during tracking
 538   if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
 539     if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
 540   if (track->GetNumberOfClusters() < kMinClusters) return 2;
 541   Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
 542   if (ratio < kMinRatio) return 3;
 543 //   Float_t mpt = track->Get1Pt();       // Get1Pt() doesn't exist any more
 544   Float_t mpt = track->GetSigned1Pt();
 545   if (TMath::Abs(mpt) > kMax1pt) return 4;
 546   //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
 547   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
 548   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
 549
 550   if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");
 551   return 0;
 552 }
 553
 554
 555 void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
 556    //
 557    // fill resolution histograms - localy - tracklet in the neighborhood
 558    // write debug information to 'TPCSelectorDebug.root'
 559    //
 560    // _ the main function, called during track analysis _
 561    //
 562    // loop over all padrows along the track
 563    // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
 564    //
 565    // loop again over all padrows along the track
 566    // fit tracklet (clusters in given padrow +- kDelta padrows)
 567    // with polynom of 2nd order and two polynoms of 1st order
 568    // take both polynoms of 1st order, calculate difference of their parameters
 569    // add covariance matrixes and calculate chi2 of this difference
 570    // if this chi2 is bigger than a given threshold, assume that the current cluster is
 571    // a kink an goto next padrow
 572    // if not:
 573    // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
 574    //
 575    // write debug information to 'TPCSelectorDebug.root'
 576    // only for every kDeltaWriteDebugStream'th padrow to reduce data volume
 577    // and to avoid redundant data
 578    //
 579
 580   static TLinearFitter fFitterParY(3,"pol2");    //
 581   static TLinearFitter fFitterParZ(3,"pol2");    //
 582   static TLinearFitter fFitterParYWeight(3,"pol2");    //
 583   static TLinearFitter fFitterParZWeight(3,"pol2");    //
 584   fFitterParY.StoreData(kTRUE);
 585   fFitterParZ.StoreData(kTRUE);
 586   fFitterParYWeight.StoreData(kTRUE);
 587   fFitterParZWeight.StoreData(kTRUE);
 588   if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
 589   const Int_t   kDelta     = 10;          // delta rows to fit
 590   const Float_t kMinRatio  = 0.75;        // minimal ratio
 591   const Float_t kChi2Cut   =  10.;          // cut chi2 - left right
 592   const Float_t kSigmaCut  = 3.;        //sigma cut
 593   const Float_t kdEdxCut=300;
 594   const Float_t kPtCut=0.300;
 595
 596   if (track->GetTPCsignal()>kdEdxCut) return;
 597   if (TMath::Abs(AliTracker::GetBz())>0.1 &&TMath::Abs(track->Pt())<kPtCut) return;
 598
 599   // estimate mean error
 600   Int_t nTrackletsAll = 0;       // number of tracklets for given track
 601   Float_t csigmaY     = 0;       // mean sigma for tracklet refit in Y direction
 602   Float_t csigmaZ     = 0;       // mean sigma for tracklet refit in Z direction
 603   Int_t nClusters     = 0;       // working variable, number of clusters per tracklet
 604   Int_t sectorG       = -1;      // working variable, sector of tracklet, has to stay constant for one tracklet
 605   Double_t refX=0;
 606   // ---------------------------------------------------------------------
 607   for (Int_t irow = 0; irow < 159; irow++){
 608     // loop over all rows along the track
 609     // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
 610     // calculate mean chi^2 for this track-fit in Y and Z direction
 611     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 612     if (!cluster0) continue;  // no cluster found
 613     Int_t sector = cluster0->GetDetector();
 614
 615     Int_t ipad = TMath::Nint(cluster0->GetPad());
 616     Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
 617     fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
 618
 619     if (sector != sectorG){
 620       // track leaves sector before it crossed enough rows to fit / initialization
 621       nClusters = 0;
 622       fFitterParY.ClearPoints();
 623       fFitterParZ.ClearPoints();
 624       sectorG = sector;
 625     }
 626     else {
 627       nClusters++;
 628       if (refX<1) refX=cluster0->GetX()+kDelta*0.5;
 629       Double_t x = cluster0->GetX()-refX;
 630       fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
 631       fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
 632       //
 633       if ( nClusters >= kDelta + 3 ){
 634         // if more than 13 (kDelta+3) clusters were added to the fitters
 635         // fit the tracklet, increase trackletCounter
 636         fFitterParY.Eval();
 637         fFitterParZ.Eval();
 638         nTrackletsAll++;
 639         csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
 640         csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
 641         nClusters = -1;
 642         fFitterParY.ClearPoints();
 643         fFitterParZ.ClearPoints();
 644         refX=0;
 645       }
 646     }
 647   }      // for (Int_t irow = 0; irow < 159; irow++)
 648   // mean chi^2 for all tracklet fits in Y and in Z direction:
 649   csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
 650   csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
 651   // ---------------------------------------------------------------------
 652   //
 653   //
 654
 655   for (Int_t irow = kDelta; irow < 159-kDelta; irow++){
 656     // loop again over all rows along the track
 657     // do analysis
 658     //
 659     Int_t nclFound = 0;  // number of clusters in the neighborhood
 660     Int_t ncl0 = 0;      // number of clusters in rows < rowOfCenterCluster
 661     Int_t ncl1 = 0;      // number of clusters in rows > rowOfCenterCluster
 662     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
 663     if (!cluster0) continue;
 664     Int_t sector = cluster0->GetDetector();
 665     Float_t xref = cluster0->GetX();
 666
 667     // Make Fit
 668     fFitterParY.ClearPoints();
 669     fFitterParZ.ClearPoints();
 670     fFitterParYWeight.ClearPoints();
 671     fFitterParZWeight.ClearPoints();
 672     // fit tracklet (clusters in given padrow +- kDelta padrows)
 673     // with polynom of 2nd order and two polynoms of 1st order
 674     // take both polynoms of 1st order, calculate difference of their parameters
 675     // add covariance matrixes and calculate chi2 of this difference
 676     // if this chi2 is bigger than a given threshold, assume that the current cluster is
 677     // a kink an goto next padrow
 678     AliRieman riemanFit(2*kDelta);
 679     AliRieman riemanFitW(2*kDelta);
 680     for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
 681       // loop over irow +- kDelta rows (neighboured rows)
 682       //
 683       //
 684       if (idelta == 0) continue;                                // don't use center cluster
 685       if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
 686       AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
 687       if (!currentCluster) continue;
 688       if (currentCluster->GetType() < 0) continue;
 689       if (currentCluster->GetDetector() != sector) continue;
 690       nclFound++;
 691       if (idelta < 0){
 692         ncl0++;
 693       }
 694       if (idelta > 0){
 695         ncl1++;
 696       }
 697       riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ);
 698       riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY*TMath::Sqrt(1+TMath::Abs(idelta)),csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta)));
 699     }  // loop over neighbourhood for fitter filling
 700     if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
 701     if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
 702     if (nclFound < kDelta * kMinRatio) continue;    // if not enough clusters (7.5) found in neighbourhood goto next padrow
 703     riemanFit.Update();
 704     riemanFitW.Update();
 705     Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5));
 706     Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5));
 707     Double_t paramYR[3], paramZR[3];
 708     Double_t paramYRW[3], paramZRW[3];
 709     //
 710     paramYR[0]    = riemanFit.GetYat(xref);
 711     paramZR[0]    = riemanFit.GetZat(xref);
 712     paramYRW[0]   = riemanFitW.GetYat(xref);
 713     paramZRW[0]   = riemanFitW.GetZat(xref);
 714     //
 715     paramYR[1]    = riemanFit.GetDYat(xref);
 716     paramZR[1]    = riemanFit.GetDZat(xref);
 717     paramYRW[1]   = riemanFitW.GetDYat(xref);
 718     paramZRW[1]   = riemanFitW.GetDZat(xref);
 719     //
 720     Int_t reject=0;
 721     if (chi2R>kChi2Cut) reject+=1;
 722     if (chi2RW>kChi2Cut) reject+=2;
 723     if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4;
 724     if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8;
 725     if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16;
 726     if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32;
 727     //
 728     TTreeSRedirector *cstream = GetDebugStreamer();
 729     // get fit parameters from pol2 fit:
 730     Double_t tracky = paramYR[0];
 731     Double_t trackz = paramZR[0];
 732     Float_t  deltay = cluster0->GetY()-tracky;
 733     Float_t  deltaz = cluster0->GetZ()-trackz;
 734     Float_t  angley = paramYR[1];
 735     Float_t  anglez = paramZR[1];
 736
 737     Int_t padSize = 0;                          // short pads
 738     if (cluster0->GetDetector() >= 36) {
 739       padSize = 1;                              // medium pads
 740       if (cluster0->GetRow() > 63) padSize = 2; // long pads
 741     }
 742     Int_t ipad = TMath::Nint(cluster0->GetPad());
 743     if (cstream){
 744       Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
 745       (*cstream)<<"Resol0"<<
 746         "run="<<fRun<<              //  run number
 747         "event="<<fEvent<<          //  event number
 748         "time="<<fTime<<            //  time stamp of event
 749         "trigger="<<fTrigger<<      //  trigger
 750         "mag="<<fMagF<<             //  magnetic field
 751         "padSize="<<padSize<<
 752         //
 753         "reject="<<reject<<
 754         "cl.="<<cluster0<<          // cluster info
 755         "xref="<<xref<<             // cluster reference X
 756         //rieman fit
 757         "yr="<<paramYR[0]<<         // track position y - rieman fit
 758         "zr="<<paramZR[0]<<         // track position z - rieman fit
 759         "yrW="<<paramYRW[0]<<         // track position y - rieman fit - weight
 760         "zrW="<<paramZRW[0]<<         // track position z - rieman fit - weight
 761         "dyr="<<paramYR[1]<<         // track position y - rieman fit
 762         "dzr="<<paramZR[1]<<         // track position z - rieman fit
 763         "dyrW="<<paramYRW[1]<<         // track position y - rieman fit - weight
 764         "dzrW="<<paramZRW[1]<<         // track position z - rieman fit - weight
 765         //
 766         "angley="<<angley<<         // angle par fit
 767         "anglez="<<anglez<<         // angle par fit
 768         "zdr="<<zdrift<<            //
 769         "dy="<<deltay<<
 770         "dz="<<deltaz<<
 771         "sy="<<csigmaY<<
 772         "sz="<<csigmaZ<<
 773         "chi2R="<<chi2R<<
 774         "chi2RW="<<chi2RW<<
 775         "\n";
 776     }
 777     if (reject) continue;
 778
 779
 780     // =========================================
 781     // wirte collected information to histograms
 782     // =========================================
 783
 784     Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
 785     fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
 786
 787
 788     TH3F * his3 = 0;
 789     his3 = (TH3F*)fRMSY->At(padSize);
 790     if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
 791     his3 = (TH3F*)fRMSZ->At(padSize);
 792     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
 793
 794     his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
 795     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
 796     his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
 797     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
 798
 799
 800     his3 = (TH3F*)fResolY->At(padSize);
 801     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
 802     his3 = (TH3F*)fResolZ->At(padSize);
 803     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
 804     his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
 805     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
 806     his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
 807     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
 808     //=============================================================================================
 809     //
 810     // Fill THN histograms
 811     //
 812     Double_t xvar[9];
 813     xvar[1]=padSize;   // pad type
 814     xvar[2]=cluster0->GetZ();  //
 815     xvar[3]=cluster0->GetMax();
 816
 817     xvar[0]=deltay;
 818     xvar[4]=cluster0->GetPad()-Int_t(cluster0->GetPad());// distance to the center of the pad
 819     xvar[5]=angley;
 820     fHisDeltaY->Fill(xvar);
 821     xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
 822     fHisRMSY->Fill(xvar);
 823
 824     xvar[0]=deltaz;
 825     xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin()); // distance to the center of the time bin
 826     xvar[5]=anglez;
 827     fHisDeltaZ->Fill(xvar);
 828     xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
 829     fHisRMSZ->Fill(xvar);
 830
 831   }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
 832 }  // FillResolutionHistoLocal(...)
 833
 834
 835
 836
 837
 838
 839
 840
 841 void  AliTPCcalibTracks::SetStyle() const {
 842    //
 843    // set style, can be called by all draw functions
 844    //
 845    gROOT->SetStyle("Plain");
 846    gStyle->SetFillColor(10);
 847    gStyle->SetPadColor(10);
 848    gStyle->SetCanvasColor(10);
 849    gStyle->SetStatColor(10);
 850    gStyle->SetPalette(1,0);
 851    gStyle->SetNumberContours(60);
 852 }
 853
 854
 855
 856 void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
 857    //
 858    // all functions are called, that produce pictures
 859    // the histograms are written to the directory 'pathName'
 860    // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
 861    // 'stat' is also the number of minEntries for MakeResPlotsQTree
 862    //
 863
 864   if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
 865   MakeResPlotsQTree(stat, pathName);
 866 }
 867
 868
 869
 870
 871 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
 872    //
 873    // Make tree with resolution parameters
 874    // the result is written to 'resol.root' in directory 'pathname'
 875    // file information are available in fileInfo
 876    // available variables in the tree 'Resol':
 877    //  Entries: number of entries for this resolution point
 878    //  nbins:   number of bins that were accumulated
 879    //  Dim:     direction, Dim==0: y-direction, Dim==1: z-direction
 880    //  Pad:     padSize; short, medium and long
 881    //  Length:  pad length, 0.75, 1, 1.5
 882    //  QMean:   mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
 883    //  Zc:      center of middle bin in drift direction
 884    //  Zm:      mean dirftlength for accumulated Delta-Histograms
 885    //  Zs:      width of driftlength bin
 886    //  AngleC:  center of middle bin in Angle-Direction
 887    //  AngleM:  mean angle for accumulated Delta-Histograms
 888    //  AngleS:  width of Angle-bin
 889    //  Resol:   sigma for gaus fit through Delta-Histograms
 890    //  Sigma:   error of sigma for gaus fit through Delta Histograms
 891    //  MeanR:   mean of the Delta-Histogram
 892    //  SigmaR:  rms of the Delta-Histogram
 893    //  RMSm:    mean of the gaus fit through RMS-Histogram
 894    //  RMS:     sigma of the gaus fit through RMS-Histogram
 895    //  RMSe0:   error of mean of gaus fit in RMS-Histogram
 896    //  RMSe1:   error of sigma of gaus fit in RMS-Histogram
 897    //
 898
 899   if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
 900   if (GetDebugLevel() > -1) cout << "    relax, the calculation will take a while..." << endl;
 901
 902    gSystem->MakeDirectory(pathName);
 903    gSystem->ChangeDirectory(pathName);
 904    TString kFileName = "resol.root";
 905    TTreeSRedirector fTreeResol(kFileName.Data());
 906
 907    TH3F *resArray[2][3][11];
 908    TH3F *rmsArray[2][3][11];
 909
 910    // load histograms from fArrayQDY and fArrayQDZ
 911    // into resArray and rmsArray
 912    // that is all we need here
 913    for (Int_t idim = 0; idim < 2; idim++){
 914       for (Int_t ipad = 0; ipad < 3; ipad++){
 915          for (Int_t iq = 0; iq <= 10; iq++){
 916             rmsArray[idim][ipad][iq]=0;
 917             resArray[idim][ipad][iq]=0;
 918             Int_t bin = GetBin(iq,ipad);
 919             TH3F *hresl = 0;
 920             if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
 921             if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
 922             if (!hresl) continue;
 923             resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
 924             resArray[idim][ipad][iq]->SetDirectory(0);
 925             TH3F * hreslRMS = 0;
 926             if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
 927             if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
 928             if (!hreslRMS) continue;
 929             rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
 930             rmsArray[idim][ipad][iq]->SetDirectory(0);
 931          }
 932       }
 933    }
 934    if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
 935
 936    //--------------------------------------------------------------------------------------------
 937
 938    char name[200];
 939    Double_t qMean;
 940    Double_t zMean, angleMean, zCenter, angleCenter;
 941    Double_t zSigma, angleSigma;
 942    TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
 943    TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
 944    TF1 *fitFunction = new TF1("fitFunction", "gaus");
 945    Float_t entriesQ = 0;
 946    Int_t loopCounter = 1;
 947
 948    for (Int_t idim = 0; idim < 2; idim++){
 949       // Loop y-z corrdinate
 950       for (Int_t ipad = 0; ipad < 3; ipad++){
 951          // loop pad type
 952          for (Int_t iq = -1; iq < 10; iq++){
 953             // LOOP Q
 954            if (GetDebugLevel() > -1)
 955                cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
 956                   << (Int_t)((loopCounter)/66.*100) << "% done), "
 957                   << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << "  \r" << std::flush;
 958             loopCounter++;
 959             TH3F *hres = 0;
 960             TH3F *hrms = 0;
 961             qMean    = 0;
 962             entriesQ = 0;
 963
 964             // calculate qMean
 965             if (iq == -1){
 966                // integrated spectra
 967                for (Int_t iql = 0; iql < 10; iql++){
 968                   Int_t bin = GetBin(iql,ipad);
 969                   TH3F *hresl = resArray[idim][ipad][iql];
 970                   TH3F *hrmsl = rmsArray[idim][ipad][iql];
 971                   if (!hresl) continue;
 972                   if (!hrmsl) continue;
 973                   entriesQ += hresl->GetEntries();
 974                   qMean += hresl->GetEntries() * GetQ(bin);
 975                   if (!hres) {
 976                      hres = (TH3F*)hresl->Clone();
 977                      hrms = (TH3F*)hrmsl->Clone();
 978                   }
 979                   else{
 980                      hres->Add(hresl);
 981                      hrms->Add(hrmsl);
 982                   }
 983                }
 984                qMean /= entriesQ;
 985                qMean *= -1.;  // integral mean charge
 986             }
 987             else {
 988                // loop over neighboured Q-bins
 989                // accumulate entries from neighboured Q-bins
 990                for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
 991                   if (iql < 0) continue;
 992                   Int_t bin = GetBin(iql,ipad);
 993                   TH3F * hresl = resArray[idim][ipad][iql];
 994                   TH3F * hrmsl = rmsArray[idim][ipad][iql];
 995                   if (!hresl) continue;
 996                   if (!hrmsl) continue;
 997                   entriesQ += hresl->GetEntries();
 998                   qMean += hresl->GetEntries() * GetQ(bin);
 999                   if (!hres) {
1000                      hres = (TH3F*) hresl->Clone();
1001                      hrms = (TH3F*) hrmsl->Clone();
1002                   }
1003                   else{
1004                      hres->Add(hresl);
1005                      hrms->Add(hrmsl);
1006                   }
1007                }
1008                qMean/=entriesQ;
1009             }
1010             if (!hres) continue;
1011             if (!hrms) continue;
1012
1013             TAxis *xAxisDriftLength = hres->GetXaxis();   // driftlength / z - axis
1014             TAxis *yAxisAngle       = hres->GetYaxis();   // angle axis
1015             TAxis *zAxisDelta       = hres->GetZaxis();   // delta axis
1016             TAxis *zAxisRms         = hrms->GetZaxis();   // rms axis
1017
1018             // loop over all angle bins
1019             for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1020                angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1021                // loop over all driftlength bins
1022                for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1023                   zCenter    = xAxisDriftLength->GetBinCenter(ibinxDL);
1024                   zSigma     = xAxisDriftLength->GetBinWidth(ibinxDL);
1025                   angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
1026                   zMean      = zCenter;      // changens, when more statistic is accumulated
1027                   angleMean  = angleCenter;  // changens, when more statistic is accumulated
1028
1029                   // create 2 1D-Histograms, projectionRes and projectionRms
1030                   // these histograms are delta histograms for given direction, padSize, chargeBin,
1031                   // angleBin and driftLengthBin
1032                   // later on they will be fitted with a gausian, its sigma is the resoltuion...
1033                   snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
1034                   // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1035                   projectionRes->SetNameTitle(name, name);
1036                   snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
1037                   // TH1D * projectionRms =  new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1038                   projectionRms->SetNameTitle(name, name);
1039
1040                   projectionRes->Reset();
1041                   projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1042                   projectionRms->Reset();
1043                   projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1044                   projectionRes->SetDirectory(0);
1045                   projectionRms->SetDirectory(0);
1046
1047                   Double_t entries = 0;
1048                   Int_t    nbins   = 0;   // counts, how many bins were accumulated
1049                   zMean     = 0;
1050                   angleMean = 0;
1051                   entries   = 0;
1052
1053                   // fill projectionRes and projectionRms for given dim, ipad and iq,
1054                   // as well as for given angleBin and driftlengthBin
1055                   // if this gives not enough statistic, include neighbourhood
1056                   // (angle and driftlength) successifely
1057                   for (Int_t dbin = 0; dbin <= 8; dbin++){              // delta-bins around centered angleBin and driftlengthBin
1058                      for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) {   // delta-bins in angle direction
1059                         for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
1060                            if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue;   // add each bin only one time !
1061                            Int_t binx2 = ibinxDL + dbinx2;                       // position variable in x (driftlength) direction
1062                            Int_t biny2 = ibinyAngle + dbiny2;                    // position variable in y (angle)  direction
1063                            if (binx2 < 1 || biny2 < 1) continue;                 // don't go out of the histogram!
1064                            if (binx2 >= xAxisDriftLength->GetNbins()) continue;  // don't go out of the histogram!
1065                            if (biny2 >= yAxisAngle->GetNbins()) continue;        // don't go out of the histogram!
1066                            nbins++;                                              // count the number of accumulated bins
1067                            // Fill resolution histo
1068                            for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
1069                               // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3);     // unused variable
1070                               projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
1071                               entries   += hres->GetBinContent(binx2, biny2, ibin3);
1072                               zMean     += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
1073                               angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
1074                            }  // ibin3 loop
1075                            // fill RMS histo
1076                            for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
1077                               projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
1078                            }
1079                         }  //dbinx2 loop
1080                         if (entries > minEntries) break; // enough statistic accumulated
1081                      }  // dbiny2 loop
1082                      if (entries > minEntries) break;    // enough statistic accumulated
1083                   }  // dbin loop
1084                   if ( entries< minEntries) continue;  // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
1085                   zMean /= entries;
1086                   angleMean /= entries;
1087
1088                   if (entries > minEntries) {
1089                      //  when enough statistic is accumulated
1090                      //  fit Delta histograms with a gausian
1091                      //  of the gausian is the resolution (resol), its fit error is sigma
1092                      //  store also mean and RMS of the histogram
1093                      Float_t xmin     = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1094                      Float_t xmax     = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1095
1096 //                      projectionRes->Fit("gaus", "q0", "", xmin, xmax);
1097 //                      Float_t resol    = projectionRes->GetFunction("gaus")->GetParameter(2);
1098 //                      Float_t sigma    = projectionRes->GetFunction("gaus")->GetParError(2);
1099                      fitFunction->SetMaximum(xmax);
1100                      fitFunction->SetMinimum(xmin);
1101                      projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
1102                      Float_t resol    = fitFunction->GetParameter(2);
1103                      Float_t sigma    = fitFunction->GetParError(2);
1104
1105                      Float_t meanR    = projectionRes->GetMean();
1106                      Float_t sigmaR   = projectionRes->GetRMS();
1107                      // fit also RMS histograms with a gausian
1108                      // store mean and sigma of the gausian in rmsMean and rmsSigma
1109                      // store also the fit errors in errorRMS and errorSigma
1110                      xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1111                      xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1112
1113 //                      projectionRms->Fit("gaus","q0","",xmin,xmax);
1114 //                      Float_t rmsMean    = projectionRms->GetFunction("gaus")->GetParameter(1);
1115 //                      Float_t rmsSigma   = projectionRms->GetFunction("gaus")->GetParameter(2);
1116 //                      Float_t errorRMS   = projectionRms->GetFunction("gaus")->GetParError(1);
1117 //                      Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
1118                      projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
1119                      Float_t rmsMean    = fitFunction->GetParameter(1);
1120                      Float_t rmsSigma   = fitFunction->GetParameter(2);
1121                      Float_t errorRMS   = fitFunction->GetParError(1);
1122                      Float_t errorSigma = fitFunction->GetParError(2);
1123
1124                      Float_t length = 0.75;
1125                      if (ipad == 1) length = 1;
1126                      if (ipad == 2) length = 1.5;
1127
1128                      fTreeResol<<"Resol"<<
1129                         "Entries="<<entries<<      // number of entries for this resolution point
1130                         "nbins="<<nbins<<          // number of bins that were accumulated
1131                         "Dim="<<idim<<             // direction, Dim==0: y-direction, Dim==1: z-direction
1132                         "Pad="<<ipad<<             // padSize; short, medium and long
1133                         "Length="<<length<<        // pad length, 0.75, 1, 1.5
1134                         "QMean="<<qMean<<          // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1135                         "Zc="<<zCenter<<           // center of middle bin in drift direction
1136                         "Zm="<<zMean<<             // mean dirftlength for accumulated Delta-Histograms
1137                         "Zs="<<zSigma<<            // width of driftlength bin
1138                         "AngleC="<<angleCenter<<   // center of middle bin in Angle-Direction
1139                         "AngleM="<<angleMean<<     // mean angle for accumulated Delta-Histograms
1140                         "AngleS="<<angleSigma<<    // width of Angle-bin
1141                         "Resol="<<resol<<          // sigma for gaus fit through Delta-Histograms
1142                         "Sigma="<<sigma<<          // error of sigma for gaus fit through Delta Histograms
1143                         "MeanR="<<meanR<<          // mean of the Delta-Histogram
1144                         "SigmaR="<<sigmaR<<        // rms of the Delta-Histogram
1145                         "RMSm="<<rmsMean<<         // mean of the gaus fit through RMS-Histogram
1146                         "RMSs="<<rmsSigma<<        // sigma of the gaus fit through RMS-Histogram
1147                         "RMSe0="<<errorRMS<<       // error of mean of gaus fit in RMS-Histogram
1148                         "RMSe1="<<errorSigma<<     // error of sigma of gaus fit in RMS-Histogram
1149                         "\n";
1150                      if (GetDebugLevel() > 5) {
1151                         projectionRes->SetDirectory(fTreeResol.GetFile());
1152                         projectionRes->Write(projectionRes->GetName());
1153                         projectionRes->SetDirectory(0);
1154                         projectionRms->SetDirectory(fTreeResol.GetFile());
1155                         projectionRms->Write(projectionRms->GetName());
1156                         projectionRes->SetDirectory(0);
1157                      }
1158                   }  // if (projectionRes->GetSum() > minEntries)
1159                }  // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
1160             }  // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
1161
1162          }  // iq-loop
1163       }  // ipad-loop
1164    }  // idim-loop
1165    delete projectionRes;
1166    delete projectionRms;
1167
1168 //    TFile resolFile(fTreeResol.GetFile());
1169    TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
1170    fileInfo.Write("fileInfo");
1171 //    resolFile.Close();
1172 //    fTreeResol.GetFile()->Close();
1173    if (GetDebugLevel() > -1) cout << endl;
1174    if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
1175    gSystem->ChangeDirectory("..");
1176 }
1177
1178
1179
1180
1181
1182 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
1183    //
1184    // function to merge several AliTPCcalibTracks objects after PROOF calculation
1185    // The object's histograms are merged via their merge functions
1186    // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
1187    //
1188
1189   if (GetDebugLevel() > 0) cout << " *****  this is AliTPCcalibTracks::Merge(TCollection *collectionList)  *****"<< endl;
1190    if (!collectionList) return 0;
1191    if (collectionList->IsEmpty()) return -1;
1192
1193    if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl;     //    REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
1194    if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
1195    collectionList->Print();
1196
1197    // create a list for each data member
1198    TList* deltaYList = new TList;
1199    TList* deltaZList = new TList;
1200    TList* arrayAmpRowList = new TList;
1201    TList* rejectedTracksList = new TList;
1202    TList* clusterCutHistoList = new TList;
1203    TList* arrayAmpList = new TList;
1204    TList* arrayQDYList = new TList;
1205    TList* arrayQDZList = new TList;
1206    TList* arrayQRMSYList = new TList;
1207    TList* arrayQRMSZList = new TList;
1208    TList* resolYList = new TList;
1209    TList* resolZList = new TList;
1210    TList* rMSYList = new TList;
1211    TList* rMSZList = new TList;
1212
1213 //    TList* nRowsList = new TList;
1214 //    TList* nSectList = new TList;
1215 //    TList* fileNoList = new TList;
1216
1217    TIterator *listIterator = collectionList->MakeIterator();
1218    AliTPCcalibTracks *calibTracks = 0;
1219    if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
1220    Int_t counter = 0;
1221    while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
1222       // loop over all entries in the collectionList and get dataMembers into lists
1223
1224       arrayQDYList->Add(calibTracks->GetfArrayQDY());
1225       arrayQDZList->Add(calibTracks->GetfArrayQDZ());
1226       arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
1227       arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
1228       resolYList->Add(calibTracks->GetfResolY());
1229       resolZList->Add(calibTracks->GetfResolZ());
1230       rMSYList->Add(calibTracks->GetfRMSY());
1231       rMSZList->Add(calibTracks->GetfRMSZ());
1232       rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
1233       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
1234       //
1235       if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
1236         fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
1237       //      fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
1238       counter++;
1239       if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
1240       AddHistos(calibTracks);
1241    }
1242
1243
1244    // merge data members
1245    if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
1246    fClusterCutHisto->Merge(clusterCutHistoList);
1247    fRejectedTracksHisto->Merge(rejectedTracksList);
1248
1249    TObjArray* objarray = 0;
1250    TH1* hist = 0;
1251    TList* histList = 0;
1252    TIterator *objListIterator = 0;
1253
1254
1255    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
1256    // merge fArrayQDY
1257    for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
1258       objListIterator = arrayQDYList->MakeIterator();
1259       histList = new TList;
1260       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1261          // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
1262          hist = (TH3F*)(objarray->At(i));
1263          histList->Add(hist);
1264       }
1265       ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
1266       delete histList;
1267       delete objListIterator;
1268    }
1269
1270    if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
1271    // merge fArrayQDZ
1272    for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1273       objListIterator = arrayQDZList->MakeIterator();
1274       histList = new TList;
1275       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1276          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1277          hist = (TH3F*)(objarray->At(i));
1278          histList->Add(hist);
1279       }
1280       ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
1281       delete histList;
1282       delete objListIterator;
1283    }
1284
1285    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
1286    // merge fArrayQRMSY
1287    for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
1288       objListIterator = arrayQRMSYList->MakeIterator();
1289       histList = new TList;
1290       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1291          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1292          //   if (!objarray) continue; // removed for coverity -> JMT
1293          hist = (TH3F*)(objarray->At(i));
1294          histList->Add(hist);
1295       }
1296       ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
1297       delete histList;
1298       delete objListIterator;
1299    }
1300
1301    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
1302    // merge fArrayQRMSZ
1303    for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1304       objListIterator = arrayQRMSZList->MakeIterator();
1305       histList = new TList;
1306       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1307          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1308          hist = (TH3F*)(objarray->At(i));
1309          histList->Add(hist);
1310       }
1311       ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
1312       delete histList;
1313       delete objListIterator;
1314    }
1315
1316
1317
1318
1319
1320
1321    if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
1322    // merge fResolY
1323    for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
1324       objListIterator = resolYList->MakeIterator();
1325       histList = new TList;
1326       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1327          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1328          hist = (TH3F*)(objarray->At(i));
1329          histList->Add(hist);
1330       }
1331       ((TH3F*)(fResolY->At(i)))->Merge(histList);
1332       delete histList;
1333       delete objListIterator;
1334    }
1335
1336    // merge fResolZ
1337    for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1338       objListIterator = resolZList->MakeIterator();
1339       histList = new TList;
1340       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1341          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1342          hist = (TH3F*)(objarray->At(i));
1343          histList->Add(hist);
1344       }
1345       ((TH3F*)(fResolZ->At(i)))->Merge(histList);
1346       delete histList;
1347       delete objListIterator;
1348    }
1349
1350    // merge fRMSY
1351    for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
1352       objListIterator = rMSYList->MakeIterator();
1353       histList = new TList;
1354       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1355          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1356          hist = (TH3F*)(objarray->At(i));
1357          histList->Add(hist);
1358       }
1359       ((TH3F*)(fRMSY->At(i)))->Merge(histList);
1360       delete histList;
1361       delete objListIterator;
1362    }
1363
1364    // merge fRMSZ
1365    for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1366       objListIterator = rMSZList->MakeIterator();
1367       histList = new TList;
1368       while (( objarray =  (TObjArray*)objListIterator->Next() )) {
1369          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1370          hist = (TH3F*)(objarray->At(i));
1371          histList->Add(hist);
1372       }
1373       ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
1374       delete histList;
1375       delete objListIterator;
1376    }
1377
1378    delete deltaYList;
1379    delete deltaZList;
1380    delete arrayAmpRowList;
1381    delete arrayAmpList;
1382    delete arrayQDYList;
1383    delete arrayQDZList;
1384    delete arrayQRMSYList;
1385    delete arrayQRMSZList;
1386    delete resolYList;
1387    delete resolZList;
1388    delete rMSYList;
1389    delete rMSZList;
1390    delete listIterator;
1391
1392    if (GetDebugLevel() > 0) cout << "merging done!" << endl;
1393
1394    return 1;
1395 }
1396
1397
1398
1399
1400 void AliTPCcalibTracks::MakeHistos(){
1401   //
1402   ////make THnSparse
1403   //
1404   //THnSparse  *fHisDeltaY;    // THnSparse - delta Y
1405   //THnSparse  *fHisDeltaZ;    // THnSparse - delta Z
1406   //THnSparse  *fHisRMSY;      // THnSparse - rms Y
1407   //THnSparse  *fHisRMSZ;      // THnSparse - rms Z
1408   //THnSparse  *fHisQmax;      // THnSparse - qmax
1409   //THnSparse  *fHisQtot;      // THnSparse - qtot
1410   // cluster  performance bins
1411   // 0 - variable of interest
1412   // 1 - pad type   - 0- short 1-medium 2-long pads
1413   // 2 - drift length - drift length -0-1
1414   // 3 - Qmax         - Qmax  - 2- 400
1415   // 4 - cog          - COG position - 0-1
1416   // 5 - tan(phi)     - local y angle
1417   // 6 - tan(theta)   - local z angle
1418   // 7 - sector       - sector number
1419   Double_t xminTrack[8], xmaxTrack[8];
1420   Int_t binsTrack[8];
1421   TString axisName[8];
1422
1423   //
1424   binsTrack[0]=200;
1425   axisName[0]  ="var";
1426
1427   binsTrack[1] =3;
1428   xminTrack[1] =0; xmaxTrack[1]=3;
1429   axisName[1]  ="pad type";
1430   //
1431   binsTrack[2] =20;
1432   xminTrack[2] =-250; xmaxTrack[2]=250;
1433   axisName[2]  ="z";
1434   //
1435   binsTrack[3] =10;
1436   xminTrack[3] =1; xmaxTrack[3]=400;
1437   axisName[3]  ="Qmax";
1438   //
1439   binsTrack[4] =20;
1440   xminTrack[4] =0; xmaxTrack[4]=1;
1441   axisName[4]  ="cog";
1442   //
1443   binsTrack[5] =15;
1444   xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
1445   axisName[5]  ="tan(angle)";
1446   //
1447   //
1448   xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1449   fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1450   xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1451   fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1452   xminTrack[0] =0.; xmaxTrack[0]=0.5;
1453   fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1454   xminTrack[0] =0.; xmaxTrack[0]=0.5;
1455   fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1456   xminTrack[0] =0.; xmaxTrack[0]=100;
1457   fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1458
1459   xminTrack[0] =0.; xmaxTrack[0]=250;
1460   fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1461
1462
1463   for (Int_t ivar=0;ivar<6;ivar++){
1464     fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1465     fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1466     fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1467     fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1468     fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1469     fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1470     fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
1471     fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1472     fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1473     fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1474     fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1475     fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1476   }
1477
1478
1479   BinLogX(fHisDeltaY,3);
1480   BinLogX(fHisDeltaZ,3);
1481   BinLogX(fHisRMSY,3);
1482   BinLogX(fHisRMSZ,3);
1483   BinLogX(fHisQmax,3);
1484   BinLogX(fHisQtot,3);
1485
1486 }
1487
1488 void    AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
1489   //
1490   // Add histograms
1491   //
1492   if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
1493   if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
1494   if (calib->fHisRMSY)   fHisRMSY->Add(calib->fHisRMSY);
1495   if (calib->fHisRMSZ)   fHisRMSZ->Add(calib->fHisRMSZ);
1496 }
1497
1498
1499
1500 void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
1501   //
1502   // Dump summary info
1503   //
1504   //  0.OBJ: TAxis     var
1505   //  1.OBJ: TAxis     pad type
1506   //  2.OBJ: TAxis     z
1507   //  3.OBJ: TAxis     Qmax
1508   //  4.OBJ: TAxis     cog
1509   //  5.OBJ: TAxis     tan(angle)
1510   //
1511   if (ptype>3) return;
1512   Int_t idim[6]={0,1,2,3,4,5};
1513   TString hname[4]={"dy","dz","rmsy","rmsz"};
1514   //
1515   Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
1516   Int_t first5=hisInput->GetAxis(5)->GetFirst();
1517   Int_t last5 =hisInput->GetAxis(5)->GetLast();
1518   //
1519   for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){   // axis 5 - local angle
1520     Bool_t bin5All=kTRUE;
1521     if (ibin5>=first5){
1522       hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
1523       if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
1524       bin5All=kFALSE;
1525     }
1526     Double_t      x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
1527     THnSparse * his5= hisInput->Projection(5,idim);         //projected histogram according selection 5
1528     Int_t nbins4=his5->GetAxis(4)->GetNbins();
1529     Int_t first4=his5->GetAxis(4)->GetFirst();
1530     Int_t last4 =his5->GetAxis(4)->GetLast();
1531
1532     for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){   // axis 4 - cog
1533       Bool_t bin4All=kTRUE;
1534       if (ibin4>=first4){
1535         his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
1536         if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
1537         bin4All=kFALSE;
1538       }
1539       Double_t      x4= his5->GetAxis(4)->GetBinCenter(ibin4);
1540       THnSparse * his4= his5->Projection(4,idim);         //projected histogram according selection 4
1541       //
1542       Int_t nbins3=his4->GetAxis(3)->GetNbins();
1543       Int_t first3=his4->GetAxis(3)->GetFirst();
1544       Int_t last3 =his4->GetAxis(3)->GetLast();
1545       //
1546       for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){   // axis 3 - Qmax
1547         Bool_t bin3All=kTRUE;
1548         if (ibin3>=first3){
1549           his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
1550           bin3All=kFALSE;
1551         }
1552         Double_t      x3= his4->GetAxis(3)->GetBinCenter(ibin3);
1553         THnSparse * his3= his4->Projection(3,idim);         //projected histogram according selection 3
1554         //
1555         Int_t nbins2    = his3->GetAxis(2)->GetNbins();
1556         Int_t first2    = his3->GetAxis(2)->GetFirst();
1557         Int_t last2     = his3->GetAxis(2)->GetLast();
1558         //
1559         for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){   // axis 2 - z
1560           Bool_t bin2All=kTRUE;
1561           Double_t      x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1562           if (ibin2>=first2){
1563             his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
1564             if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
1565             bin2All=kFALSE;
1566           }
1567           THnSparse * his2= his3->Projection(2,idim);         //projected histogram according selection 2
1568           //
1569           Int_t nbins1     = his2->GetAxis(1)->GetNbins();
1570           Int_t first1     = his2->GetAxis(1)->GetFirst();
1571           Int_t last1      = his2->GetAxis(1)->GetLast();
1572           for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){   //axis 1 - pad type
1573             Bool_t bin1All=kTRUE;
1574             if (ibin1>=first1){
1575               his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
1576               bin1All=kFALSE;
1577             }
1578             Double_t       x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
1579             TH1 * hisDelta = his2->Projection(0);
1580             Double_t entries = hisDelta->GetEntries();
1581             Double_t mean=0, rms=0;
1582             if (entries>10){
1583               mean    = hisDelta->GetMean();
1584               rms = hisDelta->GetRMS();
1585               hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
1586               mean    = hisDelta->GetMean();
1587               rms = hisDelta->GetRMS();
1588             }
1589             //
1590             (*pcstream)<<hname[ptype].Data()<<
1591               // flag - intgrated values for given bin
1592               "angleA="<<bin5All<<
1593               "cogA="<<bin4All<<
1594               "qmaxA="<<bin3All<<
1595               "zA="<<bin2All<<
1596               "ipadA="<<bin1All<<
1597               // center of bin value
1598               "angle="<<x5<<       // local angle
1599               "cog="<<x4<<         // distance cluster to center
1600               "qmax="<<x3<<        // qmax
1601               "z="<<x2<<           // z of the cluster
1602               "ipad="<<x1<<        // type of the region
1603               // mean values
1604               //
1605               "entries="<<entries<<
1606               "mean="<<mean<<
1607               "rms="<<rms<<
1608               "ptype="<<ptype<<   //
1609               "\n";
1610             delete hisDelta;
1611             printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);
1612           }//loop z
1613           delete his2;
1614         } // loop Q max
1615         delete his3;
1616       } // loop COG
1617       delete his4;
1618     }//loop local angle
1619     delete his5;
1620   }
1621 }