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 "TMath.h"
  30 #include "TClonesArray.h"
  31 #include "TH1.h"
  32 #include "TH2.h"
  33 #include "TH3.h"
  34 #include "TGraphErrors.h"
  35 #include "TF1.h"
  36 #include "TFile.h"
  37 #include "TCanvas.h"
  38 #include "TLegend.h"
  39
  40 // STEER includes
  41 #include "AliCDBManager.h"
  42 #include "AliLog.h"
  43
  44 // MUON includes
  45 #include "AliMUONCDB.h"
  46 #include "AliMUONConstants.h"
  47 #include "AliMUONTrack.h"
  48 #include "AliMUONRecoCheck.h"
  49 #include "AliMUONTrackParam.h"
  50 #include "AliMUONRecoParam.h"
  51 #include "AliMUONVTrackStore.h"
  52 #include "AliMUONVCluster.h"
  53 #include "AliMUONTrackExtrap.h"
  54
  55 Double_t langaufun(Double_t *x, Double_t *par);
  56
  57 //------------------------------------------------------------------------------------
  58 void MUONRecoCheck (Int_t nEvent = -1, const char* pathSim="./generated/", const char* esdFileName="AliESDs.root",
  59                     const char* ocdbPath = "local://$ALICE_ROOT/OCDB")
  60 {
  61
  62   AliLog::SetClassDebugLevel("AliMCEvent",-1);
  63
  64   // ###################################### define histograms ###################################### //
  65   // File for histograms and histogram booking
  66   TFile *histoFile = new TFile("MUONRecoCheck.root", "RECREATE");
  67
  68   TH1F *hReconstructible = new TH1F("hReconstructible"," Nb of reconstructible tracks ",15,-0.5,14.5);
  69   TH1F *hReco = new TH1F("hReco"," Nb of reconstructed tracks / evt",15,-0.5,14.5);
  70   TH1F *hNClusterComp = new TH1F("hNClusterComp"," Nb of compatible clusters / track ",15,-0.5,14.5);
  71   TH1F *hTrackRefID = new TH1F("hTrackRefID"," track reference ID ",100,-0.5,99.5);
  72
  73   // momentum resolution at vertex
  74   histoFile->mkdir("momentumAtVertex","momentumAtVertex");
  75   histoFile->cd("momentumAtVertex");
  76
  77   const Int_t pNBins = 30;
  78   const Double_t pEdges[2] = {0., 300.};
  79   const Int_t deltaPAtVtxNBins = 250;
  80   const Double_t deltaPAtVtxEdges[2] = {-35., 15.};
  81
  82   TH1F *hResMomVertex = new TH1F("hResMomVertex"," delta P at vertex;#Delta_{p} (GeV/c)",deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
  83
  84   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]);
  85   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]);
  86   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]);
  87   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]);
  88
  89   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)",100,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
  90   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)",100,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
  91   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)",100,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
  92
  93   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]);
  94   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]);
  95   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]);
  96
  97   TGraphErrors* gMeanResMomVertexVsMom = new TGraphErrors(pNBins);
  98   gMeanResMomVertexVsMom->SetName("gMeanResMomVertexVsMom");
  99   gMeanResMomVertexVsMom->SetTitle("<#Delta_{p}> at vertex versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
 100   TGraphErrors* gMostProbResMomVertexVsMom = new TGraphErrors(pNBins);
 101   gMostProbResMomVertexVsMom->SetName("gMostProbResMomVertexVsMom");
 102   gMostProbResMomVertexVsMom->SetTitle("Most probable #Delta_{p} at vertex versus p;p (GeV/c);Most prob. #Delta_{p} (GeV/c)");
 103   TGraphErrors* gSigmaResMomVertexVsMom = new TGraphErrors(pNBins);
 104   gSigmaResMomVertexVsMom->SetName("gSigmaResMomVertexVsMom");
 105   gSigmaResMomVertexVsMom->SetTitle("#Delta_{p}/p at vertex versus p;p (GeV/c);#sigma_{p}/p (%)");
 106   TF1 *f2 = new TF1("f2",langaufun,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1],4);
 107
 108   // momentum resolution at first cluster
 109   histoFile->mkdir("momentumAtFirstCluster","momentumAtFirstCluster");
 110   histoFile->cd("momentumAtFirstCluster");
 111
 112   const Int_t deltaPAtFirstClNBins = 500;
 113   const Double_t deltaPAtFirstClEdges[2] = {-25., 25.};
 114
 115   TH1F *hResMomFirstCluster = new TH1F("hResMomFirstCluster"," delta P at first cluster;#Delta_{p} (GeV/c)",deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
 116   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]);
 117
 118   TGraphErrors* gMeanResMomFirstClusterVsMom = new TGraphErrors(pNBins);
 119   gMeanResMomFirstClusterVsMom->SetName("gMeanResMomFirstClusterVsMom");
 120   gMeanResMomFirstClusterVsMom->SetTitle("<#Delta_{p}> at first cluster versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
 121   TGraphErrors* gSigmaResMomFirstClusterVsMom = new TGraphErrors(pNBins);
 122   gSigmaResMomFirstClusterVsMom->SetName("gSigmaResMomFirstClusterVsMom");
 123   gSigmaResMomFirstClusterVsMom->SetTitle("#Delta_{p}/p at first cluster versus p;p (GeV/c);#sigma_{p}/p (%)");
 124   TF1* f = new TF1("f","gausn");
 125
 126   // cluster resolution
 127   histoFile->mkdir("clusters","clusters");
 128   histoFile->cd("clusters");
 129
 130   TH1F* hResidualXInCh[AliMUONConstants::NTrackingCh()];
 131   TH1F* hResidualYInCh[AliMUONConstants::NTrackingCh()];
 132   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
 133     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.);
 134     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);
 135   }
 136
 137   TGraphErrors* gResidualXPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
 138   gResidualXPerChMean->SetName("gResidualXPerChMean");
 139   gResidualXPerChMean->SetTitle("cluster-track residual-X per Ch: mean;chamber ID;<#Delta_{X}> (cm)");
 140   gResidualXPerChMean->SetMarkerStyle(kFullDotLarge);
 141   TGraphErrors* gResidualYPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
 142   gResidualYPerChMean->SetName("gResidualYPerChMean");
 143   gResidualYPerChMean->SetTitle("cluster-track residual-Y per Ch: mean;chamber ID;<#Delta_{Y}> (cm)");
 144   gResidualYPerChMean->SetMarkerStyle(kFullDotLarge);
 145   TGraphErrors* gResidualXPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
 146   gResidualXPerChSigma->SetName("gResidualXPerChSigma");
 147   gResidualXPerChSigma->SetTitle("cluster-track residual-X per Ch: sigma;chamber ID;#sigma_{X} (cm)");
 148   gResidualXPerChSigma->SetMarkerStyle(kFullDotLarge);
 149   TGraphErrors* gResidualYPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
 150   gResidualYPerChSigma->SetName("gResidualYPerChSigma");
 151   gResidualYPerChSigma->SetTitle("cluster-track residual-Y per Ch: sigma;chamber ID;#sigma_{Y} (cm)");
 152   gResidualYPerChSigma->SetMarkerStyle(kFullDotLarge);
 153
 154   // ###################################### initialize ###################################### //
 155   AliMUONRecoCheck rc(esdFileName, pathSim);
 156
 157   // load necessary data from OCDB
 158   AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
 159   AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
 160   if (!AliMUONCDB::LoadField()) return;
 161   AliMUONTrackExtrap::SetField();
 162   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
 163   if (!recoParam) return;
 164
 165   // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
 166   Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
 167   // compute the mask of requested stations from recoParam
 168   UInt_t requestedStationMask = 0;
 169   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
 170   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
 171   Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
 172
 173   Int_t nevents = rc.NumberOfEvents();
 174
 175   if (nevents < nEvent || nEvent < 0) nEvent = nevents;
 176
 177   Int_t ievent;
 178   Int_t nReconstructibleTracks = 0;
 179   Int_t nReconstructedTracks = 0;
 180   Int_t nReconstructibleTracksCheck = 0;
 181   AliMUONTrackParam *trackParam;
 182   Double_t x1,y1,z1,pX1,pY1,pZ1,p1,pT1;
 183   Double_t x2,y2,z2,pX2,pY2,pZ2,p2,pT2;
 184   Double_t xAbs,yAbs,dAbs,aAbs,aMC;
 185
 186   // ###################################### fill histograms ###################################### //
 187   for (ievent=0; ievent<nEvent; ievent++)
 188   {
 189     if ((ievent+1)%100 == 0) cout<<"\rEvent processing... "<<ievent+1<<flush;
 190
 191     AliMUONVTrackStore* trackStore = rc.ReconstructedTracks(ievent, kFALSE);
 192     AliMUONVTrackStore* trackRefStore = rc.ReconstructibleTracks(ievent, requestedStationMask, request2ChInSameSt45);
 193
 194     hReconstructible->Fill(trackRefStore->GetSize());
 195     hReco->Fill(trackStore->GetSize());
 196
 197     nReconstructibleTracks += trackRefStore->GetSize();
 198     nReconstructedTracks += trackStore->GetSize();
 199
 200     // loop over trackRef
 201     TIter next(trackRefStore->CreateIterator());
 202     AliMUONTrack* trackRef;
 203     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) )
 204     {
 205
 206       hTrackRefID->Fill(trackRef->GetMCLabel());
 207
 208       AliMUONTrack* trackMatched = 0x0;
 209       Int_t nMatchClusters = 0;
 210
 211       // loop over trackReco and look for compatible track
 212       TIter next2(trackStore->CreateIterator());
 213       AliMUONTrack* trackReco;
 214       while ( ( trackReco = static_cast<AliMUONTrack*>(next2()) ) )
 215       {
 216
 217         // check if trackReco is compatible with trackRef
 218         Int_t n = 0;
 219         if (trackReco->Match(*trackRef, sigmaCut, nMatchClusters)) {
 220           trackMatched = trackReco;
 221           nMatchClusters = n;
 222           break;
 223         }
 224
 225       }
 226
 227       if (trackMatched) { // tracking requirements verified, track is found
 228         nReconstructibleTracksCheck++;
 229         hNClusterComp->Fill(nMatchClusters);
 230
 231         // compute track position at the end of the absorber
 232         AliMUONTrackParam trackParamAtAbsEnd(*((AliMUONTrackParam*)trackMatched->GetTrackParamAtCluster()->First()));
 233         AliMUONTrackExtrap::ExtrapToZ(&trackParamAtAbsEnd, AliMUONConstants::AbsZEnd());
 234         xAbs = trackParamAtAbsEnd.GetNonBendingCoor();
 235         yAbs = trackParamAtAbsEnd.GetBendingCoor();
 236         dAbs = TMath::Sqrt(xAbs*xAbs + yAbs*yAbs);
 237         aAbs = TMath::ATan(-dAbs/AliMUONConstants::AbsZEnd()) * TMath::RadToDeg();
 238
 239         trackParam = trackRef->GetTrackParamAtVertex();
 240         x1 = trackParam->GetNonBendingCoor();
 241         y1 = trackParam->GetBendingCoor();
 242         z1 = trackParam->GetZ();
 243         pX1 = trackParam->Px();
 244         pY1 = trackParam->Py();
 245         pZ1 = trackParam->Pz();
 246         p1  = trackParam->P();
 247         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
 248         aMC = TMath::ATan(-pT1/pZ1) * TMath::RadToDeg();
 249
 250         trackParam = trackMatched->GetTrackParamAtVertex();
 251         x2 = trackParam->GetNonBendingCoor();
 252         y2 = trackParam->GetBendingCoor();
 253         z2 = trackParam->GetZ();
 254         pX2 = trackParam->Px();
 255         pY2 = trackParam->Py();
 256         pZ2 = trackParam->Pz();
 257         p2  = trackParam->P();
 258         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
 259
 260         hResMomVertex->Fill(p2-p1);
 261         hResMomVertexVsMom->Fill(p1,p2-p1);
 262         hResMomVertexVsAngleVsMom->Fill(p1,aAbs,p2-p1);
 263         if (aAbs > 2. && aAbs < 3.) hResMomVertexVsMom_2_3_Deg->Fill(p1,p2-p1);
 264         else if (aAbs >= 3. && aAbs < 10.) hResMomVertexVsMom_3_10_Deg->Fill(p1,p2-p1);
 265         if (aMC < 2.) {
 266           hResMomVertexVsMom_0_2_DegMC->Fill(p1,p2-p1);
 267           hResMomVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,p2-p1);
 268         }
 269         else if (aMC >= 2. && aMC < 3) hResMomVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,p2-p1);
 270         else if (aMC >= 3. && aMC < 10.) hResMomVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,p2-p1);
 271         hResMomVertexVsAngle->Fill(aAbs,p2-p1);
 272         hResMomVertexVsMCAngle->Fill(aMC,p2-p1);
 273
 274         trackParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
 275         x1 = trackParam->GetNonBendingCoor();
 276         y1 = trackParam->GetBendingCoor();
 277         z1 = trackParam->GetZ();
 278         pX1 = trackParam->Px();
 279         pY1 = trackParam->Py();
 280         pZ1 = trackParam->Pz();
 281         p1  = trackParam->P();
 282         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
 283
 284         trackParam = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
 285         x2 = trackParam->GetNonBendingCoor();
 286         y2 = trackParam->GetBendingCoor();
 287         z2 = trackParam->GetZ();
 288         pX2 = trackParam->Px();
 289         pY2 = trackParam->Py();
 290         pZ2 = trackParam->Pz();
 291         p2  = trackParam->P();
 292         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
 293
 294         hResMomFirstCluster->Fill(p2-p1);
 295         hResMomFirstClusterVsMom->Fill(p1,p2-p1);
 296
 297         // Fill residuals
 298         // Loop over clusters of first track
 299         AliMUONTrackParam* trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
 300         while (trackParamAtCluster1) {
 301           AliMUONVCluster* cluster1 = trackParamAtCluster1->GetClusterPtr();
 302           AliMUONTrackParam* trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
 303           while (trackParamAtCluster2) {
 304             AliMUONVCluster* cluster2 = trackParamAtCluster2->GetClusterPtr();
 305             if (cluster1->GetDetElemId() == cluster2->GetDetElemId()) {
 306               hResidualXInCh[cluster1->GetChamberId()]->Fill(cluster1->GetX() - cluster2->GetX());
 307               hResidualYInCh[cluster1->GetChamberId()]->Fill(cluster1->GetY() - cluster2->GetY());
 308               break;
 309             }
 310             trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->After(trackParamAtCluster2);
 311           }
 312           trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->After(trackParamAtCluster1);
 313         }
 314
 315       }
 316
 317     } // end loop track ref.
 318
 319   } // end loop on event
 320   cout<<"\rEvent processing... "<<nevents<<" done"<<endl;
 321
 322   // ###################################### compute stuff ###################################### //
 323   // compute momentum resolution at vertex versus p
 324   Int_t rebinFactorX = TMath::Max(hResMomVertexVsMom->GetNbinsX()/pNBins, 1);
 325   for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom->GetNbinsX(); i+=rebinFactorX) {
 326     cout<<"\rFitting momentum residuals at vertex... "<<i/rebinFactorX<<"/"<<pNBins<<flush;
 327     TH1D *tmp = hResMomVertexVsMom->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
 328     Double_t p = 0.5 * (hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom->GetBinLowEdge(i+1));
 329     f2->SetParameters(0.2,0.,(Double_t)tmp->GetEntries(),1.);
 330     tmp->Fit("f2","WWNQ");
 331     Double_t fwhm = f2->GetParameter(0);
 332     Double_t sigma = f2->GetParameter(3);
 333     Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
 334     Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/tmp->GetBinWidth(1)),1);
 335     while (deltaPAtVtxNBins%rebin!=0) rebin--;
 336     tmp->Rebin(rebin);
 337     tmp->Fit("f2","NQ");
 338     fwhm = f2->GetParameter(0);
 339     sigma = f2->GetParameter(3);
 340     sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
 341     Double_t fwhmErr = f2->GetParError(0);
 342     Double_t sigmaErr = f2->GetParError(3);
 343     Double_t sigmaPErr = TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + fwhm*fwhm*fwhmErr*fwhmErr/(64.*log(2.)*log(2.))) / sigmaP;
 344     gMeanResMomVertexVsMom->SetPoint(i/rebinFactorX-1,p,tmp->GetMean());
 345     gMeanResMomVertexVsMom->SetPointError(i/rebinFactorX-1,hResMomVertexVsMom->GetBinWidth(i),tmp->GetMeanError());
 346     gMostProbResMomVertexVsMom->SetPoint(i/rebinFactorX-1,p,-f2->GetParameter(1));
 347     gMostProbResMomVertexVsMom->SetPointError(i/rebinFactorX-1,hResMomVertexVsMom->GetBinWidth(i),f2->GetParError(1));
 348     gSigmaResMomVertexVsMom->SetPoint(i/rebinFactorX-1,p,100.*sigmaP/p);
 349     gSigmaResMomVertexVsMom->SetPointError(i/rebinFactorX-1,hResMomVertexVsMom->GetBinWidth(i),100.*sigmaPErr/p);
 350     delete tmp;
 351   }
 352   cout<<"\rFitting momentum residuals at vertex... "<<pNBins<<"/"<<pNBins<<endl;
 353
 354   // compute momentum resolution at first cluster versus p
 355   rebinFactorX = TMath::Max(hResMomFirstClusterVsMom->GetNbinsX()/pNBins, 1);
 356   for (Int_t i = rebinFactorX; i <= hResMomFirstClusterVsMom->GetNbinsX(); i+=rebinFactorX) {
 357     cout<<"\rFitting momentum residuals at first cluster... "<<i/rebinFactorX<<"/"<<pNBins<<flush;
 358     TH1D *tmp = hResMomFirstClusterVsMom->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
 359     Double_t p = 0.5 * (hResMomFirstClusterVsMom->GetBinLowEdge(i-rebinFactorX+1) + hResMomFirstClusterVsMom->GetBinLowEdge(i+1));
 360     f->SetParameters(tmp->GetEntries(),0.,1.);
 361     tmp->Fit("f","WWNQ");
 362     Int_t rebin = TMath::Max(Int_t(0.5*f->GetParameter(2)/tmp->GetBinWidth(1)),1);
 363     while (deltaPAtFirstClNBins%rebin!=0) rebin--;
 364     tmp->Rebin(rebin);
 365     tmp->Fit("f","NQ");
 366     gMeanResMomFirstClusterVsMom->SetPoint(i/rebinFactorX-1,p,f->GetParameter(1));
 367     gMeanResMomFirstClusterVsMom->SetPointError(i/rebinFactorX-1,hResMomFirstClusterVsMom->GetBinWidth(i),f->GetParError(1));
 368     gSigmaResMomFirstClusterVsMom->SetPoint(i/rebinFactorX-1,p,100.*f->GetParameter(2)/p);
 369     gSigmaResMomFirstClusterVsMom->SetPointError(i/rebinFactorX-1,hResMomFirstClusterVsMom->GetBinWidth(i),100.*f->GetParError(2)/p);
 370     delete tmp;
 371   }
 372   cout<<"\rFitting momentum residuals at first cluster... "<<pNBins<<"/"<<pNBins<<endl;
 373
 374   // compute residual mean and dispersion
 375   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
 376     hResidualXInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualXInCh[i]->GetRMS(), 3.*hResidualXInCh[i]->GetRMS());
 377     gResidualXPerChMean->SetPoint(i, i+1, hResidualXInCh[i]->GetMean());
 378     gResidualXPerChMean->SetPointError(i, 0., hResidualXInCh[i]->GetMeanError());
 379     gResidualXPerChSigma->SetPoint(i, i+1, hResidualXInCh[i]->GetRMS());
 380     gResidualXPerChSigma->SetPointError(i, 0., hResidualXInCh[i]->GetRMSError());
 381     hResidualXInCh[i]->GetXaxis()->SetRange(0,0);
 382     hResidualYInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualYInCh[i]->GetRMS(), 3.*hResidualYInCh[i]->GetRMS());
 383     gResidualYPerChMean->SetPoint(i, i+1, hResidualYInCh[i]->GetMean());
 384     gResidualYPerChMean->SetPointError(i, 0., hResidualYInCh[i]->GetMeanError());
 385     gResidualYPerChSigma->SetPoint(i, i+1, hResidualYInCh[i]->GetRMS());
 386     gResidualYPerChSigma->SetPointError(i, 0., hResidualYInCh[i]->GetRMSError());
 387     hResidualYInCh[i]->GetXaxis()->SetRange(0,0);
 388   }
 389
 390   // ###################################### display histograms ###################################### //
 391   // diplay momentum residual for different angular region
 392   TCanvas cResMom("cResMom", "momentum residual at vertex in 3 angular regions");
 393   cResMom.cd();
 394   hResMomVertex->Draw();
 395   TH1D *hResMomVertex_0_2_Deg = hResMomVertexVsAngle->ProjectionY("hResMomVertex_0_2_Deg",1,2);
 396   hResMomVertex_0_2_Deg->Draw("sames");
 397   hResMomVertex_0_2_Deg->SetLineColor(2);
 398   TH1D *hResMomVertex_2_3_Deg = hResMomVertexVsAngle->ProjectionY("hResMomVertex_2_3_Deg",3,3);
 399   hResMomVertex_2_3_Deg->Draw("sames");
 400   hResMomVertex_2_3_Deg->SetLineColor(4);
 401   TH1D *hResMomVertex_3_10_Deg = hResMomVertexVsAngle->ProjectionY("hResMomVertex_3_10_Deg",4,10);
 402   hResMomVertex_3_10_Deg->Draw("sames");
 403   hResMomVertex_3_10_Deg->SetLineColor(3);
 404
 405   // diplay momentum residual for different angular region
 406   TCanvas cResMomMC("cResMomMC", "momentum residual at vertex in 3 MC angular regions");
 407   cResMomMC.cd();
 408   hResMomVertex->Draw();
 409   TH1D *hResMomVertex_0_2_DegMC = hResMomVertexVsMCAngle->ProjectionY("hResMomVertex_0_2_DegMC",1,2);
 410   hResMomVertex_0_2_DegMC->Draw("sames");
 411   hResMomVertex_0_2_DegMC->SetLineColor(2);
 412   TH1D *hResMomVertex_2_3_DegMC = hResMomVertexVsMCAngle->ProjectionY("hResMomVertex_2_3_DegMC",3,3);
 413   hResMomVertex_2_3_DegMC->Draw("sames");
 414   hResMomVertex_2_3_DegMC->SetLineColor(4);
 415   TH1D *hResMomVertex_3_10_DegMC = hResMomVertexVsMCAngle->ProjectionY("hResMomVertex_3_10_DegMC",4,10);
 416   hResMomVertex_3_10_DegMC->Draw("sames");
 417   hResMomVertex_3_10_DegMC->SetLineColor(3);
 418
 419   // diplay momentum residual versus position at absorber end for different MC angular region
 420   TCanvas cResMomVsPos("cResMomVsPos", "momentum residual at vertex versus position at absorber end in 3 MC angular regions");
 421   cResMomVsPos.cd();
 422   hResMomVertexVsPosAbsEnd_0_2_DegMC->Draw();
 423   hResMomVertexVsPosAbsEnd_0_2_DegMC->SetMarkerColor(2);
 424   hResMomVertexVsPosAbsEnd_2_3_DegMC->Draw("sames");
 425   hResMomVertexVsPosAbsEnd_2_3_DegMC->SetMarkerColor(4);
 426   hResMomVertexVsPosAbsEnd_3_10_DegMC->Draw("sames");
 427   hResMomVertexVsPosAbsEnd_3_10_DegMC->SetMarkerColor(3);
 428
 429   // diplay momentum residual of tracks between 2 and 3 deg. for different momentum values
 430   Int_t pNBinsShown = 10;
 431   TLegend lResMom_2_3_Deg(0.15,0.25,0.3,0.85);
 432   TCanvas cResMom_2_3_Deg("cResMom_2_3_Deg", "momentum residual for tracks between 2 and 3 degrees");
 433   cResMom_2_3_Deg.cd();
 434   TH1D* proj = 0x0;
 435   hResMomVertexVsMom_2_3_Deg->Sumw2();
 436   rebinFactorX = TMath::Max(hResMomVertexVsMom_2_3_Deg->GetNbinsX()/pNBinsShown, 1);
 437   for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom_2_3_Deg->GetNbinsX(); i+=rebinFactorX) {
 438     cout<<"\rFitting momentum residuals at vertex (tracks in [2,3] deg.)... "<<i/rebinFactorX<<"/"<<pNBinsShown<<flush;
 439     proj = hResMomVertexVsMom_2_3_Deg->ProjectionY(Form("hRes23_%d",i/rebinFactorX),i-rebinFactorX+1,i);
 440     if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
 441     proj->Draw((i==rebinFactorX)?"hist":"histsames");
 442     proj->SetLineColor(i/rebinFactorX);
 443     f2->SetParameters(0.2,0.,1.,1.);
 444     f2->SetLineColor(i/rebinFactorX);
 445     proj->Fit("f2","WWNQ","sames");
 446     Double_t fwhm = f2->GetParameter(0);
 447     Double_t sigma = f2->GetParameter(3);
 448     Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
 449     Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/proj->GetBinWidth(1)),1);
 450     while (deltaPAtVtxNBins%rebin!=0) rebin--;
 451     proj->Rebin(rebin);
 452     proj->Scale(1./rebin);
 453     proj->Fit("f2","Q","sames");
 454     Double_t p = 0.5 * (hResMomVertexVsMom_2_3_Deg->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom_2_3_Deg->GetBinLowEdge(i+1));
 455     lResMom_2_3_Deg.AddEntry(proj,Form("%5.1f GeV",p));
 456   }
 457   cout<<"\rFitting momentum residuals at vertex (tracks in [2,3] deg.)... "<<pNBinsShown<<"/"<<pNBinsShown<<endl;
 458   lResMom_2_3_Deg.Draw("same");
 459
 460   // diplay momentum residual of tracks between 3 and 10 deg. for different momentum values
 461   pNBinsShown = 10;
 462   TLegend lResMom_3_10_Deg(0.15,0.25,0.3,0.85);
 463   TCanvas cResMom_3_10_Deg("cResMom_3_10_Deg", "momentum residual for tracks between 3 and 10 degrees");
 464   cResMom_3_10_Deg.cd();
 465   proj = 0x0;
 466   hResMomVertexVsMom_3_10_Deg->Sumw2();
 467   rebinFactorX = TMath::Max(hResMomVertexVsMom_3_10_Deg->GetNbinsX()/pNBinsShown, 1);
 468   for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom_3_10_Deg->GetNbinsX(); i+=rebinFactorX) {
 469     cout<<"\rFitting momentum residuals at vertex (tracks in [3,10] deg.)... "<<i/rebinFactorX<<"/"<<pNBinsShown<<flush;
 470     proj = hResMomVertexVsMom_3_10_Deg->ProjectionY(Form("hRes310_%d",i/rebinFactorX),i-rebinFactorX+1,i);
 471     if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
 472     proj->Draw((i==rebinFactorX)?"hist":"histsames");
 473     proj->SetLineColor(i/rebinFactorX);
 474     f2->SetParameters(0.2,0.,1.,1.);
 475     f2->SetLineColor(i/rebinFactorX);
 476     proj->Fit("f2","WWNQ","sames");
 477     Double_t fwhm = f2->GetParameter(0);
 478     Double_t sigma = f2->GetParameter(3);
 479     Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
 480     Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/proj->GetBinWidth(1)),1);
 481     while (deltaPAtVtxNBins%rebin!=0) rebin--;
 482     proj->Rebin(rebin);
 483     proj->Scale(1./rebin);
 484     proj->Fit("f2","Q","sames");
 485     Double_t p = 0.5 * (hResMomVertexVsMom_3_10_Deg->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom_3_10_Deg->GetBinLowEdge(i+1));
 486     lResMom_3_10_Deg.AddEntry(proj,Form("%5.1f GeV",p));
 487   }
 488   cout<<"\rFitting momentum residuals at vertex (tracks in [3,10] deg.)... "<<pNBinsShown<<"/"<<pNBinsShown<<endl;
 489   lResMom_3_10_Deg.Draw("same");
 490
 491   // diplay momentum residuals of tracks with MC angle < 2 deg. for different momentum values
 492   pNBinsShown = 5;
 493   TLegend lResMom_0_2_DegMC(0.15,0.25,0.3,0.85);
 494   TCanvas cResMom_0_2_DegMC("cResMom_0_2_DegMC", "momentum residuals for tracks with MC angle < 2 degrees");
 495   cResMom_0_2_DegMC.cd();
 496   proj = 0x0;
 497   hResMomVertexVsMom_0_2_DegMC->Sumw2();
 498   rebinFactorX = TMath::Max(hResMomVertexVsMom_0_2_DegMC->GetNbinsX()/pNBinsShown, 1);
 499   for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom_0_2_DegMC->GetNbinsX(); i+=rebinFactorX) {
 500     proj = hResMomVertexVsMom_0_2_DegMC->ProjectionY(Form("hRes02_%d",i/rebinFactorX),i-rebinFactorX+1,i);
 501     if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
 502     proj->Draw((i==rebinFactorX)?"hist":"histsames");
 503     proj->SetLineColor(i/rebinFactorX);
 504     proj->SetLineWidth(2);
 505     Double_t p = 0.5 * (hResMomVertexVsMom_0_2_DegMC->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom_0_2_DegMC->GetBinLowEdge(i+1));
 506     lResMom_0_2_DegMC.AddEntry(proj,Form("%5.1f GeV",p));
 507   }
 508   lResMom_0_2_DegMC.Draw("same");
 509
 510   // ###################################### save histogram ###################################### //
 511   histoFile->Write();
 512
 513   histoFile->cd("momentumAtVertex");
 514   gMeanResMomVertexVsMom->Write();
 515   gMostProbResMomVertexVsMom->Write();
 516   gSigmaResMomVertexVsMom->Write();
 517   cResMom.Write();
 518   cResMomMC.Write();
 519   cResMomVsPos.Write();
 520   cResMom_2_3_Deg.Write();
 521   cResMom_3_10_Deg.Write();
 522   cResMom_0_2_DegMC.Write();
 523
 524   histoFile->cd("momentumAtFirstCluster");
 525   gMeanResMomFirstClusterVsMom->Write();
 526   gSigmaResMomFirstClusterVsMom->Write();
 527
 528   histoFile->cd("clusters");
 529   gResidualXPerChMean->Write();
 530   gResidualXPerChSigma->Write();
 531   gResidualYPerChMean->Write();
 532   gResidualYPerChSigma->Write();
 533
 534   histoFile->Close();
 535
 536   printf(" nb of reconstructible tracks: %d \n", nReconstructibleTracks);
 537   printf(" nb of reconstructed tracks: %d \n", nReconstructedTracks);
 538   printf(" nb of reconstructible tracks which are reconstructed: %d \n", nReconstructibleTracksCheck);
 539
 540 }
 541
 542 //------------------------------------------------------------------------------------
 543 Double_t langaufun(Double_t *x, Double_t *par) {
 544
 545   //Fit parameters:
 546   //par[0]=Width (scale) parameter of Landau density
 547   //par[1]=Most Probable (MP, location) parameter of Landau density
 548   //par[2]=Total area (integral -inf to inf, normalization constant)
 549   //par[3]=Width (sigma) of convoluted Gaussian function
 550   //
 551   //In the Landau distribution (represented by the CERNLIB approximation),
 552   //the maximum is located at x=-0.22278298 with the location parameter=0.
 553   //This shift is corrected within this function, so that the actual
 554   //maximum is identical to the MP parameter.
 555
 556   // Numeric constants
 557   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 558   Double_t mpshift  = -0.22278298;       // Landau maximum location
 559
 560   // Control constants
 561   Double_t np = 100.0; // number of convolution steps
 562   Double_t sc = 5.0;   // convolution extends to +-sc Gaussian sigmas
 563
 564   // Variables
 565   Double_t xx;
 566   Double_t mpc;
 567   Double_t fland;
 568   Double_t sum = 0.0;
 569   Double_t xlow,xupp;
 570   Double_t step;
 571   Double_t i;
 572
 573
 574   // MP shift correction
 575   mpc = par[1] - mpshift * par[0];
 576
 577   // Range of convolution integral
 578   xlow = x[0] - sc * par[3];
 579   xupp = x[0] + sc * par[3];
 580
 581   step = (xupp-xlow) / np;
 582
 583   // Convolution integral of Landau and Gaussian by sum
 584   for(i=1.0; i<=np/2; i++) {
 585     xx = xlow + (i-.5) * step;
 586     //change x -> -x because the tail of the Landau is at the left here...
 587     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
 588     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 589
 590     //change x -> -x because the tail of the Landau is at the left here...
 591     xx = xupp - (i-.5) * step;
 592     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
 593     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 594   }
 595
 596   return (par[2] * step * sum * invsq2pi / par[3]);
 597 }
 598