MUON/DIMUONFakes.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 // ROOT includes
   3 #include <TTree.h>
   4 #include <TFile.h>
   5 #include <TH1.h>
   6 #include <TCanvas.h>
   7 #include <Riostream.h>
   8 #include <TROOT.h>
   9 #include <TClonesArray.h>
  10 #include <TGeoGlobalMagField.h>
  11 #include <TLorentzVector.h>
  12
  13 // STEER includes
  14 #include "AliLog.h"
  15 #include "AliMagF.h"
  16 #include "AliESDEvent.h"
  17 #include "AliESDMuonTrack.h"
  18 #include "AliESDMuonCluster.h"
  19 #include "AliCDBPath.h"
  20 #include "AliCDBEntry.h"
  21 #include "AliCDBManager.h"
  22
  23 // MUON includes
  24 #include "AliMUONTrack.h"
  25 #include "AliMUONVTrackStore.h"
  26 #include "AliMUONTrackParam.h"
  27 #include "AliMUONTrackExtrap.h"
  28 #include "AliMUONESDInterface.h"
  29 #include "AliMUONRecoCheck.h"
  30 #include "AliMUONVCluster.h"
  31 #include "AliMUONRecoParam.h"
  32 #endif
  33
  34 /// \ingroup macros
  35 /// \file DIMUONFakes.C
  36 ///
  37 /// \author Ph. Pillot, Subatech, March. 2009
  38 ///
  39 /// Macro to study the effects of fake tracks on the dimuon spectra
  40 /// Results are saved in the root file DiFakes.root
  41 /// Results are relevent provided that you use the same recoParams as for the reconstruction
  42
  43 void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
  44 TTree* GetESDTree(TFile *esdFile);
  45 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
  46 AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
  47                                 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
  48
  49 //-----------------------------------------------------------------------
  50 void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
  51                  const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
  52 {
  53
  54   //Reset ROOT and connect tree file
  55   gROOT->Reset();
  56
  57   // File for histograms and histogram booking
  58   TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
  59
  60   TH1F *hMass = new TH1F("hMass", "Muon mass distribution (GeV/c^{2})", 100, 0., 12.);
  61   TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
  62   TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
  63   TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
  64   TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
  65   TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
  66   TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
  67   TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
  68   TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
  69   TH1F *hY = new TH1F("hY"," Muon rapidity distribution",100,-10,0);
  70   TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
  71   TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
  72   TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
  73   TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
  74   TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
  75   TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
  76   TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
  77   TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
  78
  79   // prepare for analysis
  80   AliMUONRecoParam *recoParam = 0x0;
  81   Double_t sigmaCut = -1;
  82   Prepare(recoParam, sigmaCut);
  83
  84   // link to reconstructed tracks
  85   TFile* esdFile = TFile::Open(esdFileName);
  86   TTree* esdTree = GetESDTree(esdFile);
  87   AliESDEvent* esd = new AliESDEvent();
  88   esd->ReadFromTree(esdTree);
  89
  90   // link to simulated tracks
  91   AliMUONRecoCheck rc(esdFileName, SimDir);
  92
  93   TLorentzVector vMu1, vMu2, vDiMu;
  94
  95   // Loop over ESD events
  96   FirstEvent = TMath::Max(0, FirstEvent);
  97   LastEvent = (LastEvent>=0) ? TMath::Min((Int_t)esdTree->GetEntries() - 1, LastEvent) : (Int_t)esdTree->GetEntries() - 1;
  98   for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
  99
 100     // get the ESD of current event
 101     esdTree->GetEvent(iEvent);
 102     if (!esd) {
 103       Error("CheckESD", "no ESD object found for event %d", iEvent);
 104       return;
 105     }
 106
 107     // convert TrackRef to MUON tracks
 108     AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
 109
 110     // loop over ESD tracks and flag them
 111     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
 112     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 113
 114       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 115
 116       // skip ghosts
 117       if (!muonTrack->ContainTrackerData()) continue;
 118
 119       // try to match the reconstructed track with a simulated one
 120       Float_t fractionOfMatchCluster = 0.;
 121       AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
 122
 123       // take actions according to matching result
 124       if (matchedTrackRef) {
 125
 126         // flag matched tracks
 127         muonTrack->SetLabel(matchedTrackRef->GetUniqueID());
 128
 129         // remove already matched trackRefs
 130         trackRefStore->Remove(*matchedTrackRef);
 131
 132       } else {
 133
 134         // flag fake tracks
 135         muonTrack->SetLabel(-1);
 136
 137       }
 138
 139     }
 140
 141     // double loop over ESD tracks, build pairs and fill histograms according to their label
 142     for (Int_t iTrack1 = 0; iTrack1 <  nTracks;  iTrack1++) {
 143       AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
 144
 145       // skip ghosts
 146       if (!muonTrack1->ContainTrackerData()) continue;
 147
 148       // get track info
 149       Short_t charge1 = muonTrack1->Charge();
 150       Int_t label1 = muonTrack1->GetLabel();
 151       muonTrack1->LorentzP(vMu1);
 152
 153       for (Int_t iTrack2 = iTrack1+1; iTrack2 <  nTracks;  iTrack2++) {
 154         AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
 155
 156         // skip ghosts
 157         if (!muonTrack2->ContainTrackerData()) continue;
 158
 159         // keep only opposite sign pairs
 160         Short_t charge2 = muonTrack2->Charge();
 161         if (charge1*charge2 > 0) continue;
 162
 163         // get track info
 164         Int_t label2 = muonTrack2->GetLabel();
 165         muonTrack2->LorentzP(vMu2);
 166
 167         // compute kinematics of the pair
 168         vDiMu = vMu1 + vMu2;
 169         Float_t mass = vDiMu.M();
 170         Float_t p = vDiMu.P();
 171         Float_t pt = vDiMu.Pt();
 172         Float_t y = vDiMu.Rapidity();
 173         Float_t eta = vDiMu.Eta();
 174         Float_t phi = vMu1.Phi();
 175         if (phi < 0) phi += 2.*TMath::Pi();
 176
 177         // fill global histograms
 178         hMass->Fill(mass);
 179         hP->Fill(p);
 180         hPt->Fill(pt);
 181         hY->Fill(y);
 182         hEta->Fill(eta);
 183         hPhi->Fill(phi);
 184
 185         // fill histograms according to labels
 186         if (label1 >= 0 && label2 >= 0) {
 187
 188           hMassM->Fill(mass);
 189           hPM->Fill(p);
 190           hPtM->Fill(pt);
 191           hYM->Fill(y);
 192           hEtaM->Fill(eta);
 193           hPhiM->Fill(phi);
 194
 195         } else {
 196
 197           hMassF->Fill(mass);
 198           hPF->Fill(p);
 199           hPtF->Fill(pt);
 200           hYF->Fill(y);
 201           hEtaF->Fill(eta);
 202           hPhiF->Fill(phi);
 203
 204         }
 205
 206       }
 207
 208     }
 209
 210   } // end of loop over events
 211
 212   // plot results
 213   TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
 214   cDiFakesSummary.Divide(3,2);
 215   cDiFakesSummary.cd(1);
 216   cDiFakesSummary.GetPad(1)->SetLogy();
 217   hMass->Draw();
 218   hMass->SetMinimum(0.5);
 219   hMassM->Draw("same");
 220   hMassM->SetLineColor(4);
 221   hMassF->Draw("same");
 222   hMassF->SetLineColor(2);
 223   hMassF->SetFillColor(2);
 224   hMassF->SetFillStyle(3017);
 225   cDiFakesSummary.cd(2);
 226   cDiFakesSummary.GetPad(3)->SetLogy();
 227   hP->Draw();
 228   hP->SetMinimum(0.5);
 229   hPM->Draw("same");
 230   hPM->SetLineColor(4);
 231   hPF->Draw("same");
 232   hPF->SetLineColor(2);
 233   hPF->SetFillColor(2);
 234   hPF->SetFillStyle(3017);
 235   cDiFakesSummary.cd(3);
 236   cDiFakesSummary.GetPad(4)->SetLogy();
 237   hPt->Draw();
 238   hPt->SetMinimum(0.5);
 239   hPtM->Draw("same");
 240   hPtM->SetLineColor(4);
 241   hPtF->Draw("same");
 242   hPtF->SetLineColor(2);
 243   hPtF->SetFillColor(2);
 244   hPtF->SetFillStyle(3017);
 245   cDiFakesSummary.cd(4);
 246   cDiFakesSummary.GetPad(2)->SetLogy();
 247   hY->Draw();
 248   hY->SetMinimum(0.5);
 249   hYM->Draw("same");
 250   hYM->SetLineColor(4);
 251   hYF->Draw("same");
 252   hYF->SetLineColor(2);
 253   hYF->SetFillColor(2);
 254   hYF->SetFillStyle(3017);
 255   cDiFakesSummary.cd(5);
 256   cDiFakesSummary.GetPad(5)->SetLogy();
 257   hEta->Draw();
 258   hEta->SetMinimum(0.5);
 259   hEtaM->Draw("same");
 260   hEtaM->SetLineColor(4);
 261   hEtaF->Draw("same");
 262   hEtaF->SetLineColor(2);
 263   hEtaF->SetFillColor(2);
 264   hEtaF->SetFillStyle(3017);
 265   cDiFakesSummary.cd(6);
 266   cDiFakesSummary.GetPad(6)->SetLogy();
 267   hPhi->Draw();
 268   hPhi->SetMinimum(0.5);
 269   hPhiM->Draw("same");
 270   hPhiM->SetLineColor(4);
 271   hPhiF->Draw("same");
 272   hPhiF->SetLineColor(2);
 273   hPhiF->SetFillColor(2);
 274   hPhiF->SetFillStyle(3017);
 275
 276   // save results
 277   histoFile->cd();
 278   histoFile->Write();
 279   cDiFakesSummary.Write();
 280   histoFile->Close();
 281
 282   // clear memory
 283   delete esd;
 284   esdFile->Close();
 285   delete recoParam;
 286
 287 }
 288
 289 //-----------------------------------------------------------------------
 290 void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut)
 291 {
 292   /// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef
 293
 294   // prepare OCDB access
 295   AliCDBManager* man = AliCDBManager::Instance();
 296   man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 297   man->SetRun(0);
 298
 299   // set  mag field
 300   // waiting for mag field in CDB
 301   if (!TGeoGlobalMagField::Instance()->GetField()) {
 302     printf("Loading field map...\n");
 303     AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
 304     TGeoGlobalMagField::Instance()->SetField(field);
 305   }
 306   // set the magnetic field for track extrapolations
 307   AliMUONTrackExtrap::SetField();
 308
 309   // Load initial reconstruction parameters from OCDB
 310   AliCDBPath path("MUON","Calib","RecoParam");
 311   AliCDBEntry *entry=man->Get(path.GetPath());
 312   if(entry) {
 313     recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
 314     entry->SetOwner(0);
 315     AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
 316   }
 317   if (!recoParam) {
 318     printf("Couldn't find RecoParam object in OCDB: create default one");
 319     recoParam = AliMUONRecoParam::GetLowFluxParam();
 320   }
 321
 322   Info("MUONFakes", "\n recontruction parameters:");
 323   recoParam->Print("FULL");
 324   AliMUONESDInterface::ResetTracker(recoParam);
 325
 326   // sigma cut to associate clusters with TrackRefs in case the label are not used
 327   sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
 328
 329 }
 330
 331 //-----------------------------------------------------------------------
 332 TTree* GetESDTree(TFile *esdFile)
 333 {
 334   /// Check that the file is properly open
 335   /// Return pointer to the ESD Tree
 336
 337   if (!esdFile || !esdFile->IsOpen()) {
 338     Error("GetESDTree", "opening ESD file failed");
 339     exit(-1);
 340   }
 341
 342   TTree* tree = (TTree*) esdFile->Get("esdTree");
 343   if (!tree) {
 344     Error("GetESDTree", "no ESD tree found");
 345     exit(-1);
 346   }
 347
 348   return tree;
 349
 350 }
 351
 352 //-----------------------------------------------------------------------
 353 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
 354 {
 355   /// Try to match 2 tracks
 356
 357   Bool_t compTrack[10];
 358   Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
 359   fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
 360
 361   if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
 362       (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
 363       fractionOfMatchCluster > 0.5) return kTRUE;                       // more than 50% of clusters matched
 364   else return kFALSE;
 365
 366 }
 367
 368 //-----------------------------------------------------------------------
 369 AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
 370                                 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
 371 {
 372   /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
 373
 374   AliMUONTrack *matchedTrackRef = 0x0;
 375   fractionOfMatchCluster = 0.;
 376
 377   if (useLabel) { // by using the MC label
 378
 379     // get the corresponding simulated track if any
 380     Int_t label = muonTrack.GetLabel();
 381     matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
 382
 383     // get the fraction of matched clusters
 384     if (matchedTrackRef) {
 385       Int_t nMatchClusters = 0;
 386       if (muonTrack.ClustersStored()) {
 387         AliESDMuonCluster* cluster = (AliESDMuonCluster*) muonTrack.GetClusters().First();
 388         while (cluster) {
 389           if (cluster->GetLabel() == label) nMatchClusters++;
 390           cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster);
 391         }
 392       }
 393       fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters());
 394     }
 395
 396   } else { // by comparing cluster/TrackRef positions
 397
 398     // convert ESD track to MUON track
 399     AliMUONTrack track;
 400     AliMUONESDInterface::ESDToMUON(muonTrack,track);
 401
 402     // look for the corresponding simulated track if any
 403     TIter next(trackRefStore.CreateIterator());
 404     AliMUONTrack* trackRef;
 405     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
 406
 407       // check compatibility
 408       Float_t f = 0.;
 409       if (TrackMatched(track, *trackRef, f, sigmaCut)) {
 410         matchedTrackRef = trackRef;
 411         fractionOfMatchCluster = f;
 412         break;
 413       }
 414
 415     }
 416
 417   }
 418
 419   return matchedTrackRef;
 420
 421 }
 422