MUON/DIMUONFakes.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 // ROOT includes
   3 #include <TFile.h>
   4 #include <TH1.h>
   5 #include <TCanvas.h>
   6 #include <Riostream.h>
   7 #include <TROOT.h>
   8 #include <TClonesArray.h>
   9 #include <TLorentzVector.h>
  10
  11 // STEER includes
  12 #include "AliLog.h"
  13 #include "AliESDEvent.h"
  14 #include "AliESDMuonTrack.h"
  15 #include "AliCDBManager.h"
  16
  17 // MUON includes
  18 #include "AliMUONCDB.h"
  19 #include "AliMUONTrack.h"
  20 #include "AliMUONVTrackStore.h"
  21 #include "AliMUONTrackParam.h"
  22 #include "AliMUONRecoCheck.h"
  23 #include "AliMUONVCluster.h"
  24 #include "AliMUONRecoParam.h"
  25 #endif
  26
  27 /// \ingroup macros
  28 /// \file DIMUONFakes.C
  29 ///
  30 /// \author Ph. Pillot, Subatech, March. 2009
  31 ///
  32 /// Macro to study the effects of fake tracks on the dimuon spectra
  33 /// Results are saved in the root file DiFakes.root
  34 /// Results are relevent provided that you use the same recoParams as for the reconstruction
  35
  36 Double_t sigmaCut = -1.;
  37
  38 //-----------------------------------------------------------------------
  39 void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
  40                  const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
  41                  const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
  42 {
  43
  44   //Reset ROOT and connect tree file
  45   gROOT->Reset();
  46
  47   AliLog::SetClassDebugLevel("AliMCEvent",-1);
  48
  49   // File for histograms and histogram booking
  50   TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
  51
  52   TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.);
  53   TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
  54   TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
  55   TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.);
  56   TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
  57   TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
  58   TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.);
  59   TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
  60   TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
  61   TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0);
  62   TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
  63   TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
  64   TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0);
  65   TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
  66   TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
  67   TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.);
  68   TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
  69   TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
  70
  71   // link to reconstructed and simulated tracks
  72   AliMUONRecoCheck rc(esdFileName, SimDir);
  73
  74   // load necessary data from OCDB
  75   AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
  76   AliCDBManager::Instance()->SetSpecificStorage("GRP/GRP/Data","local://.");
  77   AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
  78   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
  79   if (!recoParam) return;
  80
  81   // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
  82   sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
  83   // compute the mask of requested stations from recoParam
  84   UInt_t requestedStationMask = 0;
  85   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
  86   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
  87   Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
  88
  89   TLorentzVector vMu1, vMu2, vDiMu;
  90   Int_t nReconstructiblePairs = 0;
  91   Int_t nReconstructedPairs = 0;
  92   Int_t nMatchedPairs = 0;
  93
  94   // Loop over ESD events
  95   FirstEvent = TMath::Max(0, FirstEvent);
  96   LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
  97   for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
  98
  99     // get reconstructed and simulated tracks
 100     AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
 101     AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
 102     if (!muonTrackStore || !trackRefStore) continue;
 103
 104     // count the number of reconstructible pairs
 105     Int_t nMuPlus = 0, nMuMinus = 0;
 106     TIter next(trackRefStore->CreateIterator());
 107     AliMUONTrack* trackRef;
 108     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
 109       if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
 110         if (trackRef->GetTrackParamAtVertex()->GetCharge() > 0) nMuPlus++;
 111         else nMuMinus++;
 112       }
 113     }
 114     nReconstructiblePairs += nMuPlus*nMuMinus;
 115
 116     // loop over ESD tracks and flag them
 117     const AliESDEvent* esd = rc.GetESDEvent();
 118     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
 119     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 120
 121       AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
 122
 123       // skip ghosts
 124       if (!esdTrack->ContainTrackerData()) continue;
 125
 126       // find the corresponding MUON track
 127       AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
 128
 129       // try to match the reconstructed track with a simulated one
 130       Int_t nMatchClusters = 0;
 131       AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
 132
 133       // take actions according to matching result
 134       if (matchedTrackRef) {
 135
 136         // flag matched tracks
 137         esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
 138
 139         // remove already matched trackRefs
 140         trackRefStore->Remove(*matchedTrackRef);
 141
 142       } else {
 143
 144         // flag fake tracks
 145         esdTrack->SetLabel(-1);
 146
 147       }
 148
 149     }
 150
 151     // double loop over ESD tracks, build pairs and fill histograms according to their label
 152     for (Int_t iTrack1 = 0; iTrack1 <  nTracks;  iTrack1++) {
 153       AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
 154
 155       // skip ghosts
 156       if (!muonTrack1->ContainTrackerData()) continue;
 157
 158       // get track info
 159       Short_t charge1 = muonTrack1->Charge();
 160       Int_t label1 = muonTrack1->GetLabel();
 161       muonTrack1->LorentzP(vMu1);
 162
 163       for (Int_t iTrack2 = iTrack1+1; iTrack2 <  nTracks;  iTrack2++) {
 164         AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
 165
 166         // skip ghosts
 167         if (!muonTrack2->ContainTrackerData()) continue;
 168
 169         // keep only opposite sign pairs
 170         Short_t charge2 = muonTrack2->Charge();
 171         if (charge1*charge2 > 0) continue;
 172
 173         nReconstructedPairs++;
 174
 175         // get track info
 176         Int_t label2 = muonTrack2->GetLabel();
 177         muonTrack2->LorentzP(vMu2);
 178
 179         // compute kinematics of the pair
 180         vDiMu = vMu1 + vMu2;
 181         Float_t mass = vDiMu.M();
 182         Float_t p = vDiMu.P();
 183         Float_t pt = vDiMu.Pt();
 184         Float_t y = vDiMu.Rapidity();
 185         Float_t eta = vDiMu.Eta();
 186         Float_t phi = vDiMu.Phi();
 187         if (phi < 0) phi += 2.*TMath::Pi();
 188
 189         // fill global histograms
 190         hMass->Fill(mass);
 191         hP->Fill(p);
 192         hPt->Fill(pt);
 193         hY->Fill(y);
 194         hEta->Fill(eta);
 195         hPhi->Fill(phi);
 196
 197         // fill histograms according to labels
 198         if (label1 >= 0 && label2 >= 0) {
 199
 200           nMatchedPairs++;
 201
 202           hMassM->Fill(mass);
 203           hPM->Fill(p);
 204           hPtM->Fill(pt);
 205           hYM->Fill(y);
 206           hEtaM->Fill(eta);
 207           hPhiM->Fill(phi);
 208
 209         } else {
 210
 211           hMassF->Fill(mass);
 212           hPF->Fill(p);
 213           hPtF->Fill(pt);
 214           hYF->Fill(y);
 215           hEtaF->Fill(eta);
 216           hPhiF->Fill(phi);
 217
 218         }
 219
 220       }
 221
 222     }
 223
 224   } // end of loop over events
 225
 226   // plot results
 227   TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
 228   cDiFakesSummary.Divide(3,2);
 229   cDiFakesSummary.cd(1);
 230   cDiFakesSummary.GetPad(1)->SetLogy();
 231   hMass->Draw();
 232   hMass->SetMinimum(0.5);
 233   hMassM->Draw("same");
 234   hMassM->SetLineColor(4);
 235   hMassF->Draw("same");
 236   hMassF->SetLineColor(2);
 237   hMassF->SetFillColor(2);
 238   hMassF->SetFillStyle(3017);
 239   cDiFakesSummary.cd(2);
 240   cDiFakesSummary.GetPad(3)->SetLogy();
 241   hP->Draw();
 242   hP->SetMinimum(0.5);
 243   hPM->Draw("same");
 244   hPM->SetLineColor(4);
 245   hPF->Draw("same");
 246   hPF->SetLineColor(2);
 247   hPF->SetFillColor(2);
 248   hPF->SetFillStyle(3017);
 249   cDiFakesSummary.cd(3);
 250   cDiFakesSummary.GetPad(4)->SetLogy();
 251   hPt->Draw();
 252   hPt->SetMinimum(0.5);
 253   hPtM->Draw("same");
 254   hPtM->SetLineColor(4);
 255   hPtF->Draw("same");
 256   hPtF->SetLineColor(2);
 257   hPtF->SetFillColor(2);
 258   hPtF->SetFillStyle(3017);
 259   cDiFakesSummary.cd(4);
 260   cDiFakesSummary.GetPad(2)->SetLogy();
 261   hY->Draw();
 262   hY->SetMinimum(0.5);
 263   hYM->Draw("same");
 264   hYM->SetLineColor(4);
 265   hYF->Draw("same");
 266   hYF->SetLineColor(2);
 267   hYF->SetFillColor(2);
 268   hYF->SetFillStyle(3017);
 269   cDiFakesSummary.cd(5);
 270   cDiFakesSummary.GetPad(5)->SetLogy();
 271   hEta->Draw();
 272   hEta->SetMinimum(0.5);
 273   hEtaM->Draw("same");
 274   hEtaM->SetLineColor(4);
 275   hEtaF->Draw("same");
 276   hEtaF->SetLineColor(2);
 277   hEtaF->SetFillColor(2);
 278   hEtaF->SetFillStyle(3017);
 279   cDiFakesSummary.cd(6);
 280   cDiFakesSummary.GetPad(6)->SetLogy();
 281   hPhi->Draw();
 282   hPhi->SetMinimum(0.5);
 283   hPhiM->Draw("same");
 284   hPhiM->SetLineColor(4);
 285   hPhiF->Draw("same");
 286   hPhiF->SetLineColor(2);
 287   hPhiF->SetFillColor(2);
 288   hPhiF->SetFillStyle(3017);
 289
 290   // save results
 291   histoFile->cd();
 292   histoFile->Write();
 293   cDiFakesSummary.Write();
 294   histoFile->Close();
 295
 296   // print results
 297   printf("\n");
 298   printf("nb of reconstructible OS pairs: %d \n", nReconstructiblePairs);
 299   printf("nb of reconstructed OS pairs: %d \n", nReconstructedPairs);
 300   printf("nb of reconstructed OS pairs matched with trackRefs: %d \n", nMatchedPairs);
 301   printf("\n");
 302
 303 }
 304