ROOT logo
#if !defined(__CINT__) || defined(__MAKECINT__)
// ROOT includes
#include <TFile.h>
#include <TH1.h>
#include <TCanvas.h>
#include <Riostream.h>
#include <TROOT.h>
#include <TClonesArray.h>
#include <TLorentzVector.h>

// STEER includes
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliCDBManager.h"

// MUON includes
#include "AliMUONCDB.h"
#include "AliMUONTrack.h"
#include "AliMUONVTrackStore.h"
#include "AliMUONTrackParam.h"
#include "AliMUONRecoCheck.h"
#include "AliMUONVCluster.h"
#include "AliMUONRecoParam.h"
#endif

/// \ingroup macros
/// \file DIMUONFakes.C
///
/// \author Ph. Pillot, Subatech, March. 2009
///
/// Macro to study the effects of fake tracks on the dimuon spectra
/// Results are saved in the root file DiFakes.root
/// Results are relevent provided that you use the same recoParams as for the reconstruction

Double_t sigmaCut = -1.;

//-----------------------------------------------------------------------
void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
	         const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
		 const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
{
  
  //Reset ROOT and connect tree file
  gROOT->Reset();
  
  AliLog::SetClassDebugLevel("AliMCEvent",-1);
  
  // File for histograms and histogram booking
  TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
  
  TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.);
  TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
  TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
  TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.);
  TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
  TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
  TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.);
  TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
  TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
  TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0);
  TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
  TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
  TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0);
  TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
  TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
  TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.);
  TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
  TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
  
  // link to reconstructed and simulated tracks
  AliMUONRecoCheck rc(esdFileName, SimDir);
  
  // load necessary data from OCDB
  AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
  AliCDBManager::Instance()->SetSpecificStorage("GRP/GRP/Data","local://.");
  AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
  AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
  if (!recoParam) return;
  
  // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
  sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
  // compute the mask of requested stations from recoParam
  UInt_t requestedStationMask = 0;
  for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
  // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
  Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
  
  TLorentzVector vMu1, vMu2, vDiMu;
  Int_t nReconstructiblePairs = 0;
  Int_t nReconstructedPairs = 0;
  Int_t nMatchedPairs = 0;
  
  // Loop over ESD events
  FirstEvent = TMath::Max(0, FirstEvent);
  LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
  for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
    
    // get reconstructed and simulated tracks
    AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
    AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
    if (!muonTrackStore || !trackRefStore) continue;
    
    // count the number of reconstructible pairs
    Int_t nMuPlus = 0, nMuMinus = 0;
    TIter next(trackRefStore->CreateIterator());
    AliMUONTrack* trackRef;
    while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
      if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
	if (trackRef->GetTrackParamAtVertex()->GetCharge() > 0) nMuPlus++;
	else nMuMinus++;
      }
    }
    nReconstructiblePairs += nMuPlus*nMuMinus;
    
    // loop over ESD tracks and flag them
    const AliESDEvent* esd = rc.GetESDEvent();
    Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
    for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
      
      AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
      
      // skip ghosts
      if (!esdTrack->ContainTrackerData()) continue;
      
      // find the corresponding MUON track
      AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
      
      // try to match the reconstructed track with a simulated one
      Int_t nMatchClusters = 0;
      AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
      
      // take actions according to matching result
      if (matchedTrackRef) {
	
	// flag matched tracks
	esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
	
	// remove already matched trackRefs
	trackRefStore->Remove(*matchedTrackRef);
	
      } else {
	
	// flag fake tracks
	esdTrack->SetLabel(-1);
	
      }
      
    }
    
    // double loop over ESD tracks, build pairs and fill histograms according to their label
    for (Int_t iTrack1 = 0; iTrack1 <  nTracks;  iTrack1++) {
      AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
      
      // skip ghosts
      if (!muonTrack1->ContainTrackerData()) continue;
      
      // get track info
      Short_t charge1 = muonTrack1->Charge();
      Int_t label1 = muonTrack1->GetLabel();
      muonTrack1->LorentzP(vMu1);
      
      for (Int_t iTrack2 = iTrack1+1; iTrack2 <  nTracks;  iTrack2++) {
	AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
	
	// skip ghosts
	if (!muonTrack2->ContainTrackerData()) continue;
	
	// keep only opposite sign pairs
	Short_t charge2 = muonTrack2->Charge();
	if (charge1*charge2 > 0) continue;
	
	nReconstructedPairs++;
	
	// get track info
	Int_t label2 = muonTrack2->GetLabel();
	muonTrack2->LorentzP(vMu2);
	
	// compute kinematics of the pair
	vDiMu = vMu1 + vMu2;
	Float_t mass = vDiMu.M();
	Float_t p = vDiMu.P();
	Float_t pt = vDiMu.Pt();
	Float_t y = vDiMu.Rapidity();
	Float_t eta = vDiMu.Eta();
	Float_t phi = vDiMu.Phi();
	if (phi < 0) phi += 2.*TMath::Pi();
	
	// fill global histograms
	hMass->Fill(mass);
	hP->Fill(p);
	hPt->Fill(pt);
	hY->Fill(y);
	hEta->Fill(eta);
	hPhi->Fill(phi);
	
	// fill histograms according to labels
	if (label1 >= 0 && label2 >= 0) {
	  
	  nMatchedPairs++;
	  
	  hMassM->Fill(mass);
	  hPM->Fill(p);
	  hPtM->Fill(pt);
	  hYM->Fill(y);
	  hEtaM->Fill(eta);
	  hPhiM->Fill(phi);
	  
	} else {
	  
	  hMassF->Fill(mass);
	  hPF->Fill(p);
	  hPtF->Fill(pt);
	  hYF->Fill(y);
	  hEtaF->Fill(eta);
	  hPhiF->Fill(phi);
	  
	}
	
      }
      
    }
      
  } // end of loop over events
  
  // plot results
  TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
  cDiFakesSummary.Divide(3,2);
  cDiFakesSummary.cd(1);
  cDiFakesSummary.GetPad(1)->SetLogy();
  hMass->Draw();
  hMass->SetMinimum(0.5);
  hMassM->Draw("same");
  hMassM->SetLineColor(4);
  hMassF->Draw("same");
  hMassF->SetLineColor(2);
  hMassF->SetFillColor(2);
  hMassF->SetFillStyle(3017);
  cDiFakesSummary.cd(2);
  cDiFakesSummary.GetPad(3)->SetLogy();
  hP->Draw();
  hP->SetMinimum(0.5);
  hPM->Draw("same");
  hPM->SetLineColor(4);
  hPF->Draw("same");
  hPF->SetLineColor(2);
  hPF->SetFillColor(2);
  hPF->SetFillStyle(3017);
  cDiFakesSummary.cd(3);
  cDiFakesSummary.GetPad(4)->SetLogy();
  hPt->Draw();
  hPt->SetMinimum(0.5);
  hPtM->Draw("same");
  hPtM->SetLineColor(4);
  hPtF->Draw("same");
  hPtF->SetLineColor(2);
  hPtF->SetFillColor(2);
  hPtF->SetFillStyle(3017);
  cDiFakesSummary.cd(4);
  cDiFakesSummary.GetPad(2)->SetLogy();
  hY->Draw();
  hY->SetMinimum(0.5);
  hYM->Draw("same");
  hYM->SetLineColor(4);
  hYF->Draw("same");
  hYF->SetLineColor(2);
  hYF->SetFillColor(2);
  hYF->SetFillStyle(3017);
  cDiFakesSummary.cd(5);
  cDiFakesSummary.GetPad(5)->SetLogy();
  hEta->Draw();
  hEta->SetMinimum(0.5);
  hEtaM->Draw("same");
  hEtaM->SetLineColor(4);
  hEtaF->Draw("same");
  hEtaF->SetLineColor(2);
  hEtaF->SetFillColor(2);
  hEtaF->SetFillStyle(3017);
  cDiFakesSummary.cd(6);
  cDiFakesSummary.GetPad(6)->SetLogy();
  hPhi->Draw();
  hPhi->SetMinimum(0.5);
  hPhiM->Draw("same");
  hPhiM->SetLineColor(4);
  hPhiF->Draw("same");
  hPhiF->SetLineColor(2);
  hPhiF->SetFillColor(2);
  hPhiF->SetFillStyle(3017);
  
  // save results
  histoFile->cd();
  histoFile->Write();
  cDiFakesSummary.Write();
  histoFile->Close();
  
  // print results
  printf("\n");
  printf("nb of reconstructible OS pairs: %d \n", nReconstructiblePairs);
  printf("nb of reconstructed OS pairs: %d \n", nReconstructedPairs);
  printf("nb of reconstructed OS pairs matched with trackRefs: %d \n", nMatchedPairs);
  printf("\n");
  
}

 DIMUONFakes.C:1
 DIMUONFakes.C:2
 DIMUONFakes.C:3
 DIMUONFakes.C:4
 DIMUONFakes.C:5
 DIMUONFakes.C:6
 DIMUONFakes.C:7
 DIMUONFakes.C:8
 DIMUONFakes.C:9
 DIMUONFakes.C:10
 DIMUONFakes.C:11
 DIMUONFakes.C:12
 DIMUONFakes.C:13
 DIMUONFakes.C:14
 DIMUONFakes.C:15
 DIMUONFakes.C:16
 DIMUONFakes.C:17
 DIMUONFakes.C:18
 DIMUONFakes.C:19
 DIMUONFakes.C:20
 DIMUONFakes.C:21
 DIMUONFakes.C:22
 DIMUONFakes.C:23
 DIMUONFakes.C:24
 DIMUONFakes.C:25
 DIMUONFakes.C:26
 DIMUONFakes.C:27
 DIMUONFakes.C:28
 DIMUONFakes.C:29
 DIMUONFakes.C:30
 DIMUONFakes.C:31
 DIMUONFakes.C:32
 DIMUONFakes.C:33
 DIMUONFakes.C:34
 DIMUONFakes.C:35
 DIMUONFakes.C:36
 DIMUONFakes.C:37
 DIMUONFakes.C:38
 DIMUONFakes.C:39
 DIMUONFakes.C:40
 DIMUONFakes.C:41
 DIMUONFakes.C:42
 DIMUONFakes.C:43
 DIMUONFakes.C:44
 DIMUONFakes.C:45
 DIMUONFakes.C:46
 DIMUONFakes.C:47
 DIMUONFakes.C:48
 DIMUONFakes.C:49
 DIMUONFakes.C:50
 DIMUONFakes.C:51
 DIMUONFakes.C:52
 DIMUONFakes.C:53
 DIMUONFakes.C:54
 DIMUONFakes.C:55
 DIMUONFakes.C:56
 DIMUONFakes.C:57
 DIMUONFakes.C:58
 DIMUONFakes.C:59
 DIMUONFakes.C:60
 DIMUONFakes.C:61
 DIMUONFakes.C:62
 DIMUONFakes.C:63
 DIMUONFakes.C:64
 DIMUONFakes.C:65
 DIMUONFakes.C:66
 DIMUONFakes.C:67
 DIMUONFakes.C:68
 DIMUONFakes.C:69
 DIMUONFakes.C:70
 DIMUONFakes.C:71
 DIMUONFakes.C:72
 DIMUONFakes.C:73
 DIMUONFakes.C:74
 DIMUONFakes.C:75
 DIMUONFakes.C:76
 DIMUONFakes.C:77
 DIMUONFakes.C:78
 DIMUONFakes.C:79
 DIMUONFakes.C:80
 DIMUONFakes.C:81
 DIMUONFakes.C:82
 DIMUONFakes.C:83
 DIMUONFakes.C:84
 DIMUONFakes.C:85
 DIMUONFakes.C:86
 DIMUONFakes.C:87
 DIMUONFakes.C:88
 DIMUONFakes.C:89
 DIMUONFakes.C:90
 DIMUONFakes.C:91
 DIMUONFakes.C:92
 DIMUONFakes.C:93
 DIMUONFakes.C:94
 DIMUONFakes.C:95
 DIMUONFakes.C:96
 DIMUONFakes.C:97
 DIMUONFakes.C:98
 DIMUONFakes.C:99
 DIMUONFakes.C:100
 DIMUONFakes.C:101
 DIMUONFakes.C:102
 DIMUONFakes.C:103
 DIMUONFakes.C:104
 DIMUONFakes.C:105
 DIMUONFakes.C:106
 DIMUONFakes.C:107
 DIMUONFakes.C:108
 DIMUONFakes.C:109
 DIMUONFakes.C:110
 DIMUONFakes.C:111
 DIMUONFakes.C:112
 DIMUONFakes.C:113
 DIMUONFakes.C:114
 DIMUONFakes.C:115
 DIMUONFakes.C:116
 DIMUONFakes.C:117
 DIMUONFakes.C:118
 DIMUONFakes.C:119
 DIMUONFakes.C:120
 DIMUONFakes.C:121
 DIMUONFakes.C:122
 DIMUONFakes.C:123
 DIMUONFakes.C:124
 DIMUONFakes.C:125
 DIMUONFakes.C:126
 DIMUONFakes.C:127
 DIMUONFakes.C:128
 DIMUONFakes.C:129
 DIMUONFakes.C:130
 DIMUONFakes.C:131
 DIMUONFakes.C:132
 DIMUONFakes.C:133
 DIMUONFakes.C:134
 DIMUONFakes.C:135
 DIMUONFakes.C:136
 DIMUONFakes.C:137
 DIMUONFakes.C:138
 DIMUONFakes.C:139
 DIMUONFakes.C:140
 DIMUONFakes.C:141
 DIMUONFakes.C:142
 DIMUONFakes.C:143
 DIMUONFakes.C:144
 DIMUONFakes.C:145
 DIMUONFakes.C:146
 DIMUONFakes.C:147
 DIMUONFakes.C:148
 DIMUONFakes.C:149
 DIMUONFakes.C:150
 DIMUONFakes.C:151
 DIMUONFakes.C:152
 DIMUONFakes.C:153
 DIMUONFakes.C:154
 DIMUONFakes.C:155
 DIMUONFakes.C:156
 DIMUONFakes.C:157
 DIMUONFakes.C:158
 DIMUONFakes.C:159
 DIMUONFakes.C:160
 DIMUONFakes.C:161
 DIMUONFakes.C:162
 DIMUONFakes.C:163
 DIMUONFakes.C:164
 DIMUONFakes.C:165
 DIMUONFakes.C:166
 DIMUONFakes.C:167
 DIMUONFakes.C:168
 DIMUONFakes.C:169
 DIMUONFakes.C:170
 DIMUONFakes.C:171
 DIMUONFakes.C:172
 DIMUONFakes.C:173
 DIMUONFakes.C:174
 DIMUONFakes.C:175
 DIMUONFakes.C:176
 DIMUONFakes.C:177
 DIMUONFakes.C:178
 DIMUONFakes.C:179
 DIMUONFakes.C:180
 DIMUONFakes.C:181
 DIMUONFakes.C:182
 DIMUONFakes.C:183
 DIMUONFakes.C:184
 DIMUONFakes.C:185
 DIMUONFakes.C:186
 DIMUONFakes.C:187
 DIMUONFakes.C:188
 DIMUONFakes.C:189
 DIMUONFakes.C:190
 DIMUONFakes.C:191
 DIMUONFakes.C:192
 DIMUONFakes.C:193
 DIMUONFakes.C:194
 DIMUONFakes.C:195
 DIMUONFakes.C:196
 DIMUONFakes.C:197
 DIMUONFakes.C:198
 DIMUONFakes.C:199
 DIMUONFakes.C:200
 DIMUONFakes.C:201
 DIMUONFakes.C:202
 DIMUONFakes.C:203
 DIMUONFakes.C:204
 DIMUONFakes.C:205
 DIMUONFakes.C:206
 DIMUONFakes.C:207
 DIMUONFakes.C:208
 DIMUONFakes.C:209
 DIMUONFakes.C:210
 DIMUONFakes.C:211
 DIMUONFakes.C:212
 DIMUONFakes.C:213
 DIMUONFakes.C:214
 DIMUONFakes.C:215
 DIMUONFakes.C:216
 DIMUONFakes.C:217
 DIMUONFakes.C:218
 DIMUONFakes.C:219
 DIMUONFakes.C:220
 DIMUONFakes.C:221
 DIMUONFakes.C:222
 DIMUONFakes.C:223
 DIMUONFakes.C:224
 DIMUONFakes.C:225
 DIMUONFakes.C:226
 DIMUONFakes.C:227
 DIMUONFakes.C:228
 DIMUONFakes.C:229
 DIMUONFakes.C:230
 DIMUONFakes.C:231
 DIMUONFakes.C:232
 DIMUONFakes.C:233
 DIMUONFakes.C:234
 DIMUONFakes.C:235
 DIMUONFakes.C:236
 DIMUONFakes.C:237
 DIMUONFakes.C:238
 DIMUONFakes.C:239
 DIMUONFakes.C:240
 DIMUONFakes.C:241
 DIMUONFakes.C:242
 DIMUONFakes.C:243
 DIMUONFakes.C:244
 DIMUONFakes.C:245
 DIMUONFakes.C:246
 DIMUONFakes.C:247
 DIMUONFakes.C:248
 DIMUONFakes.C:249
 DIMUONFakes.C:250
 DIMUONFakes.C:251
 DIMUONFakes.C:252
 DIMUONFakes.C:253
 DIMUONFakes.C:254
 DIMUONFakes.C:255
 DIMUONFakes.C:256
 DIMUONFakes.C:257
 DIMUONFakes.C:258
 DIMUONFakes.C:259
 DIMUONFakes.C:260
 DIMUONFakes.C:261
 DIMUONFakes.C:262
 DIMUONFakes.C:263
 DIMUONFakes.C:264
 DIMUONFakes.C:265
 DIMUONFakes.C:266
 DIMUONFakes.C:267
 DIMUONFakes.C:268
 DIMUONFakes.C:269
 DIMUONFakes.C:270
 DIMUONFakes.C:271
 DIMUONFakes.C:272
 DIMUONFakes.C:273
 DIMUONFakes.C:274
 DIMUONFakes.C:275
 DIMUONFakes.C:276
 DIMUONFakes.C:277
 DIMUONFakes.C:278
 DIMUONFakes.C:279
 DIMUONFakes.C:280
 DIMUONFakes.C:281
 DIMUONFakes.C:282
 DIMUONFakes.C:283
 DIMUONFakes.C:284
 DIMUONFakes.C:285
 DIMUONFakes.C:286
 DIMUONFakes.C:287
 DIMUONFakes.C:288
 DIMUONFakes.C:289
 DIMUONFakes.C:290
 DIMUONFakes.C:291
 DIMUONFakes.C:292
 DIMUONFakes.C:293
 DIMUONFakes.C:294
 DIMUONFakes.C:295
 DIMUONFakes.C:296
 DIMUONFakes.C:297
 DIMUONFakes.C:298
 DIMUONFakes.C:299
 DIMUONFakes.C:300
 DIMUONFakes.C:301
 DIMUONFakes.C:302
 DIMUONFakes.C:303
 DIMUONFakes.C:304
 DIMUONFakes.C:305