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