- fix data container naming issues
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackAnalysis.cxx
CommitLineData
e806e863 1//====================================================================================================================================================
2//
3// Class for the analysis of the ALICE muon forward tracks (MUON + MFT)
4//
5// Contact author: antonio.uras@cern.ch
6//
7//====================================================================================================================================================
8
9#include "TObject.h"
10#include "TClonesArray.h"
11#include "AliMuonForwardTrack.h"
12#include "AliMuonForwardTrackPair.h"
13#include "TMatrixD.h"
14#include "TTree.h"
15#include "TH1D.h"
16#include "AliLog.h"
17#include "TFile.h"
18#include "TParticle.h"
19#include "AliMUONTrackParam.h"
20#include "AliMUONTrackExtrap.h"
21#include "TGeoManager.h"
22#include "TRandom.h"
23#include "TLorentzVector.h"
24#include "TDatabasePDG.h"
b5ab1ac4 25#include "TGraph.h"
e806e863 26#include "AliMuonForwardTrackAnalysis.h"
27
28ClassImp(AliMuonForwardTrackAnalysis)
29
30//====================================================================================================================================================
31
32AliMuonForwardTrackAnalysis::AliMuonForwardTrackAnalysis():
33 TObject(),
34 fInputDir(0),
35 fOutputDir(0),
7e3dd1af 36 fInputTreeWithBranson(0x0),
37 fInputTreeWithoutBranson(0x0),
38 fMuonForwardTracksWithBranson(0),
39 fMuonForwardTrackPairsWithBranson(0),
40 fMuonForwardTracksWithoutBranson(0),
41 fMuonForwardTrackPairsWithoutBranson(0),
42 fMFTTrackWithBranson(0),
43 fMFTTrackWithoutBranson(0),
e806e863 44 fMFTTrack(0),
7e3dd1af 45 fMFTTrackPairWithBranson(0),
46 fMFTTrackPairWithoutBranson(0),
e806e863 47 fMFTTrackPair(0),
48 fMCRefTrack(0),
49 fEv(0),
50 fFirstEvent(-1),
51 fLastEvent(-1),
52 fNTracksOfEvent(0),
53 fNTracksAnalyzedOfEvent(0),
54 fNTracksAnalyzed(0),
55 fNPairsOfEvent(0),
56 fNPairsAnalyzedOfEvent(0),
d8c2cc3e 57 fNTracksAnalyzedOfEventAfterCut(0),
58 fNPairsAnalyzedOfEventAfterCut(0),
e806e863 59 fHistOffsetSingleMuonsX(0x0),
60 fHistOffsetSingleMuonsY(0x0),
61 fHistOffsetSingleMuons(0x0),
62 fHistWOffsetSingleMuons(0x0),
63 fHistErrorSingleMuonsX(0x0),
64 fHistErrorSingleMuonsY(0x0),
d8c2cc3e 65 fHistZOriginSingleMuonsMC(0x0),
cabfe25b 66 fHistZROriginSingleMuonsMC(0x0),
67 fHistSingleMuonsPtRapidity(0x0),
b5ab1ac4 68 fHistSingleMuonsOffsetChi2(0x0),
d8c2cc3e 69 fHistRapidityPtMuonPairs(0x0),
cabfe25b 70 fHistMassMuonPairsMCVsPt(0x0),
71 fHistMassMuonPairsVsPt(0x0),
72 fHistMassMuonPairsWithoutMFTVsPt(0x0),
73 fHistMassMuonPairsVsPtLSp(0x0),
74 fHistMassMuonPairsWithoutMFTVsPtLSp(0x0),
75 fHistMassMuonPairsVsPtLSm(0x0),
76 fHistMassMuonPairsWithoutMFTVsPtLSm(0x0),
77 fHistWOffsetMuonPairsAtPrimaryVtxVsPt(0x0),
78 fHistWOffsetMuonPairsAtPCAVsPt(0x0),
79 fHistDistancePrimaryVtxPCAVsPt(0x0),
80 fHistPCAQualityVsPt(0x0),
81 fHistPseudoProperDecayLengthVsPt(0x0),
82 fHistDimuonVtxResolutionXVsPt(0x0),
83 fHistDimuonVtxResolutionYVsPt(0x0),
84 fHistDimuonVtxResolutionZVsPt(0x0),
85 fTrueMass(0.),
e806e863 86 fMassMin(0),
cabfe25b 87 fMassMax(9999),
e806e863 88 fSingleMuonAnalysis(1),
89 fMuonPairAnalysis(1),
e806e863 90 fOption(0),
cabfe25b 91 fTriggerLevel(0),
b5ab1ac4 92 fXVertResMC(50.e-4),
93 fYVertResMC(50.e-4),
94 fZVertResMC(50.e-4),
7e3dd1af 95 fPrimaryVtxX(0.),
96 fPrimaryVtxY(0.),
97 fPrimaryVtxZ(0.),
e806e863 98 fMaxNWrongClustersMC(999),
7e3dd1af 99 fPtMinSingleMuons(0),
100 fUseBransonForCut(kFALSE),
101 fUseBransonForKinematics(kFALSE),
cabfe25b 102 fCorrelateCutOnOffsetChi2(kFALSE),
103 fMaxChi2SingleMuons(1.e9),
104 fMaxOffsetSingleMuons(1.e9),
105 fMaxWOffsetMuonPairsAtPrimaryVtx(1.e9),
106 fMaxWOffsetMuonPairsAtPCA(1.e9),
107 fMaxDistancePrimaryVtxPCA(1.e9),
108 fMinPCAQuality(0.)
e806e863 109{
110
111 // default constructor
112
e806e863 113}
114
115//====================================================================================================================================================
116
117Bool_t AliMuonForwardTrackAnalysis::Init(Char_t *inputFileName) {
118
119 BookHistos();
120
121 TFile *inputFile = new TFile(Form("%s/%s",fInputDir.Data(),inputFileName));
122 if (!inputFile || !inputFile->IsOpen()) {
123 AliError(Form("Error opening file %s", inputFileName));
124 return kFALSE;
125 }
7e3dd1af 126 fInputTreeWithBranson = (TTree*) inputFile->Get("AliMuonForwardTracksWithBranson");
127 if (!fInputTreeWithBranson) {
128 AliError("Error reading input tree");
129 return kFALSE;
130 }
131 fInputTreeWithoutBranson = (TTree*) inputFile->Get("AliMuonForwardTracksWithoutBranson");
132 if (!fInputTreeWithoutBranson) {
e806e863 133 AliError("Error reading input tree");
134 return kFALSE;
135 }
136
3b0abc02 137 if (fFirstEvent>=fInputTreeWithBranson->GetEntriesFast()) return kFALSE;
138 else if (fFirstEvent<0 || fLastEvent<0 || fFirstEvent>fLastEvent || fFirstEvent>=fInputTreeWithBranson->GetEntriesFast()) {
e806e863 139 fFirstEvent = 0;
789c2e43 140 fLastEvent = fInputTreeWithBranson->GetEntriesFast()-1;
e806e863 141 }
142 else {
789c2e43 143 fLastEvent = TMath::Min(fLastEvent, Int_t(fInputTreeWithBranson->GetEntriesFast()-1));
e806e863 144 }
145
146 AliInfo(Form("Analysing events %d to %d", fFirstEvent, fLastEvent));
147
789c2e43 148 fMuonForwardTracksWithBranson = new TClonesArray("AliMuonForwardTrack",30);
7e3dd1af 149 fInputTreeWithBranson->SetBranchAddress("tracks", &fMuonForwardTracksWithBranson);
150
789c2e43 151 fMuonForwardTracksWithoutBranson = new TClonesArray("AliMuonForwardTrack",30);
7e3dd1af 152 fInputTreeWithoutBranson->SetBranchAddress("tracks", &fMuonForwardTracksWithoutBranson);
e806e863 153
154 TGeoManager::Import(Form("%s/geometry.root",fInputDir.Data()));
155
156 AliMUONTrackExtrap::SetField();
157
789c2e43 158 fMuonForwardTrackPairsWithBranson = new TClonesArray("AliMuonForwardTrackPair",10);
159 fMuonForwardTrackPairsWithoutBranson = new TClonesArray("AliMuonForwardTrackPair",10);
cabfe25b 160 fMuonForwardTrackPairsWithBranson -> SetOwner(kTRUE);
161 fMuonForwardTrackPairsWithoutBranson -> SetOwner(kTRUE);
162
e806e863 163 return kTRUE;
164
165}
166
167//====================================================================================================================================================
168
169Bool_t AliMuonForwardTrackAnalysis::LoadNextEvent() {
170
171 if (fEv>fLastEvent) return kFALSE;
172 if (fEv<fFirstEvent) { fEv++; return kTRUE; }
789c2e43 173 fMuonForwardTracksWithBranson -> Clear("");
174 fMuonForwardTracksWithoutBranson -> Clear("");
7e3dd1af 175 fInputTreeWithBranson->GetEvent(fEv);
176 fInputTreeWithoutBranson->GetEvent(fEv);
789c2e43 177 AliDebug(2,Form("**** analyzing event # %4d (%3d tracks) ****", fEv, fMuonForwardTracksWithBranson->GetEntriesFast()));
7e3dd1af 178
cabfe25b 179 AliInfo(Form("**** analyzing event # %6d of %6d ****", fEv, fLastEvent));
180
7e3dd1af 181 fPrimaryVtxX = gRandom->Gaus(0., fXVertResMC);
182 fPrimaryVtxY = gRandom->Gaus(0., fYVertResMC);
183 fPrimaryVtxZ = gRandom->Gaus(0., fZVertResMC);
e806e863 184
185 if (fSingleMuonAnalysis) {
186 fNTracksAnalyzedOfEvent = 0;
d8c2cc3e 187 fNTracksAnalyzedOfEventAfterCut = 0;
789c2e43 188 fNTracksOfEvent = fMuonForwardTracksWithBranson->GetEntriesFast();
e806e863 189 while (AnalyzeSingleMuon()) continue;
190 }
191
192 if (fMuonPairAnalysis) {
7e3dd1af 193 if (fMuonForwardTrackPairsWithBranson) {
789c2e43 194 fMuonForwardTrackPairsWithBranson->Clear("C");
195 fMuonForwardTrackPairsWithoutBranson->Clear("C");
e806e863 196 }
197 BuildMuonPairs();
198 fNPairsAnalyzedOfEvent = 0;
d8c2cc3e 199 fNPairsAnalyzedOfEventAfterCut = 0;
cabfe25b 200 fNPairsOfEvent = fMuonForwardTrackPairsWithBranson->GetEntriesFast();
e806e863 201 while (AnalyzeMuonPair()) continue;
202 }
203
789c2e43 204 AliDebug(2,Form("**** analyzed event # %4d (%3d tracks and %3d pairs analyzed) ****", fEv, fNTracksAnalyzedOfEventAfterCut, fNPairsAnalyzedOfEventAfterCut));
d8c2cc3e 205
e806e863 206 fEv++;
207
208 return kTRUE;
209
210}
211
212//====================================================================================================================================================
213
214Bool_t AliMuonForwardTrackAnalysis::AnalyzeSingleMuon() {
215
216 if (fNTracksAnalyzedOfEvent>=fNTracksOfEvent) return kFALSE;
217
7e3dd1af 218 fMFTTrackWithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(fNTracksAnalyzedOfEvent);
219 fMFTTrackWithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(fNTracksAnalyzedOfEvent);
220
e806e863 221 fNTracksAnalyzedOfEvent++;
7e3dd1af 222
223 Bool_t passedCut = kFALSE;
224 if (fUseBransonForCut) passedCut = PassedCutSingleMuon(fMFTTrackWithBranson);
225 else passedCut = PassedCutSingleMuon(fMFTTrackWithoutBranson);
226 if (!passedCut) return kTRUE;
227
228 if (fUseBransonForKinematics) fMFTTrack = fMFTTrackWithBranson;
229 else fMFTTrack = fMFTTrackWithoutBranson;
cabfe25b 230 if (fMFTTrack->GetMatchTrigger() < fTriggerLevel) return kTRUE;
e806e863 231 fMCRefTrack = fMFTTrack->GetMCTrackRef();
232
cabfe25b 233 if (fMaxNWrongClustersMC<6) {
234 if (!fMCRefTrack) return kTRUE;
235 if (fMFTTrack->GetNWrongClustersMC()>fMaxNWrongClustersMC) return kTRUE;
236 }
d8c2cc3e 237
238 if (fMCRefTrack) {
239 fHistZOriginSingleMuonsMC -> Fill(-1.*fMCRefTrack->Vz());
240 fHistZROriginSingleMuonsMC -> Fill(-1.*fMCRefTrack->Vz(), TMath::Sqrt(fMCRefTrack->Vx()*fMCRefTrack->Vx()+fMCRefTrack->Vy()*fMCRefTrack->Vy()));
241 }
e806e863 242
e806e863 243 AliMUONTrackParam *param = fMFTTrack->GetTrackParamAtMFTCluster(0);
7e3dd1af 244 AliMUONTrackExtrap::ExtrapToZCov(param, fPrimaryVtxZ);
e806e863 245
246 TLorentzVector pMu;
247 Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
248 Double_t energy = TMath::Sqrt(param->P()*param->P() + mMu*mMu);
249 pMu.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
250
e806e863 251 TMatrixD cov(5,5);
252 cov = param->GetCovariances();
253
254 fHistErrorSingleMuonsX -> Fill(1.e4*TMath::Sqrt(cov(0,0)));
255 fHistErrorSingleMuonsY -> Fill(1.e4*TMath::Sqrt(cov(2,2)));
256
7e3dd1af 257 Double_t dX = fMFTTrack->GetOffsetX(fPrimaryVtxX, fPrimaryVtxZ);
258 Double_t dY = fMFTTrack->GetOffsetY(fPrimaryVtxY, fPrimaryVtxZ);
e806e863 259
7e3dd1af 260 Double_t offset = fMFTTrack->GetOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
261 Double_t weightedOffset = fMFTTrack->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
e806e863 262
263 // AliDebug(2, Form("pdg code = %d\n", fMCRefTrack->GetPdgCode()));
264
265 fHistOffsetSingleMuonsX -> Fill(1.e4*dX);
266 fHistOffsetSingleMuonsY -> Fill(1.e4*dY);
cd2f51d2 267
cabfe25b 268 fHistSingleMuonsPtRapidity -> Fill(pMu.Rapidity(), pMu.Pt());
269 fHistOffsetSingleMuons -> Fill(1.e4*offset);
270 fHistWOffsetSingleMuons -> Fill(weightedOffset);
b5ab1ac4 271 Double_t chi2OverNdf = fMFTTrack->GetGlobalChi2()/Double_t(fMFTTrack->GetNMFTClusters()+fMFTTrack->GetNMUONClusters());
272 fHistSingleMuonsOffsetChi2 -> Fill(1.e4*offset, chi2OverNdf);
e806e863 273
274 fNTracksAnalyzed++;
d8c2cc3e 275 fNTracksAnalyzedOfEventAfterCut++;
e806e863 276
277 return kTRUE;
278
279}
280
281//====================================================================================================================================================
282
283Bool_t AliMuonForwardTrackAnalysis::AnalyzeMuonPair() {
284
285 if (fNPairsAnalyzedOfEvent>=fNPairsOfEvent) return kFALSE;
286
7e3dd1af 287 fMFTTrackPairWithBranson = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithBranson->At(fNPairsAnalyzedOfEvent);
288 fMFTTrackPairWithoutBranson = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithoutBranson->At(fNPairsAnalyzedOfEvent);
e806e863 289
7e3dd1af 290 fNPairsAnalyzedOfEvent++;
291
7e3dd1af 292 if (fUseBransonForKinematics) fMFTTrackPair = fMFTTrackPairWithBranson;
293 else fMFTTrackPair = fMFTTrackPairWithoutBranson;
e806e863 294
7e3dd1af 295 if ( fMFTTrackPair->GetTrack(0)->GetNWrongClustersMC()>fMaxNWrongClustersMC ||
296 fMFTTrackPair->GetTrack(1)->GetNWrongClustersMC()>fMaxNWrongClustersMC ) return kTRUE;
e806e863 297
cabfe25b 298 if (fOption==kResonanceOnly && !fMFTTrackPair->IsResonance()) return kTRUE;
789c2e43 299 if (fOption==kResonanceOnly && fMFTTrackPair->GetCharge() != 0) return kTRUE;
300
301 Double_t pca[3] = {0};
302 fMFTTrackPair -> GetPointOfClosestApproach(pca);
303 Double_t distancePrimaryVtxPCA = TMath::Sqrt(TMath::Power(fPrimaryVtxX-pca[0],2)+
304 TMath::Power(fPrimaryVtxY-pca[1],2)+
305 TMath::Power(fPrimaryVtxZ-pca[2],2));
306
d8c2cc3e 307 fMFTTrackPair -> SetKinem(fPrimaryVtxZ);
308
cabfe25b 309 Bool_t passedCut = kFALSE;
310 if (fUseBransonForCut) passedCut = PassedCutMuonPair(fMFTTrackPairWithBranson);
311 else passedCut = PassedCutMuonPair(fMFTTrackPairWithoutBranson);
312
313 if (!passedCut) return kTRUE;
cd2f51d2 314
cabfe25b 315 // --------------- Filling dimuon histograms --------------------------------
cd2f51d2 316
789c2e43 317 if (fMFTTrackPair->GetCharge() == 0) {
cabfe25b 318 if (fOption==kResonanceOnly) fHistMassMuonPairsMCVsPt -> Fill(fMFTTrackPair->GetMassMC(), fMFTTrackPair->GetPt());
319 fHistMassMuonPairsVsPt -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
320 fHistMassMuonPairsWithoutMFTVsPt -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
321 fHistWOffsetMuonPairsAtPrimaryVtxVsPt -> Fill(fMFTTrackPair->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
322 fHistWOffsetMuonPairsAtPCAVsPt -> Fill(fMFTTrackPair->GetWeightedOffsetAtPCA(), fMFTTrackPair->GetPt());
323 fHistDistancePrimaryVtxPCAVsPt -> Fill(distancePrimaryVtxPCA*1.e4, fMFTTrackPair->GetPt());
324 fHistPCAQualityVsPt -> Fill(fMFTTrackPair->GetPCAQuality(), fMFTTrackPair->GetPt());
325 if (fOption==kResonanceOnly) fHistPseudoProperDecayLengthVsPt->Fill(GetPseudoProperDecayLength(fMFTTrackPair, fTrueMass), fMFTTrackPair->GetPt());
326 if (fEvalDimuonVtxResolution) {
327 fHistDimuonVtxResolutionXVsPt->Fill(pca[0]*1.e4, fMFTTrackPair->GetPt());
328 fHistDimuonVtxResolutionYVsPt->Fill(pca[1]*1.e4, fMFTTrackPair->GetPt());
329 fHistDimuonVtxResolutionZVsPt->Fill(pca[2]*1.e4, fMFTTrackPair->GetPt());
330 }
3b0abc02 331 fHistRapidityPtMuonPairs -> Fill(fMFTTrackPair->GetRapidity(), fMFTTrackPair->GetPt());
789c2e43 332 }
333 else if (fMFTTrackPair->GetCharge() == -2) {
334 fHistMassMuonPairsVsPtLSm -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
335 fHistMassMuonPairsWithoutMFTVsPtLSm -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
336 }
337 else if (fMFTTrackPair->GetCharge() == 2) {
338 fHistMassMuonPairsVsPtLSp -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
339 fHistMassMuonPairsWithoutMFTVsPtLSp -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
340 }
341
d8c2cc3e 342 AliDebug(1, Form("mass = %f MC = %f", fMFTTrackPair->GetMass(), fMFTTrackPair->GetMassMC()));
343
344 fNPairsAnalyzedOfEventAfterCut++;
e806e863 345
346 return kTRUE;
347
348}
349
350//====================================================================================================================================================
351
352void AliMuonForwardTrackAnalysis::BuildMuonPairs() {
353
789c2e43 354 Int_t nMuonPairs = 0;
355
356 for (Int_t iTrack=0; iTrack<fMuonForwardTracksWithBranson->GetEntriesFast(); iTrack++) {
e806e863 357 for (Int_t jTrack=0; jTrack<iTrack; jTrack++) {
789c2e43 358
7e3dd1af 359 AliMuonForwardTrack *track0_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(iTrack);
360 AliMuonForwardTrack *track1_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(jTrack);
789c2e43 361
7e3dd1af 362 AliMuonForwardTrack *track0_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(iTrack);
363 AliMuonForwardTrack *track1_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(jTrack);
e806e863 364
cabfe25b 365 if (track0_WithBranson->GetMatchTrigger()<fTriggerLevel || track1_WithBranson->GetMatchTrigger()<fTriggerLevel) continue;
e806e863 366
789c2e43 367 new ((*fMuonForwardTrackPairsWithBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithBranson, track1_WithBranson);
368 new ((*fMuonForwardTrackPairsWithoutBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithoutBranson, track1_WithoutBranson);
369
370 AliMuonForwardTrackPair *trackPairWithBranson = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithBranson->At(nMuonPairs);
371 if (!(fOption==kResonanceOnly && !trackPairWithBranson->IsResonance())) nMuonPairs++;
372
e806e863 373 }
374 }
375
376}
377
378//====================================================================================================================================================
379
7e3dd1af 380Bool_t AliMuonForwardTrackAnalysis::PassedCutSingleMuon(AliMuonForwardTrack *track) {
381
382 AliMUONTrackParam *param = track->GetTrackParamAtMFTCluster(0);
383 AliMUONTrackExtrap::ExtrapToZCov(param, fPrimaryVtxZ);
384
385 if (track->Pt()<fPtMinSingleMuons) return kFALSE;
386
cabfe25b 387 Double_t offset = 1.e4*track->GetOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
388 Double_t chi2OverNdf = track->GetGlobalChi2() / Double_t(track->GetNMFTClusters()+track->GetNMUONClusters());
389 offset /= fMaxOffsetSingleMuons;
390 chi2OverNdf /= fMaxChi2SingleMuons;
391
392 if (fCorrelateCutOnOffsetChi2) {
393 if (TMath::Sqrt(offset*offset + chi2OverNdf*chi2OverNdf) > 1) return kFALSE;
394 }
395 else {
396 if (offset>1 || chi2OverNdf>1) return kFALSE;
7e3dd1af 397 }
398
399 return kTRUE;
400
401}
402
403//====================================================================================================================================================
404
405Bool_t AliMuonForwardTrackAnalysis::PassedCutMuonPair(AliMuonForwardTrackPair *pair) {
406
cabfe25b 407 if (!PassedCutSingleMuon(pair->GetTrack(0)) || !PassedCutSingleMuon(pair->GetTrack(1))) return kFALSE;
408
409 if (pair->GetMass()>fMassMax || pair->GetMass()<fMassMin) return kFALSE;
7e3dd1af 410
cabfe25b 411 if (pair->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ) > fMaxWOffsetMuonPairsAtPrimaryVtx) return kFALSE;
412 if (pair->GetWeightedOffsetAtPCA() > fMaxWOffsetMuonPairsAtPCA) return kFALSE;
413
414 Double_t pca[3] = {0};
415 pair -> GetPointOfClosestApproach(pca);
416 Double_t distancePrimaryVtxPCA = TMath::Sqrt(TMath::Power(fPrimaryVtxX-pca[0],2)+
417 TMath::Power(fPrimaryVtxY-pca[1],2)+
418 TMath::Power(fPrimaryVtxZ-pca[2],2));
419
420 if (distancePrimaryVtxPCA > fMaxDistancePrimaryVtxPCA) return kFALSE;
421 if (pair->GetPCAQuality() < fMinPCAQuality) return kFALSE;
422
423 return kTRUE;
424
425}
426
427//====================================================================================================================================================
428
429Double_t AliMuonForwardTrackAnalysis::GetPseudoProperDecayLength(AliMuonForwardTrackPair *pair, Double_t trueMass) {
430
431 Double_t pca[3] = {0};
432 pair -> GetPointOfClosestApproach(pca);
433
434 TVector2 vecVertexPCA(pca[0]-fPrimaryVtxX, pca[1]-fPrimaryVtxY);
435 TVector2 dimuonPt(pair->GetPx(), pair->GetPy());
436 dimuonPt.Unit();
437
438 Double_t l_xy = vecVertexPCA*dimuonPt;
439 Double_t pseudoProperDecayLength = (l_xy * trueMass / pair->GetPt()) * 1.e4 ; // in micron
440
441 return pseudoProperDecayLength;
442
7e3dd1af 443}
444
445//====================================================================================================================================================
446
e806e863 447void AliMuonForwardTrackAnalysis::Terminate(Char_t *outputFileName) {
448
e806e863 449 TFile *fileOut = new TFile(Form("%s/%s",fOutputDir.Data(),outputFileName), "recreate");
450
451 printf("Writing output objects to file %s\n", fileOut->GetName());
452
cabfe25b 453 fHistOffsetSingleMuonsX -> Write();
454 fHistOffsetSingleMuonsY -> Write();
455 fHistErrorSingleMuonsX -> Write();
456 fHistErrorSingleMuonsY -> Write();
457 fHistOffsetSingleMuons -> Write();
458 fHistWOffsetSingleMuons -> Write();
e806e863 459
cabfe25b 460 fHistSingleMuonsPtRapidity -> Write();
b5ab1ac4 461 fHistSingleMuonsOffsetChi2 -> Write();
d8c2cc3e 462 fHistZOriginSingleMuonsMC -> Write();
463 fHistZROriginSingleMuonsMC -> Write();
464
cabfe25b 465 fHistMassMuonPairsVsPt -> Write();
466 fHistMassMuonPairsWithoutMFTVsPt -> Write();
467 fHistMassMuonPairsVsPtLSp -> Write();
468 fHistMassMuonPairsWithoutMFTVsPtLSp -> Write();
469 fHistMassMuonPairsVsPtLSm -> Write();
470 fHistMassMuonPairsWithoutMFTVsPtLSm -> Write();
471 fHistMassMuonPairsMCVsPt -> Write();
472 fHistWOffsetMuonPairsAtPrimaryVtxVsPt -> Write();
473 fHistWOffsetMuonPairsAtPCAVsPt -> Write();
474 fHistDistancePrimaryVtxPCAVsPt -> Write();
475 fHistPCAQualityVsPt -> Write();
476 if (fOption==kResonanceOnly) fHistPseudoProperDecayLengthVsPt->Write();
477 if (fEvalDimuonVtxResolution) {
478 fHistDimuonVtxResolutionXVsPt -> Write();
479 fHistDimuonVtxResolutionYVsPt -> Write();
480 fHistDimuonVtxResolutionZVsPt -> Write();
cd2f51d2 481 }
d8c2cc3e 482 fHistRapidityPtMuonPairs -> Write();
e806e863 483
484 fileOut -> Close();
485
486}
487
488//====================================================================================================================================================
489
490void AliMuonForwardTrackAnalysis::BookHistos() {
491
cabfe25b 492 // -------------------- single muons
493
e806e863 494 fHistOffsetSingleMuonsX = new TH1D("fHistOffsetSingleMuonsX", "Offset for single muons along X", 200, -1000, 1000);
495 fHistOffsetSingleMuonsY = new TH1D("fHistOffsetSingleMuonsY", "Offset for single muons along Y", 200, -1000, 1000);
496 fHistErrorSingleMuonsX = new TH1D("fHistErrorSingleMuonsX", "Coordinate Error for single muons along X", 200, 0, 1000);
497 fHistErrorSingleMuonsY = new TH1D("fHistErrorSingleMuonsY", "Coordinate Error for single muons along Y", 200, 0, 1000);
b5ab1ac4 498 fHistOffsetSingleMuons = new TH1D("fHistOffsetSingleMuons", "Offset for single muons", 200, 0, 2000);
e806e863 499 fHistWOffsetSingleMuons = new TH1D("fHistWOffsetSingleMuons", "Weighted Offset for single muons", 300, 0, 15);
500
cabfe25b 501 fHistSingleMuonsPtRapidity = new TH2D("fHistSingleMuonsPtRapidity", "Phase Space for single muons", 100, -4.5, -2., 100, 0., 10.);
789c2e43 502 fHistSingleMuonsOffsetChi2 = new TH2D("fHistSingleMuonsOffsetChi2", "Offset vs #chi^{2}/ndf for single muons", 400, 0, 4000, 200, 0, 10);
d8c2cc3e 503 fHistZOriginSingleMuonsMC = new TH1D("fHistZOriginSingleMuonsMC", "Z origin for single muons (from MC)", 1000, 0., 500.);
504 fHistZROriginSingleMuonsMC = new TH2D("fHistZROriginSingleMuonsMC", "Z-R origin for single muons (from MC)", 1000, 0., 500., 1000, 0., 100.);
e806e863 505
506 fHistOffsetSingleMuonsX -> SetXTitle("Offset(X) [#mum]");
507 fHistOffsetSingleMuonsY -> SetXTitle("Offset(Y) [#mum]");
508 fHistErrorSingleMuonsX -> SetXTitle("Err. on track position at z_{vtx} (X) [#mum]");
509 fHistErrorSingleMuonsY -> SetXTitle("Err. on track position at z_{vtx} (Y) [#mum]");
510 fHistOffsetSingleMuons -> SetXTitle("Offset [#mum]");
511 fHistWOffsetSingleMuons -> SetXTitle("Weighted Offset");
512
cabfe25b 513 fHistSingleMuonsPtRapidity -> SetXTitle("y^{#mu}");
514 fHistSingleMuonsPtRapidity -> SetYTitle("p_{T}^{#mu} [GeV/c]");
b5ab1ac4 515 fHistSingleMuonsOffsetChi2 -> SetXTitle("Offset [#mum]");
516 fHistSingleMuonsOffsetChi2 -> SetYTitle("#chi^{2}/ndf");
d8c2cc3e 517
518 fHistZOriginSingleMuonsMC -> SetXTitle("Z [cm]");
519 fHistZROriginSingleMuonsMC -> SetXTitle("Z [cm]");
520 fHistZROriginSingleMuonsMC -> SetXTitle("R [cm]");
e806e863 521
522 fHistOffsetSingleMuonsX -> Sumw2();
523 fHistOffsetSingleMuonsY -> Sumw2();
524 fHistErrorSingleMuonsX -> Sumw2();
525 fHistErrorSingleMuonsY -> Sumw2();
526 fHistOffsetSingleMuons -> Sumw2();
527 fHistWOffsetSingleMuons -> Sumw2();
528
cabfe25b 529 fHistSingleMuonsPtRapidity -> Sumw2();
530 fHistSingleMuonsOffsetChi2 -> Sumw2();
d8c2cc3e 531 fHistZOriginSingleMuonsMC -> Sumw2();
532 fHistZROriginSingleMuonsMC -> Sumw2();
533
cabfe25b 534 // -------------------- muon pairs
535
536 Int_t nBinsPt = 20; Double_t ptMin = 0., ptMax = 10.; // dimuon pt
537
538 fHistMassMuonPairsVsPt = new TH2D("fHistMassMuonPairsVsPt", "Dimuon Mass (MUON+MFT) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax);
539 fHistMassMuonPairsWithoutMFTVsPt = new TH2D("fHistMassMuonPairsWithoutMFTVsPt", "Dimuon Mass (MUON only) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax);
540 fHistMassMuonPairsVsPtLSp = new TH2D("fHistMassMuonPairsVsPtLSp", "Dimuon Mass (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax);
541 fHistMassMuonPairsWithoutMFTVsPtLSp = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSp", "Dimuon Mass (MUON only) vs p_{T}^{#mu#mu} Like-Sign ++", 1000, 0, 10, nBinsPt, ptMin, ptMax);
542 fHistMassMuonPairsVsPtLSm = new TH2D("fHistMassMuonPairsVsPtLSm", "Dimuon Mass (MUON+MFT) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax);
543 fHistMassMuonPairsWithoutMFTVsPtLSm = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSm", "Dimuon Mass (MUON only) vs p_{T}^{#mu#mu} Like-Sign --", 1000, 0, 10, nBinsPt, ptMin, ptMax);
544 fHistMassMuonPairsMCVsPt = new TH2D("fHistMassMuonPairsMCVsPt", "Dimuon Mass (MC) vs p_{T}^{#mu#mu}", 1000, 0, 10, nBinsPt, ptMin, ptMax);
545 fHistWOffsetMuonPairsAtPrimaryVtxVsPt = new TH2D("fHistWOffsetMuonPairsAtPrimaryVtxVsPt", "Weighted Offset for Muon Pairs at Primary Vertex vs p_{T}^{#mu#mu}", 300, 0, 60, nBinsPt, ptMin, ptMax);
546 fHistWOffsetMuonPairsAtPCAVsPt = new TH2D("fHistWOffsetMuonPairsAtPCAVsPt", "Weighted Offset for Muon Pairs at PCA vs p_{T}^{#mu#mu}", 300, 0, 60, nBinsPt, ptMin, ptMax);
547 fHistDistancePrimaryVtxPCAVsPt = new TH2D("fHistDistancePrimaryVtxPCA_%d", "Distance between PCA and primary vertex vs p_{T}^{#mu#mu}", 1000, 0, 50000, nBinsPt, ptMin, ptMax);
548 fHistPCAQualityVsPt = new TH2D("fHistPCAQualityVsPt", "PCA Quality vs p_{T}^{#mu#mu}", 200, 0, 1, nBinsPt, ptMin, ptMax);
549 fHistPseudoProperDecayLengthVsPt = new TH2D("fHistPseudoProperDecayLengthVsPt", "Pseudo proper decay length vs p_{T}^{#mu#mu}", 1000, -5000, 5000, nBinsPt, ptMin, ptMax);
550 fHistDimuonVtxResolutionXVsPt = new TH2D("fHistDimuonVtxResolutionXVsPt", "Dimuon Vtx offset along X w.r.t. true Vtx vs p_{T}^{#mu#mu}", 2000, -1000, 1000, nBinsPt, ptMin, ptMax);
551 fHistDimuonVtxResolutionYVsPt = new TH2D("fHistDimuonVtxResolutionYVsPt", "Dimuon Vtx offset along Y w.r.t. true Vtx vs p_{T}^{#mu#mu}", 2000, -1000, 1000, nBinsPt, ptMin, ptMax);
552 fHistDimuonVtxResolutionZVsPt = new TH2D("fHistDimuonVtxResolutionZVsPt", "Dimuon Vtx offset along Z w.r.t. true Vtx vs p_{T}^{#mu#mu}", 2000, -1000, 1000, nBinsPt, ptMin, ptMax);
cd2f51d2 553
cabfe25b 554 fHistMassMuonPairsVsPt -> SetXTitle("Mass [GeV/c^{2}]");
555 fHistMassMuonPairsWithoutMFTVsPt -> SetXTitle("Mass [GeV/c^{2}]");
556 fHistMassMuonPairsVsPtLSp -> SetXTitle("Mass [GeV/c^{2}]");
557 fHistMassMuonPairsWithoutMFTVsPtLSp -> SetXTitle("Mass [GeV/c^{2}]");
558 fHistMassMuonPairsVsPtLSm -> SetXTitle("Mass [GeV/c^{2}]");
559 fHistMassMuonPairsWithoutMFTVsPtLSm -> SetXTitle("Mass [GeV/c^{2}]");
560 fHistMassMuonPairsMCVsPt -> SetXTitle("Mass [GeV/c^{2}]");
561 fHistWOffsetMuonPairsAtPrimaryVtxVsPt -> SetXTitle("Weighted Offset at Primary Vtx");
562 fHistWOffsetMuonPairsAtPCAVsPt -> SetXTitle("Weighted Offset at PCA");
563 fHistDistancePrimaryVtxPCAVsPt -> SetXTitle("PCA - Primary Vtx [#mum]");
564 fHistPCAQualityVsPt -> SetXTitle("PCA Quality");
565 fHistPseudoProperDecayLengthVsPt -> SetXTitle("l_{xy} [#mum]");
566 fHistDimuonVtxResolutionXVsPt -> SetXTitle("Offset [#mum]");
567 fHistDimuonVtxResolutionYVsPt -> SetXTitle("Offset [#mum]");
568 fHistDimuonVtxResolutionZVsPt -> SetXTitle("Offset [#mum]");
569
570 fHistMassMuonPairsVsPt -> SetYTitle("p_{T} [GeV/c]");
571 fHistMassMuonPairsWithoutMFTVsPt -> SetYTitle("p_{T} [GeV/c]");
572 fHistMassMuonPairsVsPtLSp -> SetYTitle("p_{T} [GeV/c]");
573 fHistMassMuonPairsWithoutMFTVsPtLSp -> SetYTitle("p_{T} [GeV/c]");
574 fHistMassMuonPairsVsPtLSm -> SetYTitle("p_{T} [GeV/c]");
575 fHistMassMuonPairsWithoutMFTVsPtLSm -> SetYTitle("p_{T} [GeV/c]");
576 fHistMassMuonPairsMCVsPt -> SetYTitle("p_{T} [GeV/c]");
577 fHistWOffsetMuonPairsAtPrimaryVtxVsPt -> SetYTitle("p_{T} [GeV/c]");
578 fHistWOffsetMuonPairsAtPCAVsPt -> SetYTitle("p_{T} [GeV/c]");
579 fHistDistancePrimaryVtxPCAVsPt -> SetYTitle("p_{T} [GeV/c]");
580 fHistPCAQualityVsPt -> SetYTitle("p_{T} [GeV/c]");
581 fHistPseudoProperDecayLengthVsPt -> SetYTitle("p_{T} [GeV/c]");
582 fHistDimuonVtxResolutionXVsPt -> SetYTitle("p_{T} [GeV/c]");
583 fHistDimuonVtxResolutionYVsPt -> SetYTitle("p_{T} [GeV/c]");
584 fHistDimuonVtxResolutionZVsPt -> SetYTitle("p_{T} [GeV/c]");
585
586 fHistMassMuonPairsVsPt -> Sumw2();
587 fHistMassMuonPairsWithoutMFTVsPt -> Sumw2();
588 fHistMassMuonPairsVsPtLSp -> Sumw2();
589 fHistMassMuonPairsWithoutMFTVsPtLSp -> Sumw2();
590 fHistMassMuonPairsVsPtLSm -> Sumw2();
591 fHistMassMuonPairsWithoutMFTVsPtLSm -> Sumw2();
592 fHistMassMuonPairsMCVsPt -> Sumw2();
593 fHistWOffsetMuonPairsAtPrimaryVtxVsPt -> Sumw2();
594 fHistWOffsetMuonPairsAtPCAVsPt -> Sumw2();
595 fHistDistancePrimaryVtxPCAVsPt -> Sumw2();
596 fHistPCAQualityVsPt -> Sumw2();
597 fHistPseudoProperDecayLengthVsPt -> Sumw2();
598 fHistDimuonVtxResolutionXVsPt -> Sumw2();
599 fHistDimuonVtxResolutionYVsPt -> Sumw2();
600 fHistDimuonVtxResolutionZVsPt -> Sumw2();
601
602 fHistRapidityPtMuonPairs = new TH2D("fHistRapidityPtMuonPairs", "Dimuon Phase Space (rec)", 20, -4.5, -2., 20, 0., 10.);
603 fHistRapidityPtMuonPairs -> SetXTitle("Rapidity");
604 fHistRapidityPtMuonPairs -> SetYTitle("p_{T} [GeV/c]");
605 fHistRapidityPtMuonPairs -> Sumw2();
e806e863 606
607}
608
609//====================================================================================================================================================
610