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   fHistSingleMuonsPtRapidityMC(0x0),
  68   fHistSingleMuonsOffsetChi2(0x0),
  69   fHistRapidityPtMuonPairs(0x0),
  70   fHistMassMuonPairsVsPt(0),
  71   fHistMassMuonPairsWithoutMFTVsPt(0),
  72   fHistMassMuonPairsVsPtLSp(0),
  73   fHistMassMuonPairsWithoutMFTVsPtLSp(0),
  74   fHistMassMuonPairsVsPtLSm(0),
  75   fHistMassMuonPairsWithoutMFTVsPtLSm(0),
  76   fEvalDimuonVtxResolution(kFALSE),
  77   fNMassBins(10),
  78   fNPtDimuBins(1000),
  79   fMassMin(0),
  80   fMassMax(10),
  81   fPtDimuMin(0),
  82   fPtDimuMax(5),
  83   fPtAxisDimuons(0),
  84   fSingleMuonAnalysis(1),
  85   fMuonPairAnalysis(1),
  86   fMatchTrigger(0),
  87   fOption(0),
  88   fXVertResMC(50.e-4),
  89   fYVertResMC(50.e-4),
  90   fZVertResMC(50.e-4),
  91   fPrimaryVtxX(0.),
  92   fPrimaryVtxY(0.),
  93   fPrimaryVtxZ(0.),
  94   fMaxNWrongClustersMC(999),
  95   fPtMinSingleMuons(0),
  96   fUseBransonForCut(kFALSE),
  97   fUseBransonForKinematics(kFALSE),
  98   fCutOnOffsetChi2(kFALSE),
  99   fCenterOffset(0.),
 100   fCenterChi2(0.),
 101   fScaleOffset(250.),
 102   fScaleChi2(1.5),
 103   fRadiusCut(1.)
 104 {
 105
 106   // default constructor
 107
 108   for (Int_t iPtBin=0; iPtBin<fNMaxPtBinsDimuons+1; iPtBin++) {
 109     fHistMassMuonPairs[iPtBin]                                = NULL;
 110     fHistMassMuonPairsWithoutMFT[iPtBin]                      = NULL;
 111     fHistMassMuonPairsMC[iPtBin]                              = NULL;
 112     fHistWOffsetMuonPairsAtPrimaryVtx[iPtBin]                 = NULL;
 113     fHistWOffsetMuonPairsAtPCA[iPtBin]                        = NULL;
 114     fHistDistancePrimaryVtxPCA[iPtBin]                        = NULL;
 115     fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA[iPtBin] = NULL;
 116     fHistDimuonVtxResolutionX[iPtBin]                         = NULL;
 117     fHistDimuonVtxResolutionY[iPtBin]                         = NULL;
 118     fHistDimuonVtxResolutionZ[iPtBin]                         = NULL;
 119   }
 120
 121 }
 122
 123 //====================================================================================================================================================
 124
 125 Bool_t AliMuonForwardTrackAnalysis::Init(Char_t *inputFileName) {
 126
 127   fPtAxisDimuons = new TAxis(fNPtDimuBins, fPtDimuMin, fPtDimuMax);
 128
 129   BookHistos();
 130
 131   TFile *inputFile = new TFile(Form("%s/%s",fInputDir.Data(),inputFileName));
 132   if (!inputFile || !inputFile->IsOpen()) {
 133     AliError(Form("Error opening file %s", inputFileName));
 134     return kFALSE;
 135   }
 136   fInputTreeWithBranson = (TTree*) inputFile->Get("AliMuonForwardTracksWithBranson");
 137   if (!fInputTreeWithBranson) {
 138     AliError("Error reading input tree");
 139     return kFALSE;
 140   }
 141   fInputTreeWithoutBranson = (TTree*) inputFile->Get("AliMuonForwardTracksWithoutBranson");
 142   if (!fInputTreeWithoutBranson) {
 143     AliError("Error reading input tree");
 144     return kFALSE;
 145   }
 146
 147   if (fFirstEvent<0 || fLastEvent<0 || fFirstEvent>fLastEvent || fFirstEvent>=fInputTreeWithBranson->GetEntriesFast()) {
 148     fFirstEvent = 0;
 149     fLastEvent  = fInputTreeWithBranson->GetEntriesFast()-1;
 150   }
 151   else {
 152     fLastEvent = TMath::Min(fLastEvent, Int_t(fInputTreeWithBranson->GetEntriesFast()-1));
 153   }
 154
 155   AliInfo(Form("Analysing events %d to %d", fFirstEvent, fLastEvent));
 156
 157   fMuonForwardTracksWithBranson = new TClonesArray("AliMuonForwardTrack",30);
 158   fInputTreeWithBranson->SetBranchAddress("tracks", &fMuonForwardTracksWithBranson);
 159
 160   fMuonForwardTracksWithoutBranson = new TClonesArray("AliMuonForwardTrack",30);
 161   fInputTreeWithoutBranson->SetBranchAddress("tracks", &fMuonForwardTracksWithoutBranson);
 162
 163   TGeoManager::Import(Form("%s/geometry.root",fInputDir.Data()));
 164
 165   AliMUONTrackExtrap::SetField();
 166
 167   fMuonForwardTrackPairsWithBranson    = new TClonesArray("AliMuonForwardTrackPair",10);
 168   fMuonForwardTrackPairsWithoutBranson = new TClonesArray("AliMuonForwardTrackPair",10);
 169         fMuonForwardTrackPairsWithBranson->SetOwner(kTRUE);
 170         fMuonForwardTrackPairsWithoutBranson->SetOwner(kTRUE);
 171
 172
 173   return kTRUE;
 174
 175 }
 176
 177 //====================================================================================================================================================
 178
 179 Bool_t AliMuonForwardTrackAnalysis::LoadNextEvent() {
 180
 181   if (fEv>fLastEvent) return kFALSE;
 182   if (fEv<fFirstEvent) { fEv++; return kTRUE; }
 183   fMuonForwardTracksWithBranson -> Clear("");
 184   fMuonForwardTracksWithoutBranson -> Clear("");
 185   fInputTreeWithBranson->GetEvent(fEv);
 186   fInputTreeWithoutBranson->GetEvent(fEv);
 187   AliDebug(2,Form("**** analyzing event # %4d (%3d tracks) ****", fEv, fMuonForwardTracksWithBranson->GetEntriesFast()));
 188
 189   fPrimaryVtxX = gRandom->Gaus(0., fXVertResMC);
 190   fPrimaryVtxY = gRandom->Gaus(0., fYVertResMC);
 191   fPrimaryVtxZ = gRandom->Gaus(0., fZVertResMC);
 192
 193   if (fSingleMuonAnalysis) {
 194     fNTracksAnalyzedOfEvent = 0;
 195     fNTracksAnalyzedOfEventAfterCut = 0;
 196     fNTracksOfEvent = fMuonForwardTracksWithBranson->GetEntriesFast();
 197     while (AnalyzeSingleMuon()) continue;
 198   }
 199
 200   if (fMuonPairAnalysis) {
 201     if (fMuonForwardTrackPairsWithBranson) {
 202       fMuonForwardTrackPairsWithBranson->Clear("C");
 203       fMuonForwardTrackPairsWithoutBranson->Clear("C");
 204     }
 205     BuildMuonPairs();
 206     fNPairsAnalyzedOfEvent = 0;
 207     fNPairsAnalyzedOfEventAfterCut = 0;
 208    // fNPairsOfEvent = fMuonForwardTrackPairsWithBranson->GetEntriesFast();
 209     while (AnalyzeMuonPair()) continue;
 210   }
 211
 212   AliDebug(2,Form("**** analyzed  event # %4d (%3d tracks and %3d pairs analyzed) ****", fEv, fNTracksAnalyzedOfEventAfterCut, fNPairsAnalyzedOfEventAfterCut));
 213
 214   fEv++;
 215
 216   return kTRUE;
 217
 218 }
 219
 220 //====================================================================================================================================================
 221
 222 Bool_t AliMuonForwardTrackAnalysis::AnalyzeSingleMuon() {
 223
 224   if (fNTracksAnalyzedOfEvent>=fNTracksOfEvent) return kFALSE;
 225
 226   fMFTTrackWithBranson    = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(fNTracksAnalyzedOfEvent);
 227   fMFTTrackWithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(fNTracksAnalyzedOfEvent);
 228
 229   fNTracksAnalyzedOfEvent++;
 230
 231   Bool_t passedCut = kFALSE;
 232   if (fUseBransonForCut) passedCut = PassedCutSingleMuon(fMFTTrackWithBranson);
 233   else passedCut = PassedCutSingleMuon(fMFTTrackWithoutBranson);
 234   if (!passedCut) return kTRUE;
 235
 236   if (fUseBransonForKinematics) fMFTTrack = fMFTTrackWithBranson;
 237   else fMFTTrack = fMFTTrackWithoutBranson;
 238   if (fMatchTrigger && !fMFTTrack->GetMatchTrigger()) return kTRUE;
 239   fMCRefTrack = fMFTTrack->GetMCTrackRef();
 240
 241   if (fOption!=kPionsKaons && !fMCRefTrack) return kTRUE;
 242   if (fOption!=kPionsKaons && fMFTTrack->GetNWrongClustersMC()>fMaxNWrongClustersMC) return kTRUE;
 243
 244   if (fMCRefTrack) {
 245     fHistZOriginSingleMuonsMC  -> Fill(-1.*fMCRefTrack->Vz());
 246     fHistZROriginSingleMuonsMC -> Fill(-1.*fMCRefTrack->Vz(), TMath::Sqrt(fMCRefTrack->Vx()*fMCRefTrack->Vx()+fMCRefTrack->Vy()*fMCRefTrack->Vy()));
 247   }
 248
 249   AliMUONTrackParam *param = fMFTTrack->GetTrackParamAtMFTCluster(0);
 250   AliMUONTrackExtrap::ExtrapToZCov(param, fPrimaryVtxZ);
 251
 252   TLorentzVector pMu;
 253   Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
 254   Double_t energy = TMath::Sqrt(param->P()*param->P() + mMu*mMu);
 255   pMu.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
 256
 257   TMatrixD cov(5,5);
 258   cov = param->GetCovariances();
 259
 260   fHistErrorSingleMuonsX -> Fill(1.e4*TMath::Sqrt(cov(0,0)));
 261   fHistErrorSingleMuonsY -> Fill(1.e4*TMath::Sqrt(cov(2,2)));
 262
 263   Double_t dX = fMFTTrack->GetOffsetX(fPrimaryVtxX, fPrimaryVtxZ);
 264   Double_t dY = fMFTTrack->GetOffsetY(fPrimaryVtxY, fPrimaryVtxZ);
 265
 266   Double_t offset = fMFTTrack->GetOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
 267   Double_t weightedOffset = fMFTTrack->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
 268
 269   //  AliDebug(2, Form("pdg code = %d\n", fMCRefTrack->GetPdgCode()));
 270
 271   fHistOffsetSingleMuonsX -> Fill(1.e4*dX);
 272   fHistOffsetSingleMuonsY -> Fill(1.e4*dY);
 273
 274   if (fOption!=kPionsKaons) fHistSingleMuonsPtRapidityMC-> Fill(fMCRefTrack->Y(), fMCRefTrack->Pt());
 275   fHistOffsetSingleMuons      -> Fill(1.e4*offset);
 276   fHistWOffsetSingleMuons     -> Fill(weightedOffset);
 277   Double_t chi2OverNdf = fMFTTrack->GetGlobalChi2()/Double_t(fMFTTrack->GetNMFTClusters()+fMFTTrack->GetNMUONClusters());
 278   fHistSingleMuonsOffsetChi2  -> Fill(1.e4*offset, chi2OverNdf);
 279
 280   fNTracksAnalyzed++;
 281   fNTracksAnalyzedOfEventAfterCut++;
 282
 283   return kTRUE;
 284
 285 }
 286
 287 //====================================================================================================================================================
 288
 289 Bool_t AliMuonForwardTrackAnalysis::AnalyzeMuonPair() {
 290
 291   if (fNPairsAnalyzedOfEvent>=fNPairsOfEvent) return kFALSE;
 292
 293   fMFTTrackPairWithBranson    = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithBranson->At(fNPairsAnalyzedOfEvent);
 294   fMFTTrackPairWithoutBranson = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithoutBranson->At(fNPairsAnalyzedOfEvent);
 295
 296   fNPairsAnalyzedOfEvent++;
 297
 298   Bool_t passedCut = kFALSE;
 299   if (fUseBransonForCut) passedCut = PassedCutMuonPair(fMFTTrackPairWithBranson);
 300   else passedCut = PassedCutMuonPair(fMFTTrackPairWithoutBranson);
 301
 302   if (!passedCut) return kTRUE;
 303
 304   if (fUseBransonForKinematics) fMFTTrackPair = fMFTTrackPairWithBranson;
 305   else fMFTTrackPair = fMFTTrackPairWithoutBranson;
 306
 307   if ( fMFTTrackPair->GetTrack(0)->GetNWrongClustersMC()>fMaxNWrongClustersMC ||
 308        fMFTTrackPair->GetTrack(1)->GetNWrongClustersMC()>fMaxNWrongClustersMC ) return kTRUE;
 309
 310   if (fOption==kResonanceOnly && !fMFTTrackPair->IsResonance()) return kTRUE;
 311
 312   if (fOption==kResonanceOnly && fMFTTrackPair->GetCharge() != 0) return kTRUE;
 313
 314   Double_t pca[3] = {0};
 315   fMFTTrackPair -> GetPointOfClosestApproach(pca);
 316   Double_t distancePrimaryVtxPCA = TMath::Sqrt(TMath::Power(fPrimaryVtxX-pca[0],2)+
 317                                                TMath::Power(fPrimaryVtxY-pca[1],2)+
 318                                                TMath::Power(fPrimaryVtxZ-pca[2],2));
 319
 320   fMFTTrackPair -> SetKinem(fPrimaryVtxZ);
 321
 322   Int_t ptBin = fPtAxisDimuons->FindBin(fMFTTrackPair->GetPt());
 323
 324   if (1<=ptBin && ptBin<=fNPtDimuBins) {
 325     fHistMassMuonPairs[ptBin]           -> Fill(fMFTTrackPair->GetMass());
 326     fHistMassMuonPairsWithoutMFT[ptBin] -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ));
 327     if (fOption!=kPionsKaons) fHistMassMuonPairsMC[ptBin] -> Fill(fMFTTrackPair->GetMassMC());
 328     fHistWOffsetMuonPairsAtPrimaryVtx[ptBin] -> Fill(fMFTTrackPair->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ));
 329     fHistWOffsetMuonPairsAtPCA[ptBin]        -> Fill(fMFTTrackPair->GetWeightedOffset(pca[0], pca[1], pca[2]));
 330     fHistDistancePrimaryVtxPCA[ptBin]        -> Fill(distancePrimaryVtxPCA*1.e4);
 331     fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA[ptBin] -> Fill(fMFTTrackPair->GetWeightedOffset(pca[0], pca[1], pca[2]), distancePrimaryVtxPCA*1.e4);
 332     if (fEvalDimuonVtxResolution) {
 333       fHistDimuonVtxResolutionX[ptBin]->Fill(pca[0]*1.e4);
 334       fHistDimuonVtxResolutionY[ptBin]->Fill(pca[1]*1.e4);
 335       fHistDimuonVtxResolutionZ[ptBin]->Fill(pca[2]*1.e4);
 336     }
 337   }
 338   fHistMassMuonPairs[0]           -> Fill(fMFTTrackPair->GetMass());
 339   fHistMassMuonPairsWithoutMFT[0] -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ));
 340   if (fOption!=kPionsKaons) fHistMassMuonPairsMC[0] -> Fill(fMFTTrackPair->GetMassMC());
 341   fHistWOffsetMuonPairsAtPrimaryVtx[0] -> Fill(fMFTTrackPair->GetWeightedOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ));
 342   fHistWOffsetMuonPairsAtPCA[0]        -> Fill(fMFTTrackPair->GetWeightedOffset(pca[0], pca[1], pca[2]));
 343   fHistDistancePrimaryVtxPCA[0]        -> Fill(distancePrimaryVtxPCA*1.e4);
 344   fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA[0] -> Fill(fMFTTrackPair->GetWeightedOffset(pca[0], pca[1], pca[2]), distancePrimaryVtxPCA*1.e4);
 345   if (fEvalDimuonVtxResolution) {
 346     fHistDimuonVtxResolutionX[0]->Fill(pca[0]*1.e4);
 347     fHistDimuonVtxResolutionY[0]->Fill(pca[1]*1.e4);
 348     fHistDimuonVtxResolutionZ[0]->Fill(pca[2]*1.e4);
 349   }
 350
 351   if (fMFTTrackPair->GetCharge() == 0) {
 352     fHistMassMuonPairsVsPt           -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
 353     fHistMassMuonPairsWithoutMFTVsPt -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
 354   }
 355   else if (fMFTTrackPair->GetCharge() == -2) {
 356     fHistMassMuonPairsVsPtLSm           -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
 357     fHistMassMuonPairsWithoutMFTVsPtLSm -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
 358   }
 359   else if (fMFTTrackPair->GetCharge() == 2) {
 360     fHistMassMuonPairsVsPtLSp           -> Fill(fMFTTrackPair->GetMass(), fMFTTrackPair->GetPt());
 361     fHistMassMuonPairsWithoutMFTVsPtLSp -> Fill(fMFTTrackPair->GetMassWithoutMFT(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ), fMFTTrackPair->GetPt());
 362   }
 363
 364   if (fOption!=kPionsKaons) fHistRapidityPtMuonPairs -> Fill(fMFTTrackPair->GetRapidityMC(), fMFTTrackPair->GetPt());
 365   else fHistRapidityPtMuonPairs -> Fill(fMFTTrackPair->GetRapidity(), fMFTTrackPair->GetPt());
 366
 367   AliDebug(1, Form("mass = %f   MC = %f", fMFTTrackPair->GetMass(), fMFTTrackPair->GetMassMC()));
 368
 369   fNPairsAnalyzedOfEventAfterCut++;
 370
 371   return kTRUE;
 372
 373 }
 374
 375 //====================================================================================================================================================
 376
 377 void AliMuonForwardTrackAnalysis::BuildMuonPairs() {
 378
 379   Int_t nMuonPairs = 0;
 380
 381   for (Int_t iTrack=0; iTrack<fMuonForwardTracksWithBranson->GetEntriesFast(); iTrack++) {
 382     for (Int_t jTrack=0; jTrack<iTrack; jTrack++) {
 383
 384       AliMuonForwardTrack *track0_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(iTrack);
 385       AliMuonForwardTrack *track1_WithBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithBranson->At(jTrack);
 386
 387       AliMuonForwardTrack *track0_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(iTrack);
 388       AliMuonForwardTrack *track1_WithoutBranson = (AliMuonForwardTrack*) fMuonForwardTracksWithoutBranson->At(jTrack);
 389
 390       if (fMatchTrigger) if (!track0_WithBranson->GetMatchTrigger() || !track1_WithBranson->GetMatchTrigger()) continue;
 391       if (fOption!=kPionsKaons && (!track0_WithBranson->GetMCTrackRef() || !track1_WithBranson->GetMCTrackRef())) continue;
 392
 393       new ((*fMuonForwardTrackPairsWithBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithBranson, track1_WithBranson);
 394       new ((*fMuonForwardTrackPairsWithoutBranson)[nMuonPairs]) AliMuonForwardTrackPair(track0_WithoutBranson, track1_WithoutBranson);
 395
 396       AliMuonForwardTrackPair *trackPairWithBranson = (AliMuonForwardTrackPair*) fMuonForwardTrackPairsWithBranson->At(nMuonPairs);
 397       if (!(fOption==kResonanceOnly && !trackPairWithBranson->IsResonance())) nMuonPairs++;
 398
 399     }
 400   }
 401
 402 }
 403
 404 //====================================================================================================================================================
 405
 406 Bool_t AliMuonForwardTrackAnalysis::PassedCutSingleMuon(AliMuonForwardTrack *track) {
 407
 408   AliMUONTrackParam *param = track->GetTrackParamAtMFTCluster(0);
 409   AliMUONTrackExtrap::ExtrapToZCov(param, fPrimaryVtxZ);
 410
 411   if (track->Pt()<fPtMinSingleMuons) return kFALSE;
 412
 413   if (fCutOnOffsetChi2) {
 414     Double_t offset = 1.e4*track->GetOffset(fPrimaryVtxX, fPrimaryVtxY, fPrimaryVtxZ);
 415     Double_t chi2OverNdf = track->GetGlobalChi2() / Double_t(track->GetNMFTClusters()+track->GetNMUONClusters());
 416     offset /= fScaleOffset;
 417     chi2OverNdf /= fScaleChi2;
 418     offset -= fCenterOffset;
 419     chi2OverNdf -= fCenterChi2;
 420     //    printf("cut on offset and chi2: returning %d\n", TMath::Sqrt(offset*offset + chi2OverNdf*chi2OverNdf)>fRadiusCut);
 421     if (TMath::Sqrt(offset*offset + chi2OverNdf*chi2OverNdf) > fRadiusCut) return kFALSE;
 422   }
 423
 424   return kTRUE;
 425
 426 }
 427
 428 //====================================================================================================================================================
 429
 430 Bool_t AliMuonForwardTrackAnalysis::PassedCutMuonPair(AliMuonForwardTrackPair *pair) {
 431
 432   return PassedCutSingleMuon(pair->GetTrack(0)) && PassedCutSingleMuon(pair->GetTrack(1));
 433
 434 }
 435
 436 //====================================================================================================================================================
 437
 438 void AliMuonForwardTrackAnalysis::Terminate(Char_t *outputFileName) {
 439
 440   TFile *fileOut = new TFile(Form("%s/%s",fOutputDir.Data(),outputFileName), "recreate");
 441
 442   printf("Writing output objects to file %s\n", fileOut->GetName());
 443
 444   fHistOffsetSingleMuonsX -> Write();
 445   fHistOffsetSingleMuonsY -> Write();
 446   fHistOffsetSingleMuons  -> Write();
 447   fHistWOffsetSingleMuons -> Write();
 448   fHistErrorSingleMuonsX  -> Write();
 449   fHistErrorSingleMuonsY  -> Write();
 450
 451   fHistSingleMuonsPtRapidityMC -> Write();
 452   fHistSingleMuonsOffsetChi2 -> Write();
 453
 454   fHistZOriginSingleMuonsMC  -> Write();
 455   fHistZROriginSingleMuonsMC -> Write();
 456
 457   for (Int_t iPtBin=0; iPtBin<fNPtDimuBins+1; iPtBin++) {
 458     fHistMassMuonPairs[iPtBin]                                -> Write();
 459     fHistMassMuonPairsWithoutMFT[iPtBin]                      -> Write();
 460     fHistMassMuonPairsMC[iPtBin]                              -> Write();
 461     fHistWOffsetMuonPairsAtPrimaryVtx[iPtBin]                 -> Write();
 462     fHistWOffsetMuonPairsAtPCA[iPtBin]                        -> Write();
 463     fHistDistancePrimaryVtxPCA[iPtBin]                        -> Write();
 464     fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA[iPtBin] -> Write();
 465     if (fEvalDimuonVtxResolution) {
 466       fHistDimuonVtxResolutionX[iPtBin] -> Write();
 467       fHistDimuonVtxResolutionY[iPtBin] -> Write();
 468       fHistDimuonVtxResolutionZ[iPtBin] -> Write();
 469     }
 470   }
 471
 472   fHistMassMuonPairsVsPt           -> Write();
 473   fHistMassMuonPairsWithoutMFTVsPt -> Write();
 474   fHistMassMuonPairsVsPtLSp           -> Write();
 475   fHistMassMuonPairsWithoutMFTVsPtLSp -> Write();
 476   fHistMassMuonPairsVsPtLSm           -> Write();
 477   fHistMassMuonPairsWithoutMFTVsPtLSm -> Write();
 478
 479   fHistRapidityPtMuonPairs -> Write();
 480
 481   fileOut -> Close();
 482
 483 }
 484
 485 //====================================================================================================================================================
 486
 487 void AliMuonForwardTrackAnalysis::BookHistos() {
 488
 489   fHistOffsetSingleMuonsX = new TH1D("fHistOffsetSingleMuonsX", "Offset for single muons along X",  200, -1000, 1000);
 490   fHistOffsetSingleMuonsY = new TH1D("fHistOffsetSingleMuonsY", "Offset for single muons along Y",  200, -1000, 1000);
 491   fHistErrorSingleMuonsX  = new TH1D("fHistErrorSingleMuonsX",  "Coordinate Error for single muons along X",  200, 0, 1000);
 492   fHistErrorSingleMuonsY  = new TH1D("fHistErrorSingleMuonsY",  "Coordinate Error for single muons along Y",  200, 0, 1000);
 493   fHistOffsetSingleMuons  = new TH1D("fHistOffsetSingleMuons",  "Offset for single muons",          200, 0, 2000);
 494   fHistWOffsetSingleMuons = new TH1D("fHistWOffsetSingleMuons", "Weighted Offset for single muons", 300, 0, 15);
 495
 496   fHistSingleMuonsPtRapidityMC = new TH2D("fHistSingleMuonsPtRapidityMC", "Phase Space for single muons", 100, -4.5, -2., 100, 0., 10.);
 497   fHistSingleMuonsOffsetChi2 = new TH2D("fHistSingleMuonsOffsetChi2", "Offset vs #chi^{2}/ndf for single muons", 400, 0, 4000, 200, 0, 10);
 498   fHistZOriginSingleMuonsMC  = new TH1D("fHistZOriginSingleMuonsMC",  "Z origin for single muons (from MC)",   1000, 0., 500.);
 499   fHistZROriginSingleMuonsMC = new TH2D("fHistZROriginSingleMuonsMC", "Z-R origin for single muons (from MC)", 1000, 0., 500., 1000, 0., 100.);
 500
 501   fHistOffsetSingleMuonsX -> SetXTitle("Offset(X)  [#mum]");
 502   fHistOffsetSingleMuonsY -> SetXTitle("Offset(Y)  [#mum]");
 503   fHistErrorSingleMuonsX  -> SetXTitle("Err. on track position at z_{vtx} (X)  [#mum]");
 504   fHistErrorSingleMuonsY  -> SetXTitle("Err. on track position at z_{vtx} (Y)  [#mum]");
 505   fHistOffsetSingleMuons  -> SetXTitle("Offset  [#mum]");
 506   fHistWOffsetSingleMuons -> SetXTitle("Weighted Offset");
 507
 508   fHistSingleMuonsPtRapidityMC -> SetXTitle("y^{#mu}");
 509   fHistSingleMuonsPtRapidityMC -> SetYTitle("p_{T}^{#mu}  [GeV/c]");
 510   fHistSingleMuonsOffsetChi2 -> SetXTitle("Offset  [#mum]");
 511   fHistSingleMuonsOffsetChi2 -> SetYTitle("#chi^{2}/ndf");
 512
 513   fHistZOriginSingleMuonsMC  -> SetXTitle("Z  [cm]");
 514   fHistZROriginSingleMuonsMC -> SetXTitle("Z  [cm]");
 515   fHistZROriginSingleMuonsMC -> SetXTitle("R  [cm]");
 516
 517   fHistOffsetSingleMuonsX -> Sumw2();
 518   fHistOffsetSingleMuonsY -> Sumw2();
 519   fHistErrorSingleMuonsX  -> Sumw2();
 520   fHistErrorSingleMuonsY  -> Sumw2();
 521   fHistOffsetSingleMuons  -> Sumw2();
 522   fHistWOffsetSingleMuons -> Sumw2();
 523
 524   fHistZOriginSingleMuonsMC  -> Sumw2();
 525   fHistZROriginSingleMuonsMC -> Sumw2();
 526
 527   fHistSingleMuonsPtRapidityMC -> Sumw2();
 528   fHistSingleMuonsOffsetChi2 -> Sumw2();
 529
 530   //--------------------------------------------
 531
 532   for (Int_t iPtBin=0; iPtBin<=fNPtDimuBins+1; iPtBin++) {
 533
 534     if (!iPtBin) {
 535       fHistMassMuonPairs[iPtBin]           = new TH1D(Form("fHistMassMuonPairs_%d",iPtBin),
 536                                                       "Dimuon Mass (MUON+MFT) (All p_{T}^{#mu#mu})",
 537                                                       fNMassBins, fMassMin, fMassMax);
 538       fHistMassMuonPairsWithoutMFT[iPtBin] = new TH1D(Form("fHistMassMuonPairsWithoutMFT_%d",iPtBin),
 539                                                       "Dimuon Mass (MUON only) (All p_{T}^{#mu#mu})",
 540                                                       fNMassBins, fMassMin, fMassMax);
 541       fHistMassMuonPairsMC[iPtBin]         = new TH1D(Form("fHistMassMuonPairsMC_%d",iPtBin),
 542                                                       "Dimuon Mass (MC) (All p_{T}^{#mu#mu})",
 543                                                       fNMassBins, fMassMin, fMassMax);
 544       fHistWOffsetMuonPairsAtPrimaryVtx[iPtBin] = new TH1D(Form("fHistWOffsetMuonPairsAtPrimaryVtx_%d",iPtBin),
 545                                                            "Weighted Offset for Muon Pairs at Primary Vertex (All p_{T}^{#mu#mu})",
 546                                                            300, 0, 60);
 547       fHistWOffsetMuonPairsAtPCA[iPtBin]        = new TH1D(Form("fHistWOffsetMuonPairsAtPCA_%d",iPtBin),
 548                                                            "Weighted Offset for Muon Pairs at PCA (All p_{T}^{#mu#mu})",
 549                                                            300, 0, 60);
 550       fHistDistancePrimaryVtxPCA[iPtBin]        = new TH1D(Form("fHistDistancePrimaryVtxPCA_%d",iPtBin),
 551                                                            "Distance between PCA and primary vertex (All p_{T}^{#mu#mu})",
 552                                                            1000, 0, 50000);
 553       fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA[iPtBin] = new TH2D(Form("fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA_%d",iPtBin),
 554                                                                            "DistancePrimaryVtxPCA vs WOffsetMuonPairsAtPCA (All p_{T}^{#mu#mu})",
 555                                                                            300, 0, 60, 1000, 0, 50000);
 556       if (fEvalDimuonVtxResolution) {
 557         fHistDimuonVtxResolutionX[iPtBin] = new TH1D(Form("fHistDimuonVtxResolutionX_%d",iPtBin),
 558                                                      "Dimuon Vtx offset along X (All p_{T}^{#mu#mu})",
 559                                                      2000, -1000, 1000);
 560         fHistDimuonVtxResolutionY[iPtBin] = new TH1D(Form("fHistDimuonVtxResolutionY_%d",iPtBin),
 561                                                      "Dimuon Vtx offset along Y (All p_{T}^{#mu#mu})",
 562                                                      2000, -1000, 1000);
 563         fHistDimuonVtxResolutionZ[iPtBin] = new TH1D(Form("fHistDimuonVtxResolutionZ_%d",iPtBin),
 564                                                      "Dimuon Vtx offset along Z (All p_{T}^{#mu#mu})",
 565                                                      2000, -10000, 10000);
 566       }
 567     }
 568
 569     else {
 570       Double_t ptMin = fPtAxisDimuons->GetBinLowEdge(iPtBin);
 571       Double_t ptMax = fPtAxisDimuons->GetBinUpEdge(iPtBin);
 572       fHistMassMuonPairs[iPtBin]           = new TH1D(Form("fHistMassMuonPairs_%d",iPtBin),
 573                                                       Form("Dimuon Mass (MUON+MFT) (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",ptMin,ptMax),
 574                                                       fNMassBins, fMassMin, fMassMax);
 575       fHistMassMuonPairsWithoutMFT[iPtBin] = new TH1D(Form("fHistMassMuonPairsWithoutMFT_%d",iPtBin),
 576                                                       Form("Dimuon Mass (MUON only) (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",ptMin,ptMax),
 577                                                       fNMassBins, fMassMin, fMassMax);
 578       fHistMassMuonPairsMC[iPtBin]         = new TH1D(Form("fHistMassMuonPairsMC_%d",iPtBin),
 579                                                       Form("Dimuon Mass (MC) (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",ptMin,ptMax),
 580                                                       fNMassBins, fMassMin, fMassMax);
 581       fHistWOffsetMuonPairsAtPrimaryVtx[iPtBin] = new TH1D(Form("fHistWOffsetMuonPairsAtPrimaryVtx_%d",iPtBin),
 582                                                            Form("Weighted Offset for Muon Pairs at Primary Vertex (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",
 583                                                                 ptMin,ptMax), 300, 0, 60);
 584       fHistWOffsetMuonPairsAtPCA[iPtBin]        = new TH1D(Form("fHistWOffsetMuonPairsAtPCA_%d",iPtBin),
 585                                                            Form("Weighted Offset for Muon Pairs at PCA (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",
 586                                                                 ptMin,ptMax), 300, 0, 60);
 587       fHistDistancePrimaryVtxPCA[iPtBin]        = new TH1D(Form("fHistDistancePrimaryVtxPCA_%d",iPtBin),
 588                                                            Form("Distance between PCA and primary vertex (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",
 589                                                                 ptMin,ptMax), 1000, 0, 50000);
 590       fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA[iPtBin] = new TH2D(Form("fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA_%d",iPtBin),
 591                                                                            Form("DistancePrimaryVtxPCA vs WOffsetMuonPairsAtPCA (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",
 592                                                                                 ptMin,ptMax), 300, 0, 60, 1000, 0, 50000);
 593       if (fEvalDimuonVtxResolution) {
 594         fHistDimuonVtxResolutionX[iPtBin] = new TH1D(Form("fHistDimuonVtxResolutionX_%d",iPtBin),
 595                                                      Form("Dimuon Vtx offset along X (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",ptMin,ptMax),
 596                                                      2000, -1000, 1000);
 597         fHistDimuonVtxResolutionY[iPtBin] = new TH1D(Form("fHistDimuonVtxResolutionY_%d",iPtBin),
 598                                                      Form("Dimuon Vtx offset along Y (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",ptMin,ptMax),
 599                                                      2000, -1000, 1000);
 600         fHistDimuonVtxResolutionZ[iPtBin] = new TH1D(Form("fHistDimuonVtxResolutionZ_%d",iPtBin),
 601                                                      Form("Dimuon Vtx offset along Z (%3.1f < p_{T}^{#mu#mu} < %3.1f GeV/c)",ptMin,ptMax),
 602                                                      2000, -10000, 10000);
 603       }
 604     }
 605
 606     fHistMassMuonPairs[iPtBin]                -> SetXTitle("Mass  [GeV/c^{2}]");
 607     fHistMassMuonPairsWithoutMFT[iPtBin]      -> SetXTitle("Mass  [GeV/c^{2}]");
 608     fHistMassMuonPairsMC[iPtBin]              -> SetXTitle("Mass  [GeV/c^{2}]");
 609     fHistWOffsetMuonPairsAtPrimaryVtx[iPtBin] -> SetXTitle("Weighted Offset at Primary Vtx");
 610     fHistWOffsetMuonPairsAtPCA[iPtBin]        -> SetXTitle("Weighted Offset at PCA");
 611     fHistDistancePrimaryVtxPCA[iPtBin]        -> SetXTitle("PCA - Primary Vtx [#mum]");
 612     fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA[iPtBin] -> SetXTitle("Weighted Offset at PCA");
 613     fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA[iPtBin] -> SetYTitle("PCA - Primary Vtx [#mum]");
 614     if (fEvalDimuonVtxResolution) {
 615       fHistDimuonVtxResolutionX[iPtBin] -> SetXTitle("Offset  [#mum]");
 616       fHistDimuonVtxResolutionY[iPtBin] -> SetXTitle("Offset  [#mum]");
 617       fHistDimuonVtxResolutionZ[iPtBin] -> SetXTitle("Offset  [#mum]");
 618     }
 619
 620     fHistMassMuonPairs[iPtBin]                -> Sumw2();
 621     fHistMassMuonPairsWithoutMFT[iPtBin]      -> Sumw2();
 622     fHistMassMuonPairsMC[iPtBin]              -> Sumw2();
 623     fHistWOffsetMuonPairsAtPrimaryVtx[iPtBin] -> Sumw2();
 624     fHistWOffsetMuonPairsAtPCA[iPtBin]        -> Sumw2();
 625     fHistDistancePrimaryVtxPCA[iPtBin]        -> Sumw2();
 626     fHistDistancePrimaryVtxPCAvsWOffsetMuonPairsAtPCA[iPtBin] -> Sumw2();
 627     if (fEvalDimuonVtxResolution) {
 628       fHistDimuonVtxResolutionX[iPtBin] -> Sumw2();
 629       fHistDimuonVtxResolutionY[iPtBin] -> Sumw2();
 630       fHistDimuonVtxResolutionZ[iPtBin] -> Sumw2();
 631     }
 632
 633   }
 634
 635   fHistMassMuonPairsVsPt           = new TH2D("fHistMassMuonPairsVsPt",           "Dimuon Mass (MUON+MFT) vs p_{T}",  fNMassBins, fMassMin, fMassMax, 20, 0., 10.);
 636   fHistMassMuonPairsWithoutMFTVsPt = new TH2D("fHistMassMuonPairsWithoutMFTVsPt", "Dimuon Mass (MUON only) vs p_{T}", fNMassBins, fMassMin, fMassMax, 20, 0., 10.);
 637         fHistMassMuonPairsVsPtLSp           = new TH2D("fHistMassMuonPairsVsPtLSp",           "Dimuon Mass (MUON+MFT) vs p_{T} Like sign ++",  fNMassBins, fMassMin, fMassMax, 20, 0., 10.);
 638         fHistMassMuonPairsWithoutMFTVsPtLSp = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSp", "Dimuon Mass (MUON only) vs p_{T} Like sign ++", fNMassBins, fMassMin, fMassMax, 20, 0., 10.);
 639         fHistMassMuonPairsVsPtLSm           = new TH2D("fHistMassMuonPairsVsPtLSm",           "Dimuon Mass (MUON+MFT) vs p_{T} Like sign --",  fNMassBins, fMassMin, fMassMax, 20, 0., 10.);
 640         fHistMassMuonPairsWithoutMFTVsPtLSm = new TH2D("fHistMassMuonPairsWithoutMFTVsPtLSm", "Dimuon Mass (MUON only) vs p_{T} Like sign --", fNMassBins, fMassMin, fMassMax, 20, 0., 10.);
 641
 642   fHistMassMuonPairsVsPt           -> SetXTitle("Mass  [GeV/c^{2}]");
 643   fHistMassMuonPairsWithoutMFTVsPt -> SetXTitle("Mass  [GeV/c^{2}]");
 644   fHistMassMuonPairsVsPt           -> SetYTitle("p_{T}  [GeV/c]");
 645   fHistMassMuonPairsWithoutMFTVsPt -> SetYTitle("p_{T}  [GeV/c]");
 646   fHistMassMuonPairsVsPt           -> Sumw2();
 647   fHistMassMuonPairsWithoutMFTVsPt -> Sumw2();
 648         fHistMassMuonPairsVsPtLSp           -> SetXTitle("Mass  [GeV/c^{2}]");
 649         fHistMassMuonPairsWithoutMFTVsPtLSp -> SetXTitle("Mass  [GeV/c^{2}]");
 650         fHistMassMuonPairsVsPtLSp           -> SetYTitle("p_{T}  [GeV/c]");
 651         fHistMassMuonPairsWithoutMFTVsPtLSp -> SetYTitle("p_{T}  [GeV/c]");
 652         fHistMassMuonPairsVsPtLSp           -> Sumw2();
 653         fHistMassMuonPairsWithoutMFTVsPtLSp -> Sumw2();
 654         fHistMassMuonPairsVsPtLSm           -> SetXTitle("Mass  [GeV/c^{2}]");
 655         fHistMassMuonPairsWithoutMFTVsPtLSm -> SetXTitle("Mass  [GeV/c^{2}]");
 656         fHistMassMuonPairsVsPtLSm           -> SetYTitle("p_{T}  [GeV/c]");
 657         fHistMassMuonPairsWithoutMFTVsPtLSm -> SetYTitle("p_{T}  [GeV/c]");
 658         fHistMassMuonPairsVsPtLSm           -> Sumw2();
 659         fHistMassMuonPairsWithoutMFTVsPtLSm -> Sumw2();
 660
 661   if (fOption==kPionsKaons) fHistRapidityPtMuonPairs = new TH2D("fHistRapidityPtMuonPairs", "Dimuon Phase Space (rec)", 20, -4.5, -2., 20, 0., 10.);
 662   else                      fHistRapidityPtMuonPairs = new TH2D("fHistRapidityPtMuonPairs", "Dimuon Phase Space (MC)", 100, -4.5, -2., 100, 0., 10.);
 663   fHistRapidityPtMuonPairs   -> SetXTitle("Rapidity");
 664   fHistRapidityPtMuonPairs   -> SetYTitle("p_{T}  [GeV/c]");
 665   fHistRapidityPtMuonPairs   -> Sumw2();
 666
 667 }
 668
 669 //====================================================================================================================================================
 670