MUON/MUONRecoCheck.C

   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 // $Id$
  17
  18 /// \ingroup macros
  19 /// \file MUONRecoCheck.C
  20 /// \brief Utility macro to check the muon reconstruction.
  21 ///
  22 /// Reconstructed tracks are compared to reference tracks. The reference tracks
  23 /// are built from AliTrackReference for the hit in chamber (0..9) and from
  24 /// kinematics (TreeK) for the vertex parameters.
  25 ///
  26 /// \author Jean-Pierre Cussonneau, Philippe Pillot, Subatech
  27
  28 // ROOT includes
  29 #include <Riostream.h>
  30 #include "TMath.h"
  31 #include "TClonesArray.h"
  32 #include "TH1.h"
  33 #include "TH2.h"
  34 #include "TH3.h"
  35 #include "TGraphErrors.h"
  36 #include "TGraphAsymmErrors.h"
  37 #include "TF1.h"
  38 #include "TFile.h"
  39 #include "TCanvas.h"
  40 #include "TLegend.h"
  41 #include "TGeoManager.h"
  42
  43 // STEER includes
  44 #include "AliCDBManager.h"
  45 #include "AliGeomManager.h"
  46 #include "AliLog.h"
  47
  48 // MUON includes
  49 #include "AliMUONCDB.h"
  50 #include "AliMUONConstants.h"
  51 #include "AliMUONTrack.h"
  52 #include "AliMUONRecoCheck.h"
  53 #include "AliMUONTrackParam.h"
  54 #include "AliMUONRecoParam.h"
  55 #include "AliMUONVTrackStore.h"
  56 #include "AliMUONVCluster.h"
  57 #include "AliMUONTrackExtrap.h"
  58 #include "AliMUONESDInterface.h"
  59 #include "AliMUONVTriggerTrackStore.h"
  60 #include "AliMUONTriggerTrack.h"
  61
  62 Double_t langaufun(Double_t *x, Double_t *par);
  63 void     FitGausResVsMom(TH2* h, Int_t nBins, const Double_t mean0, const Double_t sigma0, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma);
  64 void     FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gSigma);
  65 TCanvas* DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2);
  66 TCanvas* DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3);
  67 TCanvas* DrawResMomVsMom(const char* name, const char* title, TH2* h, Int_t nBins, TF1* f2 = 0x0, const char* fitting = "");
  68
  69 //------------------------------------------------------------------------------------
  70 void MUONRecoCheck (Int_t nEvent = -1, const char* pathSim="./generated/", const char* esdFileName="AliESDs.root",
  71                     const char* ocdbPath = "local://$ALICE_ROOT/OCDB", Int_t absorberRegion = -1)
  72 {
  73   /// Associate the reconstructed tracks with the simulated ones and check the quality of the reconstruction
  74   /// (tracking/trigger efficiency; momentum, slope,... resolutions at first cluster and at vertex; cluster resolution).
  75   /// You can limit the calculation of track resolution at vertex to the tracks crossing the absorber in a given region
  76   /// with the flag "absorberRegion": -1=all, 1=[2,3]deg, 2=[3,10]deg.
  77
  78   Double_t aAbsLimits[2];
  79   if (absorberRegion > -1) {
  80     if (absorberRegion == 1) {
  81       aAbsLimits[0] = 2.;
  82       aAbsLimits[1] = 3.;
  83     } else if (absorberRegion == 2) {
  84       aAbsLimits[0] = 3.;
  85       aAbsLimits[1] = 10.;
  86     } else {
  87       cout<<"Unknown absorber region. Valid choices are: -1=all, 1=[2,3]deg, 2=[3,10]deg"<<endl;
  88       return;
  89     }
  90   } else {
  91     aAbsLimits[0] = 0.;
  92     aAbsLimits[1] = 90.;
  93   }
  94
  95   AliLog::SetClassDebugLevel("AliMCEvent",-1);
  96
  97   // ###################################### define histograms ###################################### //
  98   // File for histograms and histogram booking
  99   TFile *histoFile = new TFile("MUONRecoCheck.root", "RECREATE");
 100
 101   TH1F *hReconstructible = new TH1F("hReconstructible"," Nb of reconstructible tracks / evt",15,-0.5,14.5);
 102   TH1F *hReco = new TH1F("hReco"," Nb of reconstructed tracks / evt",15,-0.5,14.5);
 103   TH1F *hNClusterComp = new TH1F("hNClusterComp"," Nb of compatible clusters / track ",15,-0.5,14.5);
 104   TH1F *hTrackRefID = new TH1F("hTrackRefID"," track reference ID ",100,-0.5,99.5);
 105   TH1F *hTriggerable = new TH1F("hTriggerable"," Nb of triggerable tracks / evt",15,-0.5,14.5);
 106   TH1F *hTriggered = new TH1F("hTriggered"," Nb of triggered tracks / evt",15,-0.5,14.5);
 107
 108   // momentum resolution at vertex
 109   histoFile->mkdir("momentumAtVertex","momentumAtVertex");
 110   histoFile->cd("momentumAtVertex");
 111
 112   const Int_t pNBins = 30;
 113   const Double_t pEdges[2] = {0., 300.};
 114   const Int_t deltaPAtVtxNBins = 250;
 115   const Double_t deltaPAtVtxEdges[2] = {-35., 15.};
 116
 117   TH1F *hResMomVertex = new TH1F("hResMomVertex"," delta P at vertex;#Delta_{p} (GeV/c)",deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 118
 119   TH2D *hResMomVertexVsMom = new TH2D("hResMomVertexVsMom","#Delta_{p} at vertex versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 120   TH2D *hResMomVertexVsMom_2_3_Deg = new TH2D("hResMomVertexVsMom_2_3_Deg","#Delta_{p} at vertex versus p for tracks between 2 and 3 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 121   TH2D *hResMomVertexVsMom_3_10_Deg = new TH2D("hResMomVertexVsMom_3_10_Deg","#Delta_{p} at vertex versus p for tracks between 3 and 10 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 122   TH2D *hResMomVertexVsMom_0_2_DegMC = new TH2D("hResMomVertexVsMom_0_2_DegMC","#Delta_{p} at vertex versus p for tracks with MC angle below 2 degrees;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins/10,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 123
 124   TH2D *hResMomVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResMomVertexVsPosAbsEnd_0_2_DegMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 125   TH2D *hResMomVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResMomVertexVsPosAbsEnd_2_3_DegMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 126   TH2D *hResMomVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResMomVertexVsPosAbsEnd_3_10_DegMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 127
 128   TH2D *hResMomVertexVsAngle = new TH2D("hResMomVertexVsAngle","#Delta_{p} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 129   TH2D *hResMomVertexVsMCAngle = new TH2D("hResMomVertexVsMCAngle","#Delta_{p} at vertex versus MC angle;MC angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 130   TH3D *hResMomVertexVsAngleVsMom = new TH3D("hResMomVertexVsAngleVsMom","#Delta_{p} at vertex versus track position at absorber end converted to degrees versus momentum;p (GeV/c);angle (Deg);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],100,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
 131
 132   TGraphAsymmErrors* gMeanResMomVertexVsMom = new TGraphAsymmErrors(pNBins);
 133   gMeanResMomVertexVsMom->SetName("gMeanResMomVertexVsMom");
 134   gMeanResMomVertexVsMom->SetTitle("<#Delta_{p}> at vertex versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
 135   TGraphAsymmErrors* gMostProbResMomVertexVsMom = new TGraphAsymmErrors(pNBins);
 136   gMostProbResMomVertexVsMom->SetName("gMostProbResMomVertexVsMom");
 137   gMostProbResMomVertexVsMom->SetTitle("Most probable #Delta_{p} at vertex versus p;p (GeV/c);Most prob. #Delta_{p} (GeV/c)");
 138   TGraphAsymmErrors* gSigmaResMomVertexVsMom = new TGraphAsymmErrors(pNBins);
 139   gSigmaResMomVertexVsMom->SetName("gSigmaResMomVertexVsMom");
 140   gSigmaResMomVertexVsMom->SetTitle("#sigma_{p}/p at vertex versus p;p (GeV/c);#sigma_{p}/p (%)");
 141
 142   // momentum resolution at first cluster
 143   histoFile->mkdir("momentumAtFirstCluster","momentumAtFirstCluster");
 144   histoFile->cd("momentumAtFirstCluster");
 145
 146   const Int_t deltaPAtFirstClNBins = 500;
 147   const Double_t deltaPAtFirstClEdges[2] = {-25., 25.};
 148
 149   TH1F *hResMomFirstCluster = new TH1F("hResMomFirstCluster"," delta P at first cluster;#Delta_{p} (GeV/c)",deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
 150   TH2D *hResMomFirstClusterVsMom = new TH2D("hResMomFirstClusterVsMom","#Delta_{p} at first cluster versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
 151
 152   TGraphAsymmErrors* gMeanResMomFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
 153   gMeanResMomFirstClusterVsMom->SetName("gMeanResMomFirstClusterVsMom");
 154   gMeanResMomFirstClusterVsMom->SetTitle("<#Delta_{p}> at first cluster versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
 155   TGraphAsymmErrors* gSigmaResMomFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
 156   gSigmaResMomFirstClusterVsMom->SetName("gSigmaResMomFirstClusterVsMom");
 157   gSigmaResMomFirstClusterVsMom->SetTitle("#sigma_{p}/p at first cluster versus p;p (GeV/c);#sigma_{p}/p (%)");
 158
 159   // angular resolution at vertex
 160   histoFile->mkdir("slopesAtVertex","slopesAtVertex");
 161   histoFile->cd("slopesAtVertex");
 162
 163   const Int_t deltaSlopeAtVtxNBins = 500;
 164   const Double_t deltaSlopeAtVtxEdges[2] = {-0.05, 0.05};
 165
 166   TH1F *hResSlopeXVertex = new TH1F("hResSlopeXVertex","#Delta_{slope_{X}} at vertex;#Delta_{slope_{X}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 167   TH1F *hResSlopeYVertex = new TH1F("hResSlopeYVertex","#Delta_{slope_{Y}} at vertex;#Delta_{slope_{Y}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 168   TH2D *hResSlopeXVertexVsMom = new TH2D("hResSlopeXVertexVsMom","#Delta_{slope_{X}} at vertex versus p;p (GeV/c);#Delta_{slope_{X}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 169   TH2D *hResSlopeYVertexVsMom = new TH2D("hResSlopeYVertexVsMom","#Delta_{slope_{Y}} at vertex versus p;p (GeV/c);#Delta_{slope_{Y}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 170
 171   TH2D *hResSlopeXVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResSlopeXVertexVsPosAbsEnd_0_2_DegMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 172   TH2D *hResSlopeYVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResSlopeYVertexVsPosAbsEnd_0_2_DegMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 173   TH2D *hResSlopeXVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResSlopeXVertexVsPosAbsEnd_2_3_DegMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 174   TH2D *hResSlopeYVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResSlopeYVertexVsPosAbsEnd_2_3_DegMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 175   TH2D *hResSlopeXVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResSlopeXVertexVsPosAbsEnd_3_10_DegMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 176   TH2D *hResSlopeYVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResSlopeYVertexVsPosAbsEnd_3_10_DegMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 177
 178   TH2D *hResSlopeXVertexVsAngle = new TH2D("hResSlopeXVertexVsAngle","#Delta_{slope_{X}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 179   TH2D *hResSlopeYVertexVsAngle = new TH2D("hResSlopeYVertexVsAngle","#Delta_{slope_{Y}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 180   TH2D *hResSlopeXVertexVsMCAngle = new TH2D("hResSlopeXVertexVsMCAngle","#Delta_{slope_{X}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 181   TH2D *hResSlopeYVertexVsMCAngle = new TH2D("hResSlopeYVertexVsMCAngle","#Delta_{slope_{Y}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
 182
 183   TGraphAsymmErrors* gMeanResSlopeXVertexVsMom = new TGraphAsymmErrors(pNBins);
 184   gMeanResSlopeXVertexVsMom->SetName("gMeanResSlopeXVertexVsMom");
 185   gMeanResSlopeXVertexVsMom->SetTitle("<#Delta_{slope_{X}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{X}}>");
 186   TGraphAsymmErrors* gSigmaResSlopeXVertexVsMom = new TGraphAsymmErrors(pNBins);
 187   gSigmaResSlopeXVertexVsMom->SetName("gSigmaResSlopeXVertexVsMom");
 188   gSigmaResSlopeXVertexVsMom->SetTitle("#sigma_{slope_{X}} at vertex versus p;p (GeV/c);#sigma_{slope_{X}}");
 189   TGraphAsymmErrors* gMeanResSlopeYVertexVsMom = new TGraphAsymmErrors(pNBins);
 190   gMeanResSlopeYVertexVsMom->SetName("gMeanResSlopeYVertexVsMom");
 191   gMeanResSlopeYVertexVsMom->SetTitle("<#Delta_{slope_{Y}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
 192   TGraphAsymmErrors* gSigmaResSlopeYVertexVsMom = new TGraphAsymmErrors(pNBins);
 193   gSigmaResSlopeYVertexVsMom->SetName("gSigmaResSlopeYVertexVsMom");
 194   gSigmaResSlopeYVertexVsMom->SetTitle("#sigma_{slope_{Y}} at vertex versus p;p (GeV/c);#sigma_{slope_{Y}}");
 195
 196   // angular resolution at first cluster
 197   histoFile->mkdir("slopesAtFirstCluster","slopesAtFirstCluster");
 198   histoFile->cd("slopesAtFirstCluster");
 199
 200   const Int_t deltaSlopeAtFirstClNBins = 500;
 201   const Double_t deltaSlopeAtFirstClEdges[2] = {-0.01, 0.01};
 202
 203   TH1F *hResSlopeXFirstCluster = new TH1F("hResSlopeXFirstCluster","#Delta_{slope_{X}} at first cluster;#Delta_{slope_{X}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
 204   TH2D *hResSlopeXFirstClusterVsMom = new TH2D("hResSlopeXFirstClusterVsMom","#Delta_{slope_{X}} at first cluster versus p;p (GeV/c);#Delta_{slope_{X}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
 205   TH1F *hResSlopeYFirstCluster = new TH1F("hResSlopeYFirstCluster","#Delta_{slope_{Y}} at first cluster;#Delta_{slope_{Y}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
 206   TH2D *hResSlopeYFirstClusterVsMom = new TH2D("hResSlopeYFirstClusterVsMom","#Delta_{slope_{Y}} at first cluster versus p;p (GeV/c);#Delta_{slope_{Y}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
 207
 208   TGraphAsymmErrors* gMeanResSlopeXFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
 209   gMeanResSlopeXFirstClusterVsMom->SetName("gMeanResSlopeXFirstClusterVsMom");
 210   gMeanResSlopeXFirstClusterVsMom->SetTitle("<#Delta_{slope_{X}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{X}}>");
 211   TGraphAsymmErrors* gSigmaResSlopeXFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
 212   gSigmaResSlopeXFirstClusterVsMom->SetName("gSigmaResSlopeXFirstClusterVsMom");
 213   gSigmaResSlopeXFirstClusterVsMom->SetTitle("#sigma_{slope_{X}} at first cluster versus p;p (GeV/c);#sigma_{slope_{X}}");
 214   TGraphAsymmErrors* gMeanResSlopeYFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
 215   gMeanResSlopeYFirstClusterVsMom->SetName("gMeanResSlopeYFirstClusterVsMom");
 216   gMeanResSlopeYFirstClusterVsMom->SetTitle("<#Delta_{slope_{Y}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
 217   TGraphAsymmErrors* gSigmaResSlopeYFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
 218   gSigmaResSlopeYFirstClusterVsMom->SetName("gSigmaResSlopeYFirstClusterVsMom");
 219   gSigmaResSlopeYFirstClusterVsMom->SetTitle("#sigma_{slope_{Y}} at first cluster versus p;p (GeV/c);#sigma_{slope_{Y}}");
 220
 221   // DCA resolution and MCS angular dispersion
 222   histoFile->mkdir("DCA","DCA");
 223   histoFile->cd("DCA");
 224
 225   const Int_t deltaPDCANBins = 500;
 226   const Double_t deltaPDCAEdges[2] = {0., 1000.};
 227   const Double_t deltaPMCSAngEdges[2] = {-0.5, 0.5};
 228
 229   TH1F *hPDCA = new TH1F("hPDCA","p #times DCA at vertex;p #times DCA (GeV #times cm)", deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 230   TH2D *hPDCAVsMom_2_3_Deg = new TH2D("hPDCAVsMom_2_3_Deg","p #times DCA versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 231   TH2D *hPDCAVsMom_3_10_Deg = new TH2D("hPDCAVsMom_3_10_Deg","p #times DCA versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 232   TH2D *hPMCSAngVsMom_2_3_Deg = new TH2D("hPMCSAngVsMom_2_3_Deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
 233   TH2D *hPMCSAngVsMom_3_10_Deg = new TH2D("hPMCSAngVsMom_3_10_Deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
 234
 235   TH2D *hPDCAVsPosAbsEnd_0_2_DegMC = new TH2D("hPDCAVsPosAbsEnd_0_2_DegMC","p #times DCA versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 236   TH2D *hPDCAVsPosAbsEnd_2_3_DegMC = new TH2D("hPDCAVsPosAbsEnd_2_3_DegMC","p #times DCA}versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 237   TH2D *hPDCAVsPosAbsEnd_3_10_DegMC = new TH2D("hPDCAVsPosAbsEnd_3_10_DegMC","p #times DCA versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 238
 239   TH2D *hPDCAVsAngle = new TH2D("hPDCAVsAngle","p #times DCA versus track position at absorber end converted to degrees;angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 240   TH2D *hPDCAVsMCAngle = new TH2D("hPDCAVsMCAngle","p #times DCA versus MC angle;MC angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
 241
 242   TGraphAsymmErrors* gSigmaPDCAVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
 243   gSigmaPDCAVsMom_2_3_Deg->SetName("gSigmaPDCAVsMom_2_3_Deg");
 244   gSigmaPDCAVsMom_2_3_Deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
 245   TGraphAsymmErrors* gSigmaPDCAVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
 246   gSigmaPDCAVsMom_3_10_Deg->SetName("gSigmaPDCAVsMom_3_10_Deg");
 247   gSigmaPDCAVsMom_3_10_Deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
 248   TGraphAsymmErrors* gMeanPMCSAngVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
 249   gMeanPMCSAngVsMom_2_3_Deg->SetName("gMeanPMCSAngVsMom_2_3_Deg");
 250   gMeanPMCSAngVsMom_2_3_Deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
 251   TGraphAsymmErrors* gSigmaPMCSAngVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
 252   gSigmaPMCSAngVsMom_2_3_Deg->SetName("gSigmaPMCSAngVsMom_2_3_Deg");
 253   gSigmaPMCSAngVsMom_2_3_Deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
 254   TGraphAsymmErrors* gMeanPMCSAngVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
 255   gMeanPMCSAngVsMom_3_10_Deg->SetName("gMeanPMCSAngVsMom_3_10_Deg");
 256   gMeanPMCSAngVsMom_3_10_Deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
 257   TGraphAsymmErrors* gSigmaPMCSAngVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
 258   gSigmaPMCSAngVsMom_3_10_Deg->SetName("gSigmaPMCSAngVsMom_3_10_Deg");
 259   gSigmaPMCSAngVsMom_3_10_Deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
 260
 261   // eta resolution at vertex
 262   histoFile->mkdir("etaAtVertex","etaAtVertex");
 263   histoFile->cd("etaAtVertex");
 264
 265   const Int_t deltaEtaAtVtxNBins = 500;
 266   const Double_t deltaEtaAtVtxEdges[2] = {-0.5, 0.5};
 267
 268   TH1F *hResEtaVertex = new TH1F("hResEtaVertex","#Delta_{eta} at vertex;#Delta_{eta}", deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 269   TH2D *hResEtaVertexVsMom = new TH2D("hResEtaVertexVsMom","#Delta_{eta} at vertex versus p;p (GeV/c);#Delta_{eta}",2*pNBins,pEdges[0],pEdges[1], deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 270
 271   TH2D *hResEtaVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResEtaVertexVsPosAbsEnd_0_2_DegMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 272   TH2D *hResEtaVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResEtaVertexVsPosAbsEnd_2_3_DegMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 273   TH2D *hResEtaVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResEtaVertexVsPosAbsEnd_3_10_DegMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 274
 275   TH2D *hResEtaVertexVsAngle = new TH2D("hResEtaVertexVsAngle","#Delta_{eta} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 276   TH2D *hResEtaVertexVsMCAngle = new TH2D("hResEtaVertexVsMCAngle","#Delta_{eta} at vertex versus MC angle;MC angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
 277
 278   TGraphAsymmErrors* gMeanResEtaVertexVsMom = new TGraphAsymmErrors(pNBins);
 279   gMeanResEtaVertexVsMom->SetName("gMeanResEtaVertexVsMom");
 280   gMeanResEtaVertexVsMom->SetTitle("<#Delta_{eta}> at vertex versus p;p (GeV/c);<#Delta_{eta}>");
 281   TGraphAsymmErrors* gSigmaResEtaVertexVsMom = new TGraphAsymmErrors(pNBins);
 282   gSigmaResEtaVertexVsMom->SetName("gSigmaResEtaVertexVsMom");
 283   gSigmaResEtaVertexVsMom->SetTitle("#sigma_{eta} at vertex versus p;p (GeV/c);#sigma_{eta}");
 284
 285   // phi resolution at vertex
 286   histoFile->mkdir("phiAtVertex","phiAtVertex");
 287   histoFile->cd("phiAtVertex");
 288
 289   const Int_t deltaPhiAtVtxNBins = 500;
 290   const Double_t deltaPhiAtVtxEdges[2] = {-0.5, 0.5};
 291
 292   TH1F *hResPhiVertex = new TH1F("hResPhiVertex","#Delta_{phi} at vertex;#Delta_{phi}", deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 293   TH2D *hResPhiVertexVsMom = new TH2D("hResPhiVertexVsMom","#Delta_{phi} at vertex versus p;p (GeV/c);#Delta_{phi}",2*pNBins,pEdges[0],pEdges[1], deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 294
 295   TH2D *hResPhiVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResPhiVertexVsPosAbsEnd_0_2_DegMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 296   TH2D *hResPhiVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResPhiVertexVsPosAbsEnd_2_3_DegMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 297   TH2D *hResPhiVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResPhiVertexVsPosAbsEnd_3_10_DegMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 298
 299   TH2D *hResPhiVertexVsAngle = new TH2D("hResPhiVertexVsAngle","#Delta_{phi} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 300   TH2D *hResPhiVertexVsMCAngle = new TH2D("hResPhiVertexVsMCAngle","#Delta_{phi} at vertex versus MC angle;MC angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
 301
 302   TGraphAsymmErrors* gMeanResPhiVertexVsMom = new TGraphAsymmErrors(pNBins);
 303   gMeanResPhiVertexVsMom->SetName("gMeanResPhiVertexVsMom");
 304   gMeanResPhiVertexVsMom->SetTitle("<#Delta_{phi}> at vertex versus p;p (GeV/c);<#Delta_{phi}>");
 305   TGraphAsymmErrors* gSigmaResPhiVertexVsMom = new TGraphAsymmErrors(pNBins);
 306   gSigmaResPhiVertexVsMom->SetName("gSigmaResPhiVertexVsMom");
 307   gSigmaResPhiVertexVsMom->SetTitle("#sigma_{phi} at vertex versus p;p (GeV/c);#sigma_{phi}");
 308
 309   // cluster resolution
 310   histoFile->mkdir("clusters","clusters");
 311   histoFile->cd("clusters");
 312
 313   TH1F* hResidualXInCh[AliMUONConstants::NTrackingCh()];
 314   TH1F* hResidualYInCh[AliMUONConstants::NTrackingCh()];
 315   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
 316     hResidualXInCh[i] = new TH1F(Form("hResidualXInCh%d",i+1), Form("cluster-track residual-X distribution in chamber %d;#Delta_{X} (cm)",i+1), 1000, -1., 1.);
 317     hResidualYInCh[i] = new TH1F(Form("hResidualYInCh%d",i+1), Form("cluster-track residual-Y distribution in chamber %d;#Delta_{Y} (cm)",i+1), 1000, -0.5, 0.5);
 318   }
 319
 320   TGraphErrors* gResidualXPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
 321   gResidualXPerChMean->SetName("gResidualXPerChMean");
 322   gResidualXPerChMean->SetTitle("cluster-trackRef residual-X per Ch: mean;chamber ID;<#Delta_{X}> (cm)");
 323   gResidualXPerChMean->SetMarkerStyle(kFullDotLarge);
 324   TGraphErrors* gResidualYPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
 325   gResidualYPerChMean->SetName("gResidualYPerChMean");
 326   gResidualYPerChMean->SetTitle("cluster-trackRef residual-Y per Ch: mean;chamber ID;<#Delta_{Y}> (cm)");
 327   gResidualYPerChMean->SetMarkerStyle(kFullDotLarge);
 328   TGraphErrors* gResidualXPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
 329   gResidualXPerChSigma->SetName("gResidualXPerChSigma");
 330   gResidualXPerChSigma->SetTitle("cluster-trackRef residual-X per Ch: sigma;chamber ID;#sigma_{X} (cm)");
 331   gResidualXPerChSigma->SetMarkerStyle(kFullDotLarge);
 332   TGraphErrors* gResidualYPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
 333   gResidualYPerChSigma->SetName("gResidualYPerChSigma");
 334   gResidualYPerChSigma->SetTitle("cluster-trackRef residual-Y per Ch: sigma;chamber ID;#sigma_{Y} (cm)");
 335   gResidualYPerChSigma->SetMarkerStyle(kFullDotLarge);
 336
 337   histoFile->mkdir("trigger");
 338   histoFile->cd("trigger");
 339   TH1F* hResidualTrigX11 = new TH1F("hResiudalTrigX11", "Residual X11", 100, -10., 10.);
 340   TH1F* hResidualTrigY11 = new TH1F("hResiudalTrigY11", "Residual Y11", 100, -10., 10.);
 341   TH1F* hResidualTrigSlopeY = new TH1F("hResiudalTrigSlopeY", "Residual Y slope", 100, -0.1, 0.1);
 342   TH1F* hTriggerableMatchFailed = new TH1F("hTriggerableMatchFailed", "Triggerable multiplicity for events with no match", 15, -0.5, 14.5);
 343
 344
 345   // ###################################### initialize ###################################### //
 346   AliMUONRecoCheck rc(esdFileName, pathSim);
 347
 348   // load necessary data from OCDB
 349   AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
 350   AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
 351   if (!AliMUONCDB::LoadField()) return;
 352   AliMUONTrackExtrap::SetField();
 353   AliGeomManager::LoadGeometry();
 354   if (!AliGeomManager::GetGeometry()) return;
 355   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
 356   if (!recoParam) return;
 357   AliMUONESDInterface::ResetTracker(recoParam);
 358
 359   // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
 360   Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
 361   // compute the mask of requested stations from recoParam
 362   UInt_t requestedStationMask = 0;
 363   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
 364   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
 365   Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
 366
 367   Int_t nevents = rc.NumberOfEvents();
 368
 369   if (nevents < nEvent || nEvent < 0) nEvent = nevents;
 370
 371   Int_t ievent;
 372   Int_t nReconstructibleTracks = 0;
 373   Int_t nReconstructedTracks = 0;
 374   Int_t nReconstructibleTracksCheck = 0;
 375   AliMUONTrackParam *trackParam;
 376   Double_t x1,y1,z1,slopex1,slopey1,pX1,pY1,pZ1,p1,pT1,eta1,phi1;
 377   Double_t x2,y2,z2,slopex2,slopey2,pX2,pY2,pZ2,p2,pT2,eta2,phi2;
 378   Double_t dPhi;
 379   Double_t xAbs,yAbs,dAbs,aAbs,aMCS,aMC;
 380   Double_t xDCA,yDCA,dca,pU;
 381   Double_t aMCSMoy = 0., aMCS2Moy = 0., dMCSMoy = 0., dMCS2Moy = 0., adMCSMoy = 0.;
 382   Int_t nMCS = 0;
 383
 384   // ###################################### fill histograms ###################################### //
 385   for (ievent=0; ievent<nEvent; ievent++)
 386   {
 387     if ((ievent+1)%100 == 0) cout<<"\rEvent processing... "<<ievent+1<<flush;
 388
 389     AliMUONVTrackStore* trackStore = rc.ReconstructedTracks(ievent, kFALSE);
 390     AliMUONVTrackStore* trackRefStore = rc.ReconstructibleTracks(ievent, requestedStationMask, request2ChInSameSt45);
 391
 392     hReconstructible->Fill(trackRefStore->GetSize());
 393     hReco->Fill(trackStore->GetSize());
 394
 395     nReconstructibleTracks += trackRefStore->GetSize();
 396     nReconstructedTracks += trackStore->GetSize();
 397
 398     AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(ievent);
 399     AliMUONVTriggerTrackStore* triggerTrackStore = rc.TriggeredTracks(ievent);
 400
 401     hTriggerable->Fill(triggerTrackRefStore->GetSize());
 402     hTriggered->Fill(triggerTrackStore->GetSize());
 403
 404     // loop over trigger trackRef
 405     TIter nextTrig(triggerTrackRefStore->CreateIterator());
 406     AliMUONTriggerTrack* triggerTrackRef;
 407     Int_t nTriggerMatches = 0;
 408     while ( ( triggerTrackRef = static_cast<AliMUONTriggerTrack*>(nextTrig()) ) )
 409     {
 410
 411       AliMUONTriggerTrack* triggerTrackMatched = 0x0;
 412
 413       // loop over trackReco and look for compatible track
 414       TIter nextTrig2(triggerTrackStore->CreateIterator());
 415       AliMUONTriggerTrack* triggerTrackReco;
 416       while ( ( triggerTrackReco = static_cast<AliMUONTriggerTrack*>(nextTrig2()) ) )
 417       {
 418
 419         // check if trackReco is compatible with trackRef
 420         if (triggerTrackReco->Match(*triggerTrackRef, sigmaCut)) {
 421           triggerTrackMatched = triggerTrackReco;
 422           nTriggerMatches++;
 423           break;
 424         }
 425       }
 426
 427       if (triggerTrackMatched) { // tracking requirements verified, track is found
 428         hResidualTrigX11->Fill( triggerTrackMatched->GetX11() - triggerTrackRef->GetX11() );
 429         hResidualTrigY11->Fill( triggerTrackMatched->GetY11() - triggerTrackRef->GetY11() );
 430         hResidualTrigSlopeY->Fill( triggerTrackMatched->GetSlopeY() - triggerTrackRef->GetSlopeY() );
 431       }
 432     } // loop on trigger track ref
 433
 434     if ( nTriggerMatches != triggerTrackStore->GetSize() )
 435       hTriggerableMatchFailed->Fill(triggerTrackRefStore->GetSize());
 436
 437     // loop over trackRef
 438     TIter next(trackRefStore->CreateIterator());
 439     AliMUONTrack* trackRef;
 440     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) )
 441     {
 442
 443       hTrackRefID->Fill(trackRef->GetUniqueID());
 444
 445       AliMUONTrack* trackMatched = 0x0;
 446       Int_t nMatchClusters = 0;
 447
 448       // loop over trackReco and look for compatible track
 449       TIter next2(trackStore->CreateIterator());
 450       AliMUONTrack* trackReco;
 451       while ( ( trackReco = static_cast<AliMUONTrack*>(next2()) ) )
 452       {
 453
 454         // check if trackReco is compatible with trackRef
 455         if (trackReco->Match(*trackRef, sigmaCut, nMatchClusters)) {
 456           trackMatched = trackReco;
 457           break;
 458         }
 459
 460       }
 461
 462       if (trackMatched) { // tracking requirements verified, track is found
 463         nReconstructibleTracksCheck++;
 464         hNClusterComp->Fill(nMatchClusters);
 465
 466         // compute track position at the end of the absorber
 467         AliMUONTrackParam trackParamAtAbsEnd(*((AliMUONTrackParam*)trackMatched->GetTrackParamAtCluster()->First()));
 468         AliMUONTrackExtrap::ExtrapToZ(&trackParamAtAbsEnd, AliMUONConstants::AbsZEnd());
 469         xAbs = trackParamAtAbsEnd.GetNonBendingCoor();
 470         yAbs = trackParamAtAbsEnd.GetBendingCoor();
 471         dAbs = TMath::Sqrt(xAbs*xAbs + yAbs*yAbs);
 472         aAbs = TMath::ATan(-dAbs/AliMUONConstants::AbsZEnd()) * TMath::RadToDeg();
 473         pX2 = trackParamAtAbsEnd.Px();
 474         pY2 = trackParamAtAbsEnd.Py();
 475         pZ2 = trackParamAtAbsEnd.Pz();
 476         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
 477         aMCS = TMath::ATan(-pT2/pZ2) * TMath::RadToDeg();
 478
 479         trackParam = trackRef->GetTrackParamAtVertex();
 480         x1 = trackParam->GetNonBendingCoor();
 481         y1 = trackParam->GetBendingCoor();
 482         z1 = trackParam->GetZ();
 483         slopex1 = trackParam->GetNonBendingSlope();
 484         slopey1 = trackParam->GetBendingSlope();
 485         pX1 = trackParam->Px();
 486         pY1 = trackParam->Py();
 487         pZ1 = trackParam->Pz();
 488         p1  = trackParam->P();
 489         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
 490         aMC = TMath::ATan(-pT1/pZ1) * TMath::RadToDeg();
 491         eta1 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT1/pZ1)));
 492         phi1 = TMath::Pi()+TMath::ATan2(-pY1, -pX1);
 493
 494         trackParam = trackMatched->GetTrackParamAtVertex();
 495         x2 = trackParam->GetNonBendingCoor();
 496         y2 = trackParam->GetBendingCoor();
 497         z2 = trackParam->GetZ();
 498         slopex2 = trackParam->GetNonBendingSlope();
 499         slopey2 = trackParam->GetBendingSlope();
 500         pX2 = trackParam->Px();
 501         pY2 = trackParam->Py();
 502         pZ2 = trackParam->Pz();
 503         p2  = trackParam->P();
 504         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
 505         eta2 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT2/pZ2)));
 506         phi2 = TMath::Pi()+TMath::ATan2(-pY2, -pX2);
 507
 508         dPhi = phi2-phi1;
 509         if (dPhi < -TMath::Pi()) dPhi += 2.*TMath::Pi();
 510         else if (dPhi > TMath::Pi()) dPhi -= 2.*TMath::Pi();
 511
 512         AliMUONTrackParam trackParamAtDCA(*((AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First()));
 513         pU = trackParamAtDCA.P();
 514         AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&trackParamAtDCA, z2);
 515         xDCA = trackParamAtDCA.GetNonBendingCoor();
 516         yDCA = trackParamAtDCA.GetBendingCoor();
 517         dca = TMath::Sqrt(xDCA*xDCA + yDCA*yDCA);
 518
 519         hResMomVertex->Fill(p2-p1);
 520         hResSlopeXVertex->Fill(slopex2-slopex1);
 521         hResSlopeYVertex->Fill(slopey2-slopey1);
 522         hPDCA->Fill(0.5*(p2+pU)*dca);
 523         hResEtaVertex->Fill(eta2-eta1);
 524         hResPhiVertex->Fill(dPhi);
 525         if (aMC >= aAbsLimits[0] && aMC <= aAbsLimits[1]) {
 526           hResMomVertexVsMom->Fill(p1,p2-p1);
 527           hResSlopeXVertexVsMom->Fill(p1,slopex2-slopex1);
 528           hResSlopeYVertexVsMom->Fill(p1,slopey2-slopey1);
 529           hResEtaVertexVsMom->Fill(p1,eta2-eta1);
 530           hResPhiVertexVsMom->Fill(p1,dPhi);
 531         }
 532         hResMomVertexVsAngleVsMom->Fill(p1,aAbs,p2-p1);
 533         if (aAbs > 2. && aAbs < 3.) {
 534           hResMomVertexVsMom_2_3_Deg->Fill(p1,p2-p1);
 535           hPDCAVsMom_2_3_Deg->Fill(p1,0.5*(p2+pU)*dca);
 536           hPMCSAngVsMom_2_3_Deg->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
 537         }
 538         else if (aAbs >= 3. && aAbs < 10.) {
 539           hResMomVertexVsMom_3_10_Deg->Fill(p1,p2-p1);
 540           hPDCAVsMom_3_10_Deg->Fill(p1,0.5*(p2+pU)*dca);
 541           hPMCSAngVsMom_3_10_Deg->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
 542           aMCSMoy += 0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad();
 543           aMCS2Moy += (0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad()) * (0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
 544           dMCSMoy += 0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd());
 545           dMCS2Moy += (0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd())) * (0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd()));
 546           adMCSMoy += (0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad()) * (0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd()));
 547           nMCS++;
 548         }
 549         if (aMC < 2.) {
 550           hResMomVertexVsMom_0_2_DegMC->Fill(p1,p2-p1);
 551           hResMomVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,p2-p1);
 552           hResSlopeXVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,slopex2-slopex1);
 553           hResSlopeYVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,slopey2-slopey1);
 554           hPDCAVsPosAbsEnd_0_2_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
 555           hResEtaVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,eta2-eta1);
 556           hResPhiVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,dPhi);
 557         }
 558         else if (aMC >= 2. && aMC < 3) {
 559           hResMomVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,p2-p1);
 560           hResSlopeXVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,slopex2-slopex1);
 561           hResSlopeYVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,slopey2-slopey1);
 562           hPDCAVsPosAbsEnd_2_3_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
 563           hResEtaVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,eta2-eta1);
 564           hResPhiVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,dPhi);
 565         }
 566         else if (aMC >= 3. && aMC < 10.) {
 567           hResMomVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,p2-p1);
 568           hResSlopeXVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,slopex2-slopex1);
 569           hResSlopeYVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,slopey2-slopey1);
 570           hPDCAVsPosAbsEnd_3_10_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
 571           hResEtaVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,eta2-eta1);
 572           hResPhiVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,dPhi);
 573         }
 574         hResMomVertexVsAngle->Fill(aAbs,p2-p1);
 575         hResSlopeXVertexVsAngle->Fill(aAbs,slopex2-slopex1);
 576         hResSlopeYVertexVsAngle->Fill(aAbs,slopey2-slopey1);
 577         hPDCAVsAngle->Fill(aAbs,0.5*(p2+pU)*dca);
 578         hResEtaVertexVsAngle->Fill(aAbs,eta2-eta1);
 579         hResPhiVertexVsAngle->Fill(aAbs,dPhi);
 580         hResMomVertexVsMCAngle->Fill(aMC,p2-p1);
 581         hResSlopeXVertexVsMCAngle->Fill(aMC,slopex2-slopex1);
 582         hResSlopeYVertexVsMCAngle->Fill(aMC,slopey2-slopey1);
 583         hPDCAVsMCAngle->Fill(aMC,0.5*(p2+pU)*dca);
 584         hResEtaVertexVsMCAngle->Fill(aMC,eta2-eta1);
 585         hResPhiVertexVsMCAngle->Fill(aMC,dPhi);
 586
 587         trackParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
 588         x1 = trackParam->GetNonBendingCoor();
 589         y1 = trackParam->GetBendingCoor();
 590         z1 = trackParam->GetZ();
 591         slopex1 = trackParam->GetNonBendingSlope();
 592         slopey1 = trackParam->GetBendingSlope();
 593         pX1 = trackParam->Px();
 594         pY1 = trackParam->Py();
 595         pZ1 = trackParam->Pz();
 596         p1  = trackParam->P();
 597         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
 598
 599         trackParam = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
 600         x2 = trackParam->GetNonBendingCoor();
 601         y2 = trackParam->GetBendingCoor();
 602         z2 = trackParam->GetZ();
 603         slopex2 = trackParam->GetNonBendingSlope();
 604         slopey2 = trackParam->GetBendingSlope();
 605         pX2 = trackParam->Px();
 606         pY2 = trackParam->Py();
 607         pZ2 = trackParam->Pz();
 608         p2  = trackParam->P();
 609         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
 610
 611         hResMomFirstCluster->Fill(p2-p1);
 612         hResMomFirstClusterVsMom->Fill(p1,p2-p1);
 613
 614         hResSlopeXFirstCluster->Fill(slopex2-slopex1);
 615         hResSlopeYFirstCluster->Fill(slopey2-slopey1);
 616         hResSlopeXFirstClusterVsMom->Fill(p1,slopex2-slopex1);
 617         hResSlopeYFirstClusterVsMom->Fill(p1,slopey2-slopey1);
 618
 619         // Fill residuals
 620         // Loop over clusters of first track
 621         AliMUONTrackParam* trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
 622         while (trackParamAtCluster1) {
 623           AliMUONVCluster* cluster1 = trackParamAtCluster1->GetClusterPtr();
 624           AliMUONTrackParam* trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
 625           while (trackParamAtCluster2) {
 626             AliMUONVCluster* cluster2 = trackParamAtCluster2->GetClusterPtr();
 627             if (cluster1->GetDetElemId() == cluster2->GetDetElemId()) {
 628               hResidualXInCh[cluster1->GetChamberId()]->Fill(cluster1->GetX() - cluster2->GetX());
 629               hResidualYInCh[cluster1->GetChamberId()]->Fill(cluster1->GetY() - cluster2->GetY());
 630               break;
 631             }
 632             trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->After(trackParamAtCluster2);
 633           }
 634           trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->After(trackParamAtCluster1);
 635         }
 636
 637       }
 638
 639     } // end loop track ref.
 640
 641   } // end loop on event
 642   cout<<"\rEvent processing... "<<nevents<<" done"<<endl;
 643
 644   // ###################################### compute stuff ###################################### //
 645   cout<<"\nWhen not specified, resolution at vertex is computed for ";
 646   if (absorberRegion == 1) cout<<"tracks in the absorber region [2,3] deg."<<endl;
 647   else if (absorberRegion == 2) cout<<"tracks in the absorber region [3,10] deg."<<endl;
 648   else cout<<"all tracks"<<endl;
 649
 650   // compute momentum resolution at vertex versus p
 651   TF1 *f2 = new TF1("f2",langaufun,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1],4);
 652   Int_t rebinFactorX = TMath::Max(hResMomVertexVsMom->GetNbinsX()/pNBins, 1);
 653   for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom->GetNbinsX(); i+=rebinFactorX) {
 654     cout<<"\rFitting momentum residuals at vertex... "<<i/rebinFactorX<<"/"<<pNBins<<flush;
 655     TH1D *tmp = hResMomVertexVsMom->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
 656     f2->SetParameters(0.2,0.,(Double_t)tmp->GetEntries(),1.);
 657     tmp->Fit("f2","WWNQ");
 658     Double_t fwhm = f2->GetParameter(0);
 659     Double_t sigma = f2->GetParameter(3);
 660     Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
 661     Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/tmp->GetBinWidth(1)),1);
 662     while (deltaPAtVtxNBins%rebin!=0) rebin--;
 663     tmp->Rebin(rebin);
 664     tmp->Fit("f2","NQ");
 665     fwhm = f2->GetParameter(0);
 666     sigma = f2->GetParameter(3);
 667     sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
 668     Double_t fwhmErr = f2->GetParError(0);
 669     Double_t sigmaErr = f2->GetParError(3);
 670     Double_t sigmaPErr = TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + fwhm*fwhm*fwhmErr*fwhmErr/(64.*log(2.)*log(2.))) / sigmaP;
 671     hResMomVertexVsMom->GetXaxis()->SetRange(i-rebinFactorX+1,i);
 672     Double_t p = (tmp->GetEntries() > 0) ? hResMomVertexVsMom->GetMean() : 0.5 * (hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom->GetBinLowEdge(i+1));
 673     hResMomVertexVsMom->GetXaxis()->SetRange();
 674     Double_t pErr[2] = {p-hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1), hResMomVertexVsMom->GetBinLowEdge(i+1)-p};
 675     gMeanResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, tmp->GetMean());
 676     gMeanResMomVertexVsMom->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], tmp->GetMeanError(), tmp->GetMeanError());
 677     gMostProbResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, -f2->GetParameter(1));
 678     gMostProbResMomVertexVsMom->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], f2->GetParError(1), f2->GetParError(1));
 679     gSigmaResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, 100.*sigmaP/p);
 680     gSigmaResMomVertexVsMom->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], 100.*sigmaPErr/p, 100.*sigmaPErr/p);
 681     delete tmp;
 682   }
 683   cout<<"\rFitting momentum residuals at vertex... "<<pNBins<<"/"<<pNBins<<endl;
 684
 685   // compute momentum relative resolution at first cluster versus p
 686   FitGausResVsMom(hResMomFirstClusterVsMom, pNBins, 0., 1., "momentum residuals at first cluster", gMeanResMomFirstClusterVsMom, gSigmaResMomFirstClusterVsMom);
 687   rebinFactorX = TMath::Max(hResMomFirstClusterVsMom->GetNbinsX()/pNBins, 1);
 688   for (Int_t i = rebinFactorX; i <= hResMomFirstClusterVsMom->GetNbinsX(); i+=rebinFactorX) {
 689     Double_t x,y;
 690     gSigmaResMomFirstClusterVsMom->GetPoint(i/rebinFactorX-1, x, y);
 691     gSigmaResMomFirstClusterVsMom->SetPoint(i/rebinFactorX-1, x, 100.*y/x);
 692     gSigmaResMomFirstClusterVsMom->SetPointEYlow(i/rebinFactorX-1, 100.*gSigmaResMomFirstClusterVsMom->GetErrorYlow(i/rebinFactorX-1)/x);
 693     gSigmaResMomFirstClusterVsMom->SetPointEYhigh(i/rebinFactorX-1, 100.*gSigmaResMomFirstClusterVsMom->GetErrorYhigh(i/rebinFactorX-1)/x);
 694   }
 695
 696   // compute slopeX resolution at vertex versus p
 697   FitGausResVsMom(hResSlopeXVertexVsMom, pNBins, 0., 2.e-3, "slopeX residuals at vertex", gMeanResSlopeXVertexVsMom, gSigmaResSlopeXVertexVsMom);
 698
 699   // compute slopeY resolution at vertex versus p
 700   FitGausResVsMom(hResSlopeYVertexVsMom, pNBins, 0., 2.e-3, "slopeY residuals at vertex", gMeanResSlopeYVertexVsMom, gSigmaResSlopeYVertexVsMom);
 701
 702   // compute slopeX resolution at first cluster versus p
 703   FitGausResVsMom(hResSlopeXFirstClusterVsMom, pNBins, 0., 3.e-4, "slopeX residuals at first cluster", gMeanResSlopeXFirstClusterVsMom, gSigmaResSlopeXFirstClusterVsMom);
 704
 705   // compute slopeY resolution at first cluster versus p
 706   FitGausResVsMom(hResSlopeYFirstClusterVsMom, pNBins, 0., 2.e-4, "slopeY residuals at first cluster", gMeanResSlopeYFirstClusterVsMom, gSigmaResSlopeYFirstClusterVsMom);
 707
 708   // compute p*DCA resolution in the region [2,3] deg at absorber end
 709   FitPDCAVsMom(hPDCAVsMom_2_3_Deg, pNBins, "p*DCA (tracks in [2,3] deg.)", gSigmaPDCAVsMom_2_3_Deg);
 710
 711   // compute p*DCA resolution in the region [3,10] deg at absorber end
 712   FitPDCAVsMom(hPDCAVsMom_3_10_Deg, pNBins, "p*DCA (tracks in [3,10] deg.)", gSigmaPDCAVsMom_3_10_Deg);
 713
 714   // compute MCS angular dispersion in the region [2,3] deg at absorber end
 715   FitGausResVsMom(hPMCSAngVsMom_2_3_Deg, pNBins, 0., 2.e-3, "p*MCSAngle (tracks in [2,3] deg.)", gMeanPMCSAngVsMom_2_3_Deg, gSigmaPMCSAngVsMom_2_3_Deg);
 716
 717   // compute MCS angular dispersion in the region [3,10] deg at absorber end
 718   FitGausResVsMom(hPMCSAngVsMom_3_10_Deg, pNBins, 0., 2.e-3, "p*MCSAngle (tracks in [3,10] deg.)", gMeanPMCSAngVsMom_3_10_Deg, gSigmaPMCSAngVsMom_3_10_Deg);
 719
 720   // compute eta resolution at vertex versus p
 721   FitGausResVsMom(hResEtaVertexVsMom, pNBins, 0., 0.1, "eta residuals at vertex", gMeanResEtaVertexVsMom, gSigmaResEtaVertexVsMom);
 722
 723   // compute phi resolution at vertex versus p
 724   FitGausResVsMom(hResPhiVertexVsMom, pNBins, 0., 0.01, "phi residuals at vertex", gMeanResPhiVertexVsMom, gSigmaResPhiVertexVsMom);
 725
 726   // compute cluster-track residual mean and dispersion
 727   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
 728     hResidualXInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualXInCh[i]->GetRMS(), 3.*hResidualXInCh[i]->GetRMS());
 729     gResidualXPerChMean->SetPoint(i, i+1, hResidualXInCh[i]->GetMean());
 730     gResidualXPerChMean->SetPointError(i, 0., hResidualXInCh[i]->GetMeanError());
 731     gResidualXPerChSigma->SetPoint(i, i+1, hResidualXInCh[i]->GetRMS());
 732     gResidualXPerChSigma->SetPointError(i, 0., hResidualXInCh[i]->GetRMSError());
 733     hResidualXInCh[i]->GetXaxis()->SetRange(0,0);
 734     hResidualYInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualYInCh[i]->GetRMS(), 3.*hResidualYInCh[i]->GetRMS());
 735     gResidualYPerChMean->SetPoint(i, i+1, hResidualYInCh[i]->GetMean());
 736     gResidualYPerChMean->SetPointError(i, 0., hResidualYInCh[i]->GetMeanError());
 737     gResidualYPerChSigma->SetPoint(i, i+1, hResidualYInCh[i]->GetRMS());
 738     gResidualYPerChSigma->SetPointError(i, 0., hResidualYInCh[i]->GetRMSError());
 739     hResidualYInCh[i]->GetXaxis()->SetRange(0,0);
 740   }
 741
 742   // ###################################### display histograms ###################################### //
 743   // diplay momentum residuals
 744   TCanvas* cResMom = DrawVsAng("cResMom", "momentum residual at vertex in 3 angular regions", hResMomVertex, hResMomVertexVsAngle);
 745   TCanvas* cResMomMC = DrawVsAng("cResMomMC", "momentum residual at vertex in 3 MC angular regions", hResMomVertex, hResMomVertexVsMCAngle);
 746   TCanvas* cResMomVsPos = DrawVsPos("cResMomVsPos", "momentum residual at vertex versus position at absorber end in 3 MC angular regions",
 747                                     hResMomVertexVsPosAbsEnd_0_2_DegMC, hResMomVertexVsPosAbsEnd_2_3_DegMC, hResMomVertexVsPosAbsEnd_3_10_DegMC);
 748   TCanvas* cResMom_2_3_Deg = DrawResMomVsMom("cResMom_2_3_Deg", "momentum residual for tracks between 2 and 3 degrees",
 749                                              hResMomVertexVsMom_2_3_Deg, 10, f2, "momentum residuals at vertex (tracks in [2,3] deg.)");
 750   TCanvas* cResMom_3_10_Deg = DrawResMomVsMom("cResMom_3_10_Deg", "momentum residual for tracks between 3 and 10 degrees",
 751                                               hResMomVertexVsMom_3_10_Deg, 10, f2, "momentum residuals at vertex (tracks in [3,10] deg.)");
 752   TCanvas* cResMom_0_2_DegMC = DrawResMomVsMom("cResMom_0_2_DegMC", "momentum residuals for tracks with MC angle < 2 degrees", hResMomVertexVsMom_0_2_DegMC, 5);
 753
 754   // diplay slopeX residuals
 755   TCanvas* cResSlopeX = DrawVsAng("cResSlopeX", "slope_{X} residual at vertex in 3 angular regions", hResSlopeXVertex, hResSlopeXVertexVsAngle);
 756   TCanvas* cResSlopeXMC = DrawVsAng("cResSlopeXMC", "slope_{X} residual at vertex in 3 MC angular regions", hResSlopeXVertex, hResSlopeXVertexVsMCAngle);
 757   TCanvas* cResSlopeXVsPos = DrawVsPos("cResSlopeXVsPos", "slope_{X} residual at vertex versus position at absorber end in 3 MC angular regions",
 758                                        hResSlopeXVertexVsPosAbsEnd_0_2_DegMC, hResSlopeXVertexVsPosAbsEnd_2_3_DegMC, hResSlopeXVertexVsPosAbsEnd_3_10_DegMC);
 759
 760   // diplay slopeY residuals
 761   TCanvas* cResSlopeY = DrawVsAng("cResSlopeY", "slope_{Y} residual at vertex in 3 angular regions", hResSlopeYVertex, hResSlopeYVertexVsAngle);
 762   TCanvas* cResSlopeYMC = DrawVsAng("cResSlopeYMC", "slope_{Y} residual at vertex in 3 MC angular regions", hResSlopeYVertex, hResSlopeYVertexVsMCAngle);
 763   TCanvas* cResSlopeYVsPos = DrawVsPos("cResSlopeYVsPos", "slope_{Y} residual at vertex versus position at absorber end in 3 MC angular regions",
 764                                        hResSlopeYVertexVsPosAbsEnd_0_2_DegMC, hResSlopeYVertexVsPosAbsEnd_2_3_DegMC, hResSlopeYVertexVsPosAbsEnd_3_10_DegMC);
 765
 766   // diplay P*DCA
 767   TCanvas* cPDCA = DrawVsAng("cPDCA", "p #times DCA in 3 angular regions", hPDCA, hPDCAVsAngle);
 768   TCanvas* cPDCAMC = DrawVsAng("cPDCAMC", "p #times DCA in 3 MC angular regions", hPDCA, hPDCAVsMCAngle);
 769   TCanvas* cPDCAVsPos = DrawVsPos("cPDCAVsPos", "p #times DCA versus position at absorber end in 3 MC angular regions",
 770                                   hPDCAVsPosAbsEnd_0_2_DegMC, hPDCAVsPosAbsEnd_2_3_DegMC, hPDCAVsPosAbsEnd_3_10_DegMC);
 771
 772   // diplay eta residuals
 773   TCanvas* cResEta = DrawVsAng("cResEta", "eta residual at vertex in 3 angular regions", hResEtaVertex, hResEtaVertexVsAngle);
 774   TCanvas* cResEtaMC = DrawVsAng("cResEtaMC", "eta residual at vertex in 3 MC angular regions", hResEtaVertex, hResEtaVertexVsMCAngle);
 775   TCanvas* cResEtaVsPos = DrawVsPos("cResEtaVsPos", "eta residual at vertex versus position at absorber end in 3 MC angular regions",
 776                                     hResEtaVertexVsPosAbsEnd_0_2_DegMC, hResEtaVertexVsPosAbsEnd_2_3_DegMC, hResEtaVertexVsPosAbsEnd_3_10_DegMC);
 777
 778   // diplay phi residuals
 779   TCanvas* cResPhi = DrawVsAng("cResPhi", "phi residual at vertex in 3 angular regions", hResPhiVertex, hResPhiVertexVsAngle);
 780   TCanvas* cResPhiMC = DrawVsAng("cResPhiMC", "phi residual at vertex in 3 MC angular regions", hResPhiVertex, hResPhiVertexVsMCAngle);
 781   TCanvas* cResPhiVsPos = DrawVsPos("cResPhiVsPos", "phi residual at vertex versus position at absorber end in 3 MC angular regions",
 782                                     hResPhiVertexVsPosAbsEnd_0_2_DegMC, hResPhiVertexVsPosAbsEnd_2_3_DegMC, hResPhiVertexVsPosAbsEnd_3_10_DegMC);
 783
 784   // ###################################### save histogram ###################################### //
 785   histoFile->Write();
 786
 787   histoFile->cd("momentumAtVertex");
 788   gMeanResMomVertexVsMom->Write();
 789   gMostProbResMomVertexVsMom->Write();
 790   gSigmaResMomVertexVsMom->Write();
 791   cResMom->Write();
 792   cResMomMC->Write();
 793   cResMomVsPos->Write();
 794   cResMom_2_3_Deg->Write();
 795   cResMom_3_10_Deg->Write();
 796   cResMom_0_2_DegMC->Write();
 797
 798   histoFile->cd("slopesAtVertex");
 799   gMeanResSlopeXVertexVsMom->Write();
 800   gMeanResSlopeYVertexVsMom->Write();
 801   gSigmaResSlopeXVertexVsMom->Write();
 802   gSigmaResSlopeYVertexVsMom->Write();
 803   cResSlopeX->Write();
 804   cResSlopeY->Write();
 805   cResSlopeXMC->Write();
 806   cResSlopeYMC->Write();
 807   cResSlopeXVsPos->Write();
 808   cResSlopeYVsPos->Write();
 809
 810   histoFile->cd("DCA");
 811   gSigmaPDCAVsMom_2_3_Deg->Write();
 812   gSigmaPDCAVsMom_3_10_Deg->Write();
 813   gMeanPMCSAngVsMom_2_3_Deg->Write();
 814   gSigmaPMCSAngVsMom_2_3_Deg->Write();
 815   gMeanPMCSAngVsMom_3_10_Deg->Write();
 816   gSigmaPMCSAngVsMom_3_10_Deg->Write();
 817   cPDCA->Write();
 818   cPDCAMC->Write();
 819   cPDCAVsPos->Write();
 820
 821   histoFile->cd("etaAtVertex");
 822   gMeanResEtaVertexVsMom->Write();
 823   gSigmaResEtaVertexVsMom->Write();
 824   cResEta->Write();
 825   cResEtaMC->Write();
 826   cResEtaVsPos->Write();
 827
 828   histoFile->cd("phiAtVertex");
 829   gMeanResPhiVertexVsMom->Write();
 830   gSigmaResPhiVertexVsMom->Write();
 831   cResPhi->Write();
 832   cResPhiMC->Write();
 833   cResPhiVsPos->Write();
 834
 835   histoFile->cd("momentumAtFirstCluster");
 836   gMeanResMomFirstClusterVsMom->Write();
 837   gSigmaResMomFirstClusterVsMom->Write();
 838
 839   histoFile->cd("slopesAtFirstCluster");
 840   gMeanResSlopeXFirstClusterVsMom->Write();
 841   gMeanResSlopeYFirstClusterVsMom->Write();
 842   gSigmaResSlopeXFirstClusterVsMom->Write();
 843   gSigmaResSlopeYFirstClusterVsMom->Write();
 844
 845   histoFile->cd("clusters");
 846   gResidualXPerChMean->Write();
 847   gResidualXPerChSigma->Write();
 848   gResidualYPerChMean->Write();
 849   gResidualYPerChSigma->Write();
 850
 851   histoFile->Close();
 852
 853   // ###################################### clean memory ###################################### //
 854   delete cResMom;
 855   delete cResMomMC;
 856   delete cResMomVsPos;
 857   delete cResMom_2_3_Deg;
 858   delete cResMom_3_10_Deg;
 859   delete cResMom_0_2_DegMC;
 860   delete cResSlopeX;
 861   delete cResSlopeY;
 862   delete cResSlopeXMC;
 863   delete cResSlopeYMC;
 864   delete cResSlopeXVsPos;
 865   delete cResSlopeYVsPos;
 866   delete cPDCA;
 867   delete cPDCAMC;
 868   delete cPDCAVsPos;
 869   delete cResEta;
 870   delete cResEtaMC;
 871   delete cResEtaVsPos;
 872   delete cResPhi;
 873   delete cResPhiMC;
 874   delete cResPhiVsPos;
 875
 876   // ###################################### print statistics ###################################### //
 877   printf("\n");
 878   printf("nb of reconstructible tracks: %d \n", nReconstructibleTracks);
 879   printf("nb of reconstructed tracks: %d \n", nReconstructedTracks);
 880   printf("nb of reconstructible tracks which are reconstructed: %d \n", nReconstructibleTracksCheck);
 881
 882   aMCSMoy /= (Double_t) nMCS;
 883   aMCS2Moy /= (Double_t) nMCS;
 884   dMCSMoy /= (Double_t) nMCS;
 885   dMCS2Moy /= (Double_t) nMCS;
 886   adMCSMoy /= (Double_t) nMCS;
 887   Double_t sigma2_ThetaMCS = aMCS2Moy - aMCSMoy*aMCSMoy;
 888   Double_t sigma2_PosMCS = dMCS2Moy - dMCSMoy*dMCSMoy;
 889   Double_t cov_ThetaPosMCS = - (adMCSMoy - aMCSMoy*dMCSMoy);
 890   printf("\nmultiple scattering of tracks between 3 and 10 deg. at absorber end:\n");
 891   printf(" sigma_ThetaMCS = %f\n", TMath::Sqrt(sigma2_ThetaMCS));
 892   printf(" sigma_PosMCS = %f\n", TMath::Sqrt(sigma2_PosMCS));
 893   printf(" cov_ThetaPosMCS = %f\n", cov_ThetaPosMCS);
 894   printf(" --> sigma_DCA = %f\n", TMath::Sqrt(AliMUONConstants::AbsZEnd()*AliMUONConstants::AbsZEnd()*sigma2_ThetaMCS
 895                                               - 2.*AliMUONConstants::AbsZEnd()*cov_ThetaPosMCS + sigma2_PosMCS));
 896   printf("\n");
 897 }
 898
 899 //------------------------------------------------------------------------------------
 900 Double_t langaufun(Double_t *x, Double_t *par) {
 901
 902   //Fit parameters:
 903   //par[0]=Width (scale) parameter of Landau density
 904   //par[1]=Most Probable (MP, location) parameter of Landau density
 905   //par[2]=Total area (integral -inf to inf, normalization constant)
 906   //par[3]=Width (sigma) of convoluted Gaussian function
 907   //
 908   //In the Landau distribution (represented by the CERNLIB approximation),
 909   //the maximum is located at x=-0.22278298 with the location parameter=0.
 910   //This shift is corrected within this function, so that the actual
 911   //maximum is identical to the MP parameter.
 912
 913   // Numeric constants
 914   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 915   Double_t mpshift  = -0.22278298;       // Landau maximum location
 916
 917   // Control constants
 918   Double_t np = 100.0; // number of convolution steps
 919   Double_t sc = 5.0;   // convolution extends to +-sc Gaussian sigmas
 920
 921   // Variables
 922   Double_t xx;
 923   Double_t mpc;
 924   Double_t fland;
 925   Double_t sum = 0.0;
 926   Double_t xlow,xupp;
 927   Double_t step;
 928   Double_t i;
 929
 930
 931   // MP shift correction
 932   mpc = par[1] - mpshift * par[0];
 933
 934   // Range of convolution integral
 935   xlow = x[0] - sc * par[3];
 936   xupp = x[0] + sc * par[3];
 937
 938   step = (xupp-xlow) / np;
 939
 940   // Convolution integral of Landau and Gaussian by sum
 941   for(i=1.0; i<=np/2; i++) {
 942     xx = xlow + (i-.5) * step;
 943     //change x -> -x because the tail of the Landau is at the left here...
 944     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
 945     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 946
 947     //change x -> -x because the tail of the Landau is at the left here...
 948     xx = xupp - (i-.5) * step;
 949     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
 950     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 951   }
 952
 953   return (par[2] * step * sum * invsq2pi / par[3]);
 954 }
 955
 956 //------------------------------------------------------------------------------------
 957 void FitGausResVsMom(TH2* h, Int_t nBins, const Double_t mean0, const Double_t sigma0,
 958                      const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
 959 {
 960   /// generic function to fit residuals versus momentum with a gaussian
 961   static TF1* fGaus = 0x0;
 962   if (!fGaus) fGaus = new TF1("fGaus","gaus");
 963
 964   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
 965   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
 966     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
 967     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
 968     fGaus->SetParameters(tmp->GetEntries(), mean0, sigma0);
 969     tmp->Fit("fGaus","WWNQ");
 970     Int_t rebin = TMath::Max(Int_t(0.5*fGaus->GetParameter(2)/tmp->GetBinWidth(1)),1);
 971     while (tmp->GetNbinsX()%rebin!=0) rebin--;
 972     tmp->Rebin(rebin);
 973     tmp->Fit("fGaus","NQ");
 974     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
 975     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
 976     h->GetXaxis()->SetRange();
 977     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
 978     gMean->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(1));
 979     gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(1), fGaus->GetParError(1));
 980     gSigma->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(2));
 981     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(2), fGaus->GetParError(2));
 982     delete tmp;
 983   }
 984   cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
 985 }
 986
 987 //------------------------------------------------------------------------------------
 988 void FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gSigma)
 989 {
 990   /// generic function to fit p*DCA distributions
 991   static TF1* fPGaus = 0x0;
 992   if (!fPGaus) fPGaus = new TF1("fPGaus","x*gaus");
 993
 994   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
 995   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
 996     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
 997     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
 998     fPGaus->SetParameters(1.,-100.,100.);
 999     Int_t rebin = 50.*(tmp->GetNbinsX()/(tmp->GetBinLowEdge(tmp->GetNbinsX()+1)-tmp->GetBinLowEdge(1)));
1000     while (tmp->GetNbinsX()%rebin!=0) rebin--;
1001     tmp->Rebin(rebin);
1002     tmp->Fit("fPGaus","NQ");
1003     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
1004     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1005     h->GetXaxis()->SetRange();
1006     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
1007     gSigma->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(2));
1008     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(2), fPGaus->GetParError(2));
1009     delete tmp;
1010   }
1011   cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
1012 }
1013
1014 //------------------------------------------------------------------------------------
1015 TCanvas* DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2)
1016 {
1017   /// generic function to draw histograms versus absorber angular region
1018   TCanvas* c = new TCanvas(name, title);
1019   c->cd();
1020   h1->Draw();
1021   TH1D *proj1 = h2->ProjectionY(Form("%s_proj_0_2",h2->GetName()),1,2);
1022   proj1->Draw("sames");
1023   proj1->SetLineColor(2);
1024   TH1D *proj2 = h2->ProjectionY(Form("%s_proj_2_3",h2->GetName()),3,3);
1025   proj2->Draw("sames");
1026   proj2->SetLineColor(4);
1027   TH1D *proj3 = h2->ProjectionY(Form("%s__proj_3_10",h2->GetName()),4,10);
1028   proj3->Draw("sames");
1029   proj3->SetLineColor(3);
1030   return c;
1031 }
1032
1033 //------------------------------------------------------------------------------------
1034 TCanvas* DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3)
1035 {
1036   /// generic function to draw histograms versus position at absorber end
1037   TCanvas* c = new TCanvas(name, title);
1038   c->cd();
1039   h1->Draw();
1040   h1->SetMarkerColor(2);
1041   h2->Draw("sames");
1042   h2->SetMarkerColor(4);
1043   h3->Draw("sames");
1044   h3->SetMarkerColor(3);
1045   return c;
1046 }
1047
1048 //------------------------------------------------------------------------------------
1049 TCanvas* DrawResMomVsMom(const char* name, const char* title, TH2* h, Int_t nBins, TF1* f2, const char* fitting)
1050 {
1051   /// generic function to draw and eventually fit momentum residuals versus momentum
1052   TLegend* l = new TLegend(0.15,0.25,0.3,0.85);
1053   TCanvas* c = new TCanvas(name, title);
1054   c->cd();
1055   TH1D* proj = 0x0;
1056   h->Sumw2();
1057   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
1058   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1059     if (f2) cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
1060     proj = h->ProjectionY(Form("%s_%d",h->GetName(),i/rebinFactorX),i-rebinFactorX+1,i);
1061     if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
1062     proj->Draw((i==rebinFactorX)?"hist":"histsames");
1063     proj->SetLineColor(i/rebinFactorX);
1064     if (f2) {
1065       f2->SetParameters(0.2,0.,1.,1.);
1066       f2->SetLineColor(i/rebinFactorX);
1067       proj->Fit("f2","WWNQ","sames");
1068       Double_t fwhm = f2->GetParameter(0);
1069       Double_t sigma = f2->GetParameter(3);
1070       Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
1071       Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/proj->GetBinWidth(1)),1);
1072       while (proj->GetNbinsX()%rebin!=0) rebin--;
1073       proj->Rebin(rebin);
1074       proj->Scale(1./rebin);
1075       proj->Fit("f2","Q","sames");
1076     } else proj->SetLineWidth(2);
1077     Double_t p = 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1078     l->AddEntry(proj,Form("%5.1f GeV",p));
1079   }
1080   if (f2) cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
1081   l->Draw("same");
1082   return c;
1083 }
1084