MFT/AliMuonForwardTrackAnalysis.cxx

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