MUON/macros/MUONFakes.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 <TObjArray.h>
   9 #include <TArrayI.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 "AliMUONESDInterface.h"
  23 #include "AliMUONRecoCheck.h"
  24 #include "AliMUONVCluster.h"
  25 #include "AliMUONRecoParam.h"
  26 #endif
  27
  28 /// \ingroup macros
  29 /// \file MUONFakes.C
  30 ///
  31 /// \author Ph. Pillot, Subatech, March. 2009
  32 ///
  33 /// Macro to study fake tracks by comparing reconstructed tracks with TrackRefs
  34 /// Results are saved in the root file Fakes.root
  35 /// Results are relevent provided that you use the same recoParams as for the reconstruction
  36
  37 UInt_t requestedStationMask = 0;
  38 Bool_t request2ChInSameSt45 = kFALSE;
  39 Double_t sigmaCut = -1.;
  40
  41 Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
  42                            Bool_t useLabel, TH1F &hFractionOfConnectedClusters);
  43
  44 //-----------------------------------------------------------------------
  45 void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
  46                const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
  47                const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
  48 {
  49
  50   // Disable printout of AliMCEvent
  51   AliLog::SetClassDebugLevel("AliMCEvent",-1);
  52
  53   //Reset ROOT and connect tree file
  54   gROOT->Reset();
  55
  56   // File for histograms and histogram booking
  57   TFile *histoFile = new TFile("Fakes.root", "RECREATE");
  58
  59   TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks","nb of tracks /evt",20,0.,20.);
  60   TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks","nb of fake - nb of missing track",20,0.,20.);
  61
  62   TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters","nb of clusters /track",20,0.,20.);
  63   TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM","nb of clusters /matched track",20,0.,20.);
  64   TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF","nb of clusters /fake track",20,0.,20.);
  65   TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC","nb of clusters /MC track",20,0.,20.);
  66   TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters","nb of matched clusters / nb of clusters",110,0.,1.1);
  67   TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters","nb of connected clusters / nb of clusters in fake tracks",110,0.,1.1);
  68
  69   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 100, 0., 20.);
  70   TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 100, 0., 20.);
  71   TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 100, 0., 20.);
  72   TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
  73   TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
  74   TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
  75   TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
  76   TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
  77   TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
  78   TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
  79   TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
  80   TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
  81   TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
  82   TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
  83   TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
  84   TH1F *hNumberOfChamberHit = new TH1F("hNumberOfChamberHit","nb of chambers hit /track",20,0.,20.);
  85   TH1F *hNumberOfChamberHitM = new TH1F("hNumberOfChamberHitM","nb of chambers hit /matched track",20,0.,20.);
  86   TH1F *hNumberOfChamberHitF = new TH1F("hNumberOfChamberHitF","nb of chambers hit /fake track",20,0.,20.);
  87
  88   // link to reconstructed and simulated tracks
  89   AliMUONRecoCheck rc(esdFileName, SimDir);
  90
  91   // load necessary data from OCDB
  92   AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
  93   AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
  94   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
  95   if (!recoParam) return;
  96
  97   // get sigma cut from recoParam to associate clusters with TrackRefs in case the labels are not used
  98   sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
  99   // compute the mask of requested stations from recoParam
 100   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
 101   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
 102   request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
 103
 104   // initialize global counters
 105   Int_t nReconstructibleTracks = 0;
 106   Int_t nReconstructedTracks = 0;
 107   Int_t nEventsWithTrackReconstructedYet = 0;
 108   Int_t nEventsWithFake = 0;
 109   Int_t nEventsWithAdditionalFake = 0;
 110   Int_t nTotMatchedTracks = 0;
 111   Int_t nTotTracksReconstructedYet = 0;
 112   Int_t nTotFakeTracks = 0;
 113   Int_t nTotConnectedTracks = 0;
 114   Int_t nTotAdditionalTracks = 0;
 115   Bool_t trackReconstructedYet;
 116   TArrayI eventsWithTrackReconstructedYet(10);
 117   TArrayI eventsWithFake(10);
 118   TArrayI eventsWithAdditionalFake(10);
 119
 120   // Loop over ESD events
 121   FirstEvent = TMath::Max(0, FirstEvent);
 122   LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
 123   for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
 124
 125     // get reconstructed and simulated tracks
 126     AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
 127     AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
 128     if (!muonTrackStore || !trackRefStore) continue;
 129
 130     // count the number of reconstructible tracks
 131     TIter next(trackRefStore->CreateIterator());
 132     AliMUONTrack* trackRef;
 133     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
 134       if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) nReconstructibleTracks++;
 135     }
 136
 137     // loop over ESD tracks
 138     Int_t nTrackerTracks = 0;
 139     trackReconstructedYet = kFALSE;
 140     AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
 141     const AliESDEvent* esd = rc.GetESDEvent();
 142     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
 143     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 144
 145       AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
 146
 147       // skip ghosts
 148       if (!esdTrack->ContainTrackerData()) continue;
 149       nTrackerTracks++;
 150
 151       // find the corresponding MUON track
 152       AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
 153
 154       // get track info
 155       Int_t nClusters = esdTrack->GetNClusters();
 156       Double_t normalizedChi2 = esdTrack->GetChi2() / (2. * esdTrack->GetNHit() - 5);
 157       Double_t p = esdTrack->P();
 158       Double_t pT = esdTrack->Pt();
 159       Double_t eta = esdTrack->Eta();
 160       Double_t phi = esdTrack->Phi();
 161       Int_t nChamberHit = 0;
 162       for (Int_t ich=0; ich<10; ich++) if (esdTrack->IsInMuonClusterMap(ich)) nChamberHit++;
 163
 164       // fill global histograms
 165       hNumberOfClusters->Fill(nClusters);
 166       hChi2PerDof->Fill(normalizedChi2);
 167       hP->Fill(p);
 168       hPt->Fill(pT);
 169       hEta->Fill(eta);
 170       hPhi->Fill(phi);
 171       hNumberOfChamberHit->Fill(nChamberHit);
 172
 173       // try to match the reconstructed track with a simulated one
 174       Int_t nMatchClusters = 0;
 175       AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
 176
 177       // take actions according to matching result
 178       if (matchedTrackRef) {
 179
 180         // global counter
 181         nTotMatchedTracks++;
 182         if (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
 183           trackReconstructedYet = kTRUE;
 184           nTotTracksReconstructedYet++;
 185         }
 186
 187         // fill histograms
 188         hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
 189         hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
 190         hNumberOfClustersM->Fill(nClusters);
 191         hChi2PerDofM->Fill(normalizedChi2);
 192         hPM->Fill(p);
 193         hPtM->Fill(pT);
 194         hEtaM->Fill(eta);
 195         hPhiM->Fill(phi);
 196         hNumberOfChamberHitM->Fill(nChamberHit);
 197
 198         // remove already matched trackRefs
 199         trackRefStore->Remove(*matchedTrackRef);
 200
 201       } else {
 202
 203         // global counter
 204         nTotFakeTracks++;
 205
 206         // fill histograms
 207         hNumberOfClustersF->Fill(nClusters);
 208         hChi2PerDofF->Fill(normalizedChi2);
 209         hPF->Fill(p);
 210         hPtF->Fill(pT);
 211         hEtaF->Fill(eta);
 212         hPhiF->Fill(phi);
 213         hNumberOfChamberHitF->Fill(nChamberHit);
 214
 215         // store fake tracks
 216         fakeTrackStore->Add(*muonTrack);
 217
 218       }
 219
 220     } // end of loop over ESD tracks
 221
 222     // fill histograms
 223     hNumberOfTracks->Fill(nTrackerTracks);
 224     nReconstructedTracks += nTrackerTracks;
 225     if (trackReconstructedYet) {
 226       if (nEventsWithTrackReconstructedYet >= eventsWithTrackReconstructedYet.GetSize())
 227         eventsWithTrackReconstructedYet.Set(10*eventsWithTrackReconstructedYet.GetSize());
 228       eventsWithTrackReconstructedYet[nEventsWithTrackReconstructedYet++] = iEvent;
 229     }
 230
 231     // count the number the additional fake tracks
 232     if (fakeTrackStore->GetSize() > 0) {
 233
 234       // remove the most connected fake tracks
 235       Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore,  useLabel, *hFractionOfConnectedClusters);
 236
 237       // remove the remaining free reconstructible tracks
 238       Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
 239
 240       // fill histograms
 241       if (nEventsWithFake >= eventsWithFake.GetSize()) eventsWithFake.Set(10*eventsWithFake.GetSize());
 242       eventsWithFake[nEventsWithFake++] = iEvent;
 243       if (nAdditionalTracks > 0) {
 244         if (nEventsWithAdditionalFake >= eventsWithAdditionalFake.GetSize()) eventsWithAdditionalFake.Set(10*eventsWithAdditionalFake.GetSize());
 245         eventsWithAdditionalFake[nEventsWithAdditionalFake++] = iEvent;
 246         nTotAdditionalTracks += nAdditionalTracks;
 247         hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
 248       }
 249
 250     }
 251
 252     delete fakeTrackStore;
 253
 254   } // end of loop over events
 255
 256   // total number of connected tracks
 257   nTotConnectedTracks = hFractionOfConnectedClusters->GetEntries();
 258
 259   // plot results
 260   TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
 261   cFakesSummary.Divide(3,2);
 262   cFakesSummary.cd(1);
 263   TVirtualPad* pad1 = cFakesSummary.GetPad(1);
 264   pad1->Divide(0,2);
 265   pad1->cd(1);
 266   pad1->GetPad(1)->SetLogy();
 267   hNumberOfClusters->Draw();
 268   hNumberOfClusters->SetMinimum(0.5);
 269   hNumberOfClustersM->Draw("same");
 270   hNumberOfClustersM->SetLineColor(4);
 271   hNumberOfClustersF->Draw("same");
 272   hNumberOfClustersF->SetLineColor(2);
 273   hNumberOfClustersF->SetFillColor(2);
 274   hNumberOfClustersF->SetFillStyle(3017);
 275   pad1->cd(2);
 276   pad1->GetPad(2)->SetLogy();
 277   hNumberOfChamberHit->Draw();
 278   hNumberOfChamberHit->SetMinimum(0.5);
 279   hNumberOfChamberHitM->Draw("same");
 280   hNumberOfChamberHitM->SetLineColor(4);
 281   hNumberOfChamberHitF->Draw("same");
 282   hNumberOfChamberHitF->SetLineColor(2);
 283   hNumberOfChamberHitF->SetFillColor(2);
 284   hNumberOfChamberHitF->SetFillStyle(3017);
 285   cFakesSummary.cd(2);
 286   cFakesSummary.GetPad(2)->SetLogy();
 287   hChi2PerDof->Draw();
 288   hChi2PerDof->SetMinimum(0.5);
 289   hChi2PerDofM->Draw("same");
 290   hChi2PerDofM->SetLineColor(4);
 291   hChi2PerDofF->Draw("same");
 292   hChi2PerDofF->SetLineColor(2);
 293   hChi2PerDofF->SetFillColor(2);
 294   hChi2PerDofF->SetFillStyle(3017);
 295   cFakesSummary.cd(3);
 296   cFakesSummary.GetPad(3)->SetLogy();
 297   hP->Draw();
 298   hP->SetMinimum(0.5);
 299   hPM->Draw("same");
 300   hPM->SetLineColor(4);
 301   hPF->Draw("same");
 302   hPF->SetLineColor(2);
 303   hPF->SetFillColor(2);
 304   hPF->SetFillStyle(3017);
 305   cFakesSummary.cd(4);
 306   cFakesSummary.GetPad(4)->SetLogy();
 307   hPt->Draw();
 308   hPt->SetMinimum(0.5);
 309   hPtM->Draw("same");
 310   hPtM->SetLineColor(4);
 311   hPtF->Draw("same");
 312   hPtF->SetLineColor(2);
 313   hPtF->SetFillColor(2);
 314   hPtF->SetFillStyle(3017);
 315   cFakesSummary.cd(5);
 316   cFakesSummary.GetPad(5)->SetLogy();
 317   hEta->Draw();
 318   hEta->SetMinimum(0.5);
 319   hEtaM->Draw("same");
 320   hEtaM->SetLineColor(4);
 321   hEtaF->Draw("same");
 322   hEtaF->SetLineColor(2);
 323   hEtaF->SetFillColor(2);
 324   hEtaF->SetFillStyle(3017);
 325   cFakesSummary.cd(6);
 326   cFakesSummary.GetPad(6)->SetLogy();
 327   hPhi->Draw();
 328   hPhi->SetMinimum(0.5);
 329   hPhiM->Draw("same");
 330   hPhiM->SetLineColor(4);
 331   hPhiF->Draw("same");
 332   hPhiF->SetLineColor(2);
 333   hPhiF->SetFillColor(2);
 334   hPhiF->SetFillStyle(3017);
 335
 336   // save results
 337   histoFile->cd();
 338   histoFile->Write();
 339   cFakesSummary.Write();
 340   histoFile->Close();
 341
 342   // print results
 343   cout << endl;
 344   cout << "- Number of reconstructible tracks: " << nReconstructibleTracks << endl;
 345   cout << "- Number of reconstructed tracks: " << nReconstructedTracks << endl;
 346   cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
 347   cout << "  (including " << nTotTracksReconstructedYet << " track(s) matched with a TrackRef that is not reconstructible";
 348   if (nTotTracksReconstructedYet > 0) {
 349     for(Int_t i=0; i<nEventsWithTrackReconstructedYet; i++){
 350       if (i==0) cout << " (eventID = " << eventsWithTrackReconstructedYet[i];
 351       else cout << ", " << eventsWithTrackReconstructedYet[i];
 352     }
 353     cout << "))" << endl;
 354   } else cout << ")" << endl;
 355   cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
 356   cout << "  (including " << nTotConnectedTracks << " track(s) still connected to a reconstructible one)" << endl;
 357   cout << "  (including " << nTotAdditionalTracks << " additional track(s) (compared to the number of expected ones))" << endl;
 358   cout << "- Number of events with fake track(s): " << nEventsWithFake;
 359   if (nEventsWithFake > 0) {
 360     for(Int_t i=0; i<nEventsWithFake; i++){
 361       if (i==0) cout << " (eventID = " << eventsWithFake[i];
 362       else cout << ", " << eventsWithFake[i];
 363     }
 364     cout << ")" << endl;
 365   } else cout << endl;
 366   cout << "  (including " << nEventsWithAdditionalFake << " events with additional track(s)";
 367   if (nEventsWithAdditionalFake > 0) {
 368     for(Int_t i=0; i<nEventsWithAdditionalFake; i++){
 369       if (i==0) cout << " (eventID = " << eventsWithAdditionalFake[i];
 370       else cout << ", " << eventsWithAdditionalFake[i];
 371     }
 372     cout << "))" << endl;
 373   } else cout << ")" << endl;
 374   cout << endl;
 375   cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
 376   cout << endl;
 377
 378 }
 379
 380 //-----------------------------------------------------------------------
 381 Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
 382                            Bool_t useLabel, TH1F &hFractionOfConnectedClusters)
 383 {
 384   /// loop over reconstructible TrackRef not associated with reconstructed track:
 385   /// for each of them, find and remove the most connected the fake track, if any,
 386   /// and fill the histograms with the fraction of connected clusters.
 387   /// Return the number of reconstructible track not connected to any fake
 388
 389   Int_t nFreeMissingTracks = 0;
 390
 391   // loop over trackRefs
 392   TIter next(trackRefStore.CreateIterator());
 393   AliMUONTrack* trackRef;
 394   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
 395
 396     // skip not reconstructible trackRefs
 397     if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue;
 398
 399     Int_t label = trackRef->GetUniqueID();
 400
 401     // look for the most connected fake track
 402     AliMUONTrack *connectedFake = 0x0;
 403     Double_t fractionOfConnectedClusters = 0.;
 404     TIter next2(fakeTrackStore.CreateIterator());
 405     AliMUONTrack* fakeTrack;
 406     while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
 407
 408       // get the number of connected clusters
 409       Int_t nConnectedClusters = 0;
 410       if (useLabel) { // by using the MC label
 411         for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
 412           if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
 413             nConnectedClusters++;
 414       } else { // by comparing cluster/TrackRef positions
 415         Bool_t compTrack[10];
 416         nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, sigmaCut, compTrack);
 417       }
 418
 419       // skip non-connected fake tracks
 420       if (nConnectedClusters == 0) continue;
 421
 422       // check if it is the most connected fake track
 423       Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters());
 424       if (f > fractionOfConnectedClusters) {
 425         connectedFake = fakeTrack;
 426         fractionOfConnectedClusters = f;
 427       }
 428
 429     }
 430
 431     // remove the most connected fake track
 432     if (connectedFake) {
 433       hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters);
 434       fakeTrackStore.Remove(*connectedFake);
 435     } else nFreeMissingTracks++;
 436
 437   }
 438
 439   return nFreeMissingTracks;
 440
 441 }
 442