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