PWG3/muondep/AliAnalysisTaskMuonFakes.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // ROOT includes
  19 #include <TObjArray.h>
  20 #include <TClonesArray.h>
  21 #include <TH1F.h>
  22 #include <TH2F.h>
  23 #include <TCanvas.h>
  24
  25 // STEER includes
  26 #include "AliLog.h"
  27 #include "AliESDEvent.h"
  28 #include "AliESDMuonTrack.h"
  29 #include "AliESDInputHandler.h"
  30 #include "AliMCEventHandler.h"
  31 #include "AliCDBManager.h"
  32 #include "AliCentrality.h"
  33
  34 // ANALYSIS includes
  35 #include "AliAnalysisDataSlot.h"
  36 #include "AliAnalysisDataContainer.h"
  37 #include "AliAnalysisManager.h"
  38
  39 // MUON includes
  40 #include "AliMUONCDB.h"
  41 #include "AliMUONRecoParam.h"
  42 #include "AliMUONRecoCheck.h"
  43 #include "AliMUONVCluster.h"
  44 #include "AliMUONVTrackStore.h"
  45 #include "AliMUONVTriggerTrackStore.h"
  46 #include "AliMUONTrack.h"
  47 #include "AliMUONTrackParam.h"
  48 #include "AliMUONESDInterface.h"
  49
  50 #include "AliAnalysisTaskMuonFakes.h"
  51 #include "AliCounterCollection.h"
  52
  53 ClassImp(AliAnalysisTaskMuonFakes)
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskMuonFakes::AliAnalysisTaskMuonFakes() :
  57 AliAnalysisTaskSE(),
  58 fList(0x0),
  59 fList2(0x0),
  60 fCanvases(0x0),
  61 fTrackCounters(0x0),
  62 fFakeTrackCounters(0x0),
  63 fMatchedTrackCounters(0x0),
  64 fEventCounters(0x0),
  65 fPairCounters(0x0),
  66 fCurrentFileName(""),
  67 fRequestedStationMask(0),
  68 fRequest2ChInSameSt45(kFALSE),
  69 fSigmaCut(-1.),
  70 fUseLabel(kFALSE),
  71 fMatchTrig(kFALSE),
  72 fApplyAccCut(kFALSE),
  73 fRecoParamLocation("alien://folder=/alice/simulation/2008/v4-15-Release/Full")
  74 {
  75   /// Default constructor.
  76 }
  77
  78 //________________________________________________________________________
  79 AliAnalysisTaskMuonFakes::AliAnalysisTaskMuonFakes(const char *name) :
  80 AliAnalysisTaskSE(name),
  81 fList(0x0),
  82 fList2(0x0),
  83 fCanvases(0x0),
  84 fTrackCounters(0x0),
  85 fFakeTrackCounters(0x0),
  86 fMatchedTrackCounters(0x0),
  87 fEventCounters(0x0),
  88 fPairCounters(0x0),
  89 fCurrentFileName(""),
  90 fRequestedStationMask(0),
  91 fRequest2ChInSameSt45(kFALSE),
  92 fSigmaCut(-1.),
  93 fUseLabel(kFALSE),
  94 fMatchTrig(kFALSE),
  95 fApplyAccCut(kFALSE),
  96 fRecoParamLocation("alien://folder=/alice/simulation/2008/v4-15-Release/Full")
  97 {
  98   /// Constructor.
  99   // Output slot #1 writes into a TObjArray container
 100   DefineOutput(1,TObjArray::Class());
 101   // Output slot #2 writes into an AliCounterCollection container
 102   DefineOutput(2,AliCounterCollection::Class());
 103   // Output slot #3 writes into an AliCounterCollection container
 104   DefineOutput(3,AliCounterCollection::Class());
 105   // Output slot #4 writes into an AliCounterCollection container
 106   DefineOutput(4,AliCounterCollection::Class());
 107   // Output slot #5 writes into an AliCounterCollection container
 108   DefineOutput(5,AliCounterCollection::Class());
 109   // Output slot #6 writes into a TObjArray container
 110   DefineOutput(6,TObjArray::Class());
 111   // Output slot #7 writes into an AliCounterCollection container
 112   DefineOutput(7,AliCounterCollection::Class());
 113 }
 114
 115 //________________________________________________________________________
 116 AliAnalysisTaskMuonFakes::~AliAnalysisTaskMuonFakes()
 117 {
 118   /// Destructor.
 119   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 120     delete fList;
 121     delete fList2;
 122     delete fTrackCounters;
 123     delete fFakeTrackCounters;
 124     delete fMatchedTrackCounters;
 125     delete fEventCounters;
 126     delete fPairCounters;
 127   }
 128   delete fCanvases;
 129 }
 130
 131 //___________________________________________________________________________
 132 void AliAnalysisTaskMuonFakes::UserCreateOutputObjects()
 133 {
 134   /// Create histograms and counters.
 135
 136   // single track histograms
 137   fList = new TObjArray(100);
 138   fList->SetOwner();
 139
 140   // number of tracks
 141   TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks", "nb of tracks /evt", 21, -0.5, 20.5);
 142   fList->AddAtAndExpand(hNumberOfTracks, kNumberOfTracks);
 143   TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks", "nb of fake - nb of missing track", 21, -0.5, 20.5);
 144   fList->AddAtAndExpand(hNumberOfAdditionalTracks, kNumberOfAdditionalTracks);
 145
 146   // number of clusters
 147   TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters", "nb of clusters /track", 21, -0.5, 20.5);
 148   fList->AddAtAndExpand(hNumberOfClusters, kNumberOfClusters);
 149   TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM", "nb of clusters /matched track", 21, -0.5, 20.5);
 150   fList->AddAtAndExpand(hNumberOfClustersM, kNumberOfClustersM);
 151   TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF", "nb of clusters /fake track", 21, -0.5, 20.5);
 152   fList->AddAtAndExpand(hNumberOfClustersF, kNumberOfClustersF);
 153   TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC", "nb of clusters /MC track", 21, -0.5, 20.5);
 154   fList->AddAtAndExpand(hNumberOfClustersMC, kNumberOfClustersMC);
 155   TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters", "nb of matched clusters / nb of clusters", 110, 0., 1.1);
 156   fList->AddAtAndExpand(hFractionOfMatchedClusters, kFractionOfMatchedClusters);
 157   TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters", "nb of connected clusters / nb of clusters in fake tracks", 110, 0., 1.1);
 158   fList->AddAtAndExpand(hFractionOfConnectedClusters, kFractionOfConnectedClusters);
 159
 160   // number of fired chambers
 161   TH1F *hNumberOfChamberHit = new TH1F("hNumberOfChamberHit", "nb of chambers hit /track", 16, -0.5, 15.5);
 162   fList->AddAtAndExpand(hNumberOfChamberHit, kNumberOfChamberHit);
 163   TH1F *hNumberOfChamberHitM = new TH1F("hNumberOfChamberHitM", "nb of chambers hit /matched track", 16, -0.5, 15.5);
 164   fList->AddAtAndExpand(hNumberOfChamberHitM, kNumberOfChamberHitM);
 165   TH1F *hNumberOfChamberHitF = new TH1F("hNumberOfChamberHitF", "nb of chambers hit /fake track", 16, -0.5, 15.5);
 166   fList->AddAtAndExpand(hNumberOfChamberHitF, kNumberOfChamberHitF);
 167
 168   // chi2
 169   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 200, 0., 20.);
 170   fList->AddAtAndExpand(hChi2PerDof, kChi2PerDof);
 171   TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 200, 0., 20.);
 172   fList->AddAtAndExpand(hChi2PerDofM, kChi2PerDofM);
 173   TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 200, 0., 20.);
 174   fList->AddAtAndExpand(hChi2PerDofF, kChi2PerDofF);
 175
 176   // chi2 versus number of clusters
 177   TH2F *hChi2PerDofVsNClusters = new TH2F("hChi2PerDofVsNClusters", "track chi2/d.o.f. versus nb of clusters", 21, -0.5, 20.5, 100, 0., 20.);
 178   fList->AddAtAndExpand(hChi2PerDofVsNClusters, kChi2PerDofVsNClusters);
 179   TH2F *hChi2PerDofVsNClustersM = new TH2F("hChi2PerDofVsNClustersM", "matched track chi2/d.o.f. versus nb of clusters", 21, -0.5, 20.5, 100, 0., 20.);
 180   fList->AddAtAndExpand(hChi2PerDofVsNClustersM, kChi2PerDofVsNClustersM);
 181   TH2F *hChi2PerDofVsNClustersF = new TH2F("hChi2PerDofVsNClustersF", "fake track chi2/d.o.f. versus nb of clusters", 21, -0.5, 20.5, 100, 0., 20.);
 182   fList->AddAtAndExpand(hChi2PerDofVsNClustersF, kChi2PerDofVsNClustersF);
 183
 184   // chi2 versus number of fired chambers
 185   TH2F *hChi2PerDofVsNChamberHit = new TH2F("hChi2PerDofVsNChamberHit", "track chi2/d.o.f. versus nb of fired chambers", 16, -0.5, 15.5, 100, 0., 20.);
 186   fList->AddAtAndExpand(hChi2PerDofVsNChamberHit, kChi2PerDofVsNChamberHit);
 187   TH2F *hChi2PerDofVsNChamberHitM = new TH2F("hChi2PerDofVsNChamberHitM", "matched track chi2/d.o.f. versus nb of fired chambers", 16, -0.5, 15.5, 100, 0., 20.);
 188   fList->AddAtAndExpand(hChi2PerDofVsNChamberHitM, kChi2PerDofVsNChamberHitM);
 189   TH2F *hChi2PerDofVsNChamberHitF = new TH2F("hChi2PerDofVsNChamberHitF", "fake track chi2/d.o.f. versus nb of fired chambers", 16, -0.5, 15.5, 100, 0., 20.);
 190   fList->AddAtAndExpand(hChi2PerDofVsNChamberHitF, kChi2PerDofVsNChamberHitF);
 191
 192   // physics quantities
 193   TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
 194   fList->AddAtAndExpand(hP, kP);
 195   TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
 196   fList->AddAtAndExpand(hPM, kPM);
 197   TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
 198   fList->AddAtAndExpand(hPF, kPF);
 199   TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
 200   fList->AddAtAndExpand(hPt, kPt);
 201   TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
 202   fList->AddAtAndExpand(hPtM, kPtM);
 203   TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
 204   fList->AddAtAndExpand(hPtF, kPtF);
 205   TH1F *hEta = new TH1F("hEta", "Muon pseudo-rapidity distribution", 200, -10., 0.);
 206   fList->AddAtAndExpand(hEta , kEta );
 207   TH1F *hEtaM = new TH1F("hEtaM", "matched track pseudo-rapidity distribution", 200, -10., 0.);
 208   fList->AddAtAndExpand(hEtaM, kEtaM);
 209   TH1F *hEtaF = new TH1F("hEtaF", "fake track pseudo-rapidity distribution", 200, -10., 0.);
 210   fList->AddAtAndExpand(hEtaF, kEtaF);
 211   TH1F *hPhi = new TH1F("hPhi", "Muon phi distribution", 100, -1., 9.);
 212   fList->AddAtAndExpand(hPhi, kPhi);
 213   TH1F *hPhiM = new TH1F("hPhiM", "matched track phi distribution", 100, -1., 9.);
 214   fList->AddAtAndExpand(hPhiM, kPhiM);
 215   TH1F *hPhiF = new TH1F("hPhiF", "fake track phi distribution", 100, -1., 9.);
 216   fList->AddAtAndExpand(hPhiF, kPhiF);
 217   TH1F *hDCA = new TH1F("hDCA", "Muon DCA distribution", 250, 0., 500.);
 218   fList->AddAtAndExpand(hDCA, kDCA);
 219   TH1F *hDCAM = new TH1F("hDCAM", "matched track DCA distribution", 250, 0., 500.);
 220   fList->AddAtAndExpand(hDCAM, kDCAM);
 221   TH1F *hDCAF = new TH1F("hDCAF", "fake track DCA distribution", 250, 0., 500.);
 222   fList->AddAtAndExpand(hDCAF, kDCAF);
 223   TH1F *hRAbs = new TH1F("hRAbs", "Muon R_{Abs} distribution", 200, 0., 100.);
 224   fList->AddAtAndExpand(hRAbs, kRAbs);
 225   TH1F *hRAbsM = new TH1F("hRAbsM", "matched track R_{Abs} distribution", 300, 0., 150.);
 226   fList->AddAtAndExpand(hRAbsM, kRAbsM);
 227   TH1F *hRAbsF = new TH1F("hRAbsF", "fake track R_{Abs} distribution", 200, 0., 100.);
 228   fList->AddAtAndExpand(hRAbsF, kRAbsF);
 229
 230   // track pair histograms
 231   fList2 = new TObjArray(100);
 232   fList2->SetOwner();
 233
 234   // physics quantities of opposite-sign track pairs
 235   TH1F* h = 0x0;
 236   TString nameSuffix[4] = {"", "M", "F1", "F2"};
 237   TString titlePrefix[4] = {"Dimuon", "matched-matched pair", "matched-fake pair", "fake-fake pair"};
 238   for (Int_t i = 0; i < 4; i++) {
 239     h = new TH1F(Form("h2Mass%s",nameSuffix[i].Data()), Form("%s mass distribution (GeV/c^{2})",titlePrefix[i].Data()), 300, 0., 15.);
 240     fList2->AddAtAndExpand(h, k2Mass+i);
 241     h = new TH1F(Form("h2P%s",nameSuffix[i].Data()), Form("%s P distribution (GeV/c)",titlePrefix[i].Data()), 100, 0., 200.);
 242     fList2->AddAtAndExpand(h, k2P+i);
 243     h = new TH1F(Form("h2Pt%s",nameSuffix[i].Data()), Form("%s Pt distribution (GeV/c)",titlePrefix[i].Data()), 100, 0., 20.);
 244     fList2->AddAtAndExpand(h, k2Pt+i);
 245     h = new TH1F(Form("h2Y%s",nameSuffix[i].Data()), Form("%s rapidity distribution",titlePrefix[i].Data()), 200, -10., 0.);
 246     fList2->AddAtAndExpand(h, k2Y+i);
 247     h = new TH1F(Form("h2Eta%s",nameSuffix[i].Data()), Form("%s pseudo-rapidity distribution",titlePrefix[i].Data()), 200, -10., 0.);
 248     fList2->AddAtAndExpand(h, k2Eta+i);
 249     h = new TH1F(Form("h2Phi%s",nameSuffix[i].Data()), Form("%s phi distribution",titlePrefix[i].Data()), 100, -1., 9.);
 250     fList2->AddAtAndExpand(h, k2Phi+i);
 251   }
 252
 253   // global counters of tracks:
 254   // - reconstructible = number of reconstructible tracks
 255   // - reconstructed   = number of reconstructed tracks
 256   // - matched         = number of reconstructed tracks matched with a simulated one (reconstructible or not)
 257   // - matchedyet      = number of reconstructed tracks matched with a simulated one that is not reconstructible
 258   // - fake            = number of fake tracks
 259   // - connected       = number of fake tracks connected to a reconstructible simulated track
 260   // - additional      = number of additional (fake) tracks compared to the number of reconstructible ones
 261   fTrackCounters = new AliCounterCollection(GetOutputSlot(2)->GetContainer()->GetName());
 262   fTrackCounters->AddRubric("track", "reconstructible/reconstructed/matched/matchedyet/fake/connected/additional");
 263   fTrackCounters->AddRubric("run", 1000000);
 264   fTrackCounters->AddRubric("trig", "yes/no/unknown");
 265   fTrackCounters->AddRubric("selected", "yes/no");
 266   fTrackCounters->AddRubric("acc", "in/out/unknown");
 267   TString centralityClasses = "5/10/15/20/25/30/35/40/45/50/55/60/65/70/75/80/85/90/95/100";
 268   fTrackCounters->AddRubric("cent", centralityClasses.Data());
 269   fTrackCounters->Init();
 270
 271   // detailled counters of fake tracks:
 272   fFakeTrackCounters = new AliCounterCollection(GetOutputSlot(3)->GetContainer()->GetName());
 273   fFakeTrackCounters->AddRubric("track", "fake/connected/additional/matchedyet/fake?");
 274   fFakeTrackCounters->AddRubric("run", 1000000);
 275   fFakeTrackCounters->AddRubric("file", 1000000);
 276   fFakeTrackCounters->AddRubric("event", 1000000);
 277   fFakeTrackCounters->AddRubric("trig", "yes/no/unknown");
 278   fFakeTrackCounters->AddRubric("selected", "yes/no");
 279   fFakeTrackCounters->AddRubric("acc", "in/out/unknown");
 280   fFakeTrackCounters->AddRubric("cent", centralityClasses.Data());
 281   fFakeTrackCounters->Init();
 282
 283   // counters of tracks matched by position or by using MC labels
 284   fMatchedTrackCounters = new AliCounterCollection(GetOutputSlot(4)->GetContainer()->GetName());
 285   fMatchedTrackCounters->AddRubric("position", "match/not match");
 286   fMatchedTrackCounters->AddRubric("label", "match/not match/match other");
 287   fMatchedTrackCounters->AddRubric("run", 1000000);
 288   fMatchedTrackCounters->AddRubric("trig", "yes/no");
 289   fMatchedTrackCounters->AddRubric("selected", "yes/no");
 290   fMatchedTrackCounters->AddRubric("acc", "in/out");
 291   fMatchedTrackCounters->AddRubric("cent", centralityClasses.Data());
 292   fMatchedTrackCounters->Init();
 293
 294   // global counters of events
 295   // - any             = total number of events with reconstructed tracks
 296   // - fake            = number of events with fake track(s)
 297   // - notconnected    = number of events with fake tracks that are not connected to a reconstructible simulated track
 298   // - additional      = number of events with additional (fake) tracks compared to the number of reconstructible ones
 299   // - matchedyet      = number of events with reconstructed tracks matched with a simulated one that is not reconstructible
 300   // if trig = yes: only the tracks matched with the trigger are considered in the above logic
 301   fEventCounters = new AliCounterCollection(GetOutputSlot(5)->GetContainer()->GetName());
 302   fEventCounters->AddRubric("event", "any/fake/notconnected/additional/matchedyet");
 303   fEventCounters->AddRubric("run", 1000000);
 304   fEventCounters->AddRubric("trig", "any/yes");
 305   fEventCounters->AddRubric("selected", "yes/no");
 306   fEventCounters->AddRubric("cent", centralityClasses.Data());
 307   fEventCounters->Init();
 308
 309   // global counters of track pairs:
 310   // - reconstructible = number of reconstructible track pairs
 311   // - reconstructed   = number of reconstructed track pairs
 312   // - matched         = number of reconstructed track pairs fully matched with a simulated one (reconstructible or not)
 313   // - 1fake           = number of reconstructed track pairs made of one matched track and one fake
 314   // - 2fakes          = number of reconstructed track pairs fully fake
 315   fPairCounters = new AliCounterCollection(GetOutputSlot(7)->GetContainer()->GetName());
 316   fPairCounters->AddRubric("pair", "reconstructible/reconstructed/matched/1fake/2fakes");
 317   fPairCounters->AddRubric("run", 1000000);
 318   fPairCounters->AddRubric("trig", "0/1/2");
 319   fPairCounters->AddRubric("selected", "yes/no");
 320   fPairCounters->AddRubric("acc", "in/out/unknown");
 321   fPairCounters->AddRubric("cent", centralityClasses.Data());
 322   fPairCounters->Init();
 323
 324   // Disable printout of AliMCEvent
 325   AliLog::SetClassDebugLevel("AliMCEvent",-1);
 326
 327   // Post data at least once per task to ensure data synchronisation (required for merging)
 328   PostData(1, fList);
 329   PostData(2, fTrackCounters);
 330   PostData(3, fFakeTrackCounters);
 331   PostData(4, fMatchedTrackCounters);
 332   PostData(5, fEventCounters);
 333   PostData(6, fList2);
 334   PostData(7, fPairCounters);
 335 }
 336
 337 //________________________________________________________________________
 338 void AliAnalysisTaskMuonFakes::UserExec(Option_t *)
 339 {
 340   /// Process event: looks for fakes...
 341
 342   // check that reconstructions parameters for that run have been properly set
 343   if (fSigmaCut < 0) return;
 344
 345   // check physics selection
 346   TString selected = (fInputHandler && fInputHandler->IsEventSelected() != 0) ? "selected:yes" : "selected:no";
 347
 348   // current file name
 349   fCurrentFileName = CurrentFileName();
 350   fCurrentFileName.ReplaceAll("alien://","");
 351   fCurrentFileName.ReplaceAll("/","\\");
 352   fCurrentFileName.ReplaceAll(":",";");
 353
 354   // Load ESD event
 355   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 356   if (!esd) {
 357     AliError("Cannot get input event");
 358     return;
 359   }
 360
 361   // current centrality class
 362   TString centrality = "cent:";
 363   Double_t centralityValue = esd->GetCentrality()->GetCentralityPercentile("V0M");
 364   TObjArray* centralylimits = fTrackCounters->GetKeyWords("cent").Tokenize(",");
 365   TObjString* limit = 0x0;
 366   TIter nextLimit(centralylimits);
 367   while ((limit = static_cast<TObjString*>(nextLimit()))) {
 368     if (centralityValue < limit->String().Atoi()) {
 369       centrality += limit->GetName();
 370       break;
 371     }
 372   }
 373   if (!limit) centrality += static_cast<TObjString*>(centralylimits->Last())->GetName();
 374   delete centralylimits;
 375
 376   // Load MC event
 377   AliMCEventHandler *mcH = 0;
 378   if(MCEvent()) mcH = static_cast<AliMCEventHandler*>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 379
 380   // get reconstructed and simulated tracks
 381   AliMUONRecoCheck rc(esd,mcH);
 382   AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(-1, kFALSE);
 383   AliMUONVTrackStore* trackRefStore = rc.TrackRefs(-1);
 384   AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(-1);
 385   if (!muonTrackStore || !trackRefStore) return;
 386
 387   // loop over trackRefs
 388   Int_t nMuPlus[2] = {0, 0};
 389   Int_t nMuMinus[2] = {0, 0};
 390   TIter next(trackRefStore->CreateIterator());
 391   AliMUONTrack* trackRef;
 392   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
 393
 394     // skip trackRefs that are not reconstructible
 395     if (!trackRef->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) continue;
 396
 397     // trigger condition
 398     Bool_t trigger = (triggerTrackRefStore->FindObject(trackRef->GetUniqueID()));
 399     Int_t iTrig = trigger ? 1 : 0;
 400     TString trig = trigger ? "trig:yes" : "trig:no";
 401
 402     // count muons
 403     if (trackRef->GetTrackParamAtVertex()->GetCharge() > 0) nMuPlus[iTrig]++;
 404     else nMuMinus[iTrig]++;
 405
 406     // fill global counters
 407     fTrackCounters->Count(Form("track:reconstructible/run:%d/%s/%s/acc:unknown/%s", fCurrentRunNumber, trig.Data(), selected.Data(), centrality.Data()));
 408
 409   }
 410
 411   // fill global counters
 412   fPairCounters->Count(Form("pair:reconstructible/run:%d/trig:0/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nMuPlus[0]*nMuMinus[0]);
 413   fPairCounters->Count(Form("pair:reconstructible/run:%d/trig:1/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nMuPlus[1]*nMuMinus[0]+nMuPlus[0]*nMuMinus[1]);
 414   fPairCounters->Count(Form("pair:reconstructible/run:%d/trig:2/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nMuPlus[1]*nMuMinus[1]);
 415
 416   // loop over ESD tracks
 417   Int_t nTrackerTracks = 0;
 418   Bool_t containTrack[2] = {kFALSE, kFALSE};
 419   Bool_t containFakeTrack[2] = {kFALSE, kFALSE};
 420   Bool_t containMatchedYetTrack[2] = {kFALSE, kFALSE};
 421   AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
 422   Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
 423   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 424
 425     AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
 426
 427     // skip ghosts
 428     if (!esdTrack->ContainTrackerData()) continue;
 429     containTrack[0] = kTRUE;
 430
 431     // trigger condition
 432     Bool_t trigger = esdTrack->ContainTriggerData();
 433     TString trig = trigger ? "trig:yes" : "trig:no";
 434     if (trigger) containTrack[1] = kTRUE;
 435
 436     // acceptance condition
 437     Double_t rAbs = esdTrack->GetRAtAbsorberEnd();
 438     Double_t thetaTrackAbsEnd = TMath::ATan(rAbs/505.) * TMath::RadToDeg();
 439     Double_t eta = esdTrack->Eta();
 440     Bool_t inAcc = (thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 10. && eta >= -4. && eta <= -2.5);
 441     TString acc = inAcc ? "acc:in" : "acc:out";
 442
 443     // fill global counters
 444     if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) nTrackerTracks++;
 445     fTrackCounters->Count(Form("track:reconstructed/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 446
 447     // find the corresponding MUON track
 448     AliMUONTrack* muonTrack = static_cast<AliMUONTrack*>(muonTrackStore->FindObject(esdTrack->GetUniqueID()));
 449
 450     // get track info
 451     Int_t nClusters = esdTrack->GetNClusters();
 452     Int_t nChamberHit = 0;
 453     for (Int_t ich=0; ich<10; ich++) if (esdTrack->IsInMuonClusterMap(ich)) nChamberHit++;
 454     Double_t normalizedChi2 = esdTrack->GetChi2() / (2. * esdTrack->GetNHit() - 5);
 455     Double_t p = esdTrack->P();
 456     Double_t pT = esdTrack->Pt();
 457     Double_t phi = esdTrack->Phi();
 458     Double_t dca = esdTrack->GetDCA();
 459
 460     // fill global histograms
 461     if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) {
 462       ((TH1F*)fList->UncheckedAt(kNumberOfClusters))->Fill(nClusters);
 463       ((TH1F*)fList->UncheckedAt(kNumberOfChamberHit))->Fill(nChamberHit);
 464       ((TH1F*)fList->UncheckedAt(kChi2PerDof))->Fill(normalizedChi2);
 465       ((TH1F*)fList->UncheckedAt(kP))->Fill(p);
 466       ((TH1F*)fList->UncheckedAt(kPt))->Fill(pT);
 467       ((TH1F*)fList->UncheckedAt(kEta))->Fill(eta);
 468       ((TH1F*)fList->UncheckedAt(kPhi))->Fill(phi);
 469       ((TH1F*)fList->UncheckedAt(kDCA))->Fill(dca);
 470       ((TH1F*)fList->UncheckedAt(kRAbs))->Fill(rAbs);
 471       ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNClusters))->Fill(nClusters,normalizedChi2);
 472       ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNChamberHit))->Fill(nChamberHit,normalizedChi2);
 473     }
 474
 475     // try to match, by position, the reconstructed track with a simulated one
 476     Int_t nMatchClustersByPosition = 0;
 477     AliMUONTrack* matchedTrackRefByPosition = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClustersByPosition, kFALSE, fSigmaCut);
 478     Int_t MCLabelByPosition = (matchedTrackRefByPosition) ? static_cast<Int_t>(matchedTrackRefByPosition->GetUniqueID()) : -1;
 479
 480     // try to match, by using MC labels, the reconstructed track with a simulated one
 481     Int_t nMatchClustersByLabel = 0;
 482     AliMUONTrack* matchedTrackRefByLabel = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClustersByLabel, kTRUE, fSigmaCut);
 483     Int_t MCLabelByLabel = (matchedTrackRefByLabel) ? static_cast<Int_t>(matchedTrackRefByLabel->GetUniqueID()) : -1;
 484
 485     // fill global counters
 486     TString positionCase = (MCLabelByPosition >= 0) ? "position:match" : "position:not match";
 487     TString labelCase = "label:";
 488     if (MCLabelByLabel >= 0 && MCLabelByPosition >= 0 && MCLabelByLabel != MCLabelByPosition) labelCase += "match other";
 489     else if (MCLabelByLabel >= 0) labelCase += "match";
 490     else labelCase += "not match";
 491     fMatchedTrackCounters->Count(Form("%s/%s/run:%d/%s/%s/%s/%s", positionCase.Data(), labelCase.Data(), fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 492     if ((MCLabelByLabel >= 0 && MCLabelByPosition < 0) || (MCLabelByLabel < 0 && MCLabelByPosition >= 0))
 493       fFakeTrackCounters->Count(Form("track:fake?/run:%d/file:%s/event:%d/%s/%s/%s/%s", fCurrentRunNumber, fCurrentFileName.Data(),
 494                                      esd->GetEventNumberInFile(), trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 495
 496     // take actions according to the matching result we are interested in
 497     Int_t nMatchClusters = (fUseLabel) ? nMatchClustersByLabel : nMatchClustersByPosition;
 498     AliMUONTrack* matchedTrackRef = (fUseLabel) ? matchedTrackRefByLabel : matchedTrackRefByPosition;
 499     if (matchedTrackRef) {
 500
 501       // fill global counters
 502       fTrackCounters->Count(Form("track:matched/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 503
 504       // track matched with a trackRef that is not reconstructible
 505       if (!matchedTrackRef->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) {
 506
 507         containMatchedYetTrack[0] = kTRUE;
 508         if (trigger) containMatchedYetTrack[1] = kTRUE;
 509
 510         // fill global counters
 511         fTrackCounters->Count(Form("track:matchedyet/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 512         fFakeTrackCounters->Count(Form("track:matchedyet/run:%d/file:%s/event:%d/%s/%s/%s/%s", fCurrentRunNumber, fCurrentFileName.Data(),
 513                                        esd->GetEventNumberInFile(), trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 514       }
 515
 516       // fill histograms
 517       if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) {
 518         if (nClusters > 0) ((TH1F*)fList->UncheckedAt(kFractionOfMatchedClusters))->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
 519         ((TH1F*)fList->UncheckedAt(kNumberOfClustersMC))->Fill(matchedTrackRef->GetNClusters());
 520         ((TH1F*)fList->UncheckedAt(kNumberOfClustersM))->Fill(nClusters);
 521         ((TH1F*)fList->UncheckedAt(kNumberOfChamberHitM))->Fill(nChamberHit);
 522         ((TH1F*)fList->UncheckedAt(kChi2PerDofM))->Fill(normalizedChi2);
 523         ((TH1F*)fList->UncheckedAt(kPM))->Fill(p);
 524         ((TH1F*)fList->UncheckedAt(kPtM))->Fill(pT);
 525         ((TH1F*)fList->UncheckedAt(kEtaM))->Fill(eta);
 526         ((TH1F*)fList->UncheckedAt(kPhiM))->Fill(phi);
 527         ((TH1F*)fList->UncheckedAt(kDCAM))->Fill(dca);
 528         ((TH1F*)fList->UncheckedAt(kRAbsM))->Fill(rAbs);
 529         ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNClustersM))->Fill(nClusters,normalizedChi2);
 530         ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNChamberHitM))->Fill(nChamberHit,normalizedChi2);
 531       }
 532
 533       // flag matched tracks
 534       esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
 535
 536       // remove already matched trackRefs
 537       trackRefStore->Remove(*matchedTrackRef);
 538
 539     } else {
 540
 541       containFakeTrack[0] = kTRUE;
 542       if (trigger) containFakeTrack[1] = kTRUE;
 543
 544       // fill global counters
 545       fTrackCounters->Count(Form("track:fake/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 546       fFakeTrackCounters->Count(Form("track:fake/run:%d/file:%s/event:%d/%s/%s/%s/%s", fCurrentRunNumber, fCurrentFileName.Data(),
 547                                      esd->GetEventNumberInFile(), trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 548
 549       // fill histograms
 550       if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) {
 551         ((TH1F*)fList->UncheckedAt(kNumberOfClustersF))->Fill(nClusters);
 552         ((TH1F*)fList->UncheckedAt(kNumberOfChamberHitF))->Fill(nChamberHit);
 553         ((TH1F*)fList->UncheckedAt(kChi2PerDofF))->Fill(normalizedChi2);
 554         ((TH1F*)fList->UncheckedAt(kPF))->Fill(p);
 555         ((TH1F*)fList->UncheckedAt(kPtF))->Fill(pT);
 556         ((TH1F*)fList->UncheckedAt(kEtaF))->Fill(eta);
 557         ((TH1F*)fList->UncheckedAt(kPhiF))->Fill(phi);
 558         ((TH1F*)fList->UncheckedAt(kDCAF))->Fill(dca);
 559         ((TH1F*)fList->UncheckedAt(kRAbsF))->Fill(rAbs);
 560         ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNClustersF))->Fill(nClusters,normalizedChi2);
 561         ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNChamberHitF))->Fill(nChamberHit,normalizedChi2);
 562       }
 563
 564       // flag fake tracks
 565       esdTrack->SetLabel(-1);
 566
 567       // store fake tracks
 568       fakeTrackStore->Add(*muonTrack);
 569
 570     }
 571
 572   } // end of loop over ESD tracks
 573
 574   // fill histogram and global counters
 575   ((TH1F*)fList->UncheckedAt(kNumberOfTracks))->Fill(nTrackerTracks);
 576   if (containTrack[0]) fEventCounters->Count(Form("event:any/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 577   if (containTrack[1]) fEventCounters->Count(Form("event:any/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 578   if (containFakeTrack[0]) fEventCounters->Count(Form("event:fake/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 579   if (containFakeTrack[1]) fEventCounters->Count(Form("event:fake/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 580   if (containMatchedYetTrack[0]) fEventCounters->Count(Form("event:matchedyet/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 581   if (containMatchedYetTrack[1]) fEventCounters->Count(Form("event:matchedyet/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 582
 583   // count the number of not connected and additional fake tracks
 584   if (fakeTrackStore->GetSize() > 0) {
 585
 586     // remove the most connected fake tracks
 587     Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, selected, centrality);
 588
 589     if (fakeTrackStore->GetSize() > 0) {
 590
 591       // fill global counters
 592       fEventCounters->Count(Form("event:notconnected/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 593
 594       // check status of remaining fakes with respect to the matching with trigger
 595       Bool_t containMatchedFake = kFALSE;
 596       Bool_t containUnmatchedFake = kFALSE;
 597       AliMUONTrack* fakeTrack = 0x0;
 598       TIter next3(fakeTrackStore->CreateIterator());
 599       while ( ( fakeTrack = static_cast<AliMUONTrack*>(next3()) ) ) {
 600         if (fakeTrack->GetMatchTrigger() > 0) containMatchedFake = kTRUE;
 601         else containUnmatchedFake = kTRUE;
 602       }
 603
 604       // fill global counters
 605       if (containMatchedFake) fEventCounters->Count(Form("event:notconnected/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 606
 607       // remove the remaining free reconstructible tracks
 608       Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
 609
 610       if (nAdditionalTracks > 0) {
 611
 612         // fill histogram and global counters
 613         ((TH1F*)fList->UncheckedAt(kNumberOfAdditionalTracks))->Fill(nAdditionalTracks);
 614         fEventCounters->Count(Form("event:additional/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 615         if (!containUnmatchedFake) { // all matched
 616           fTrackCounters->Count(Form("track:additional/run:%d/trig:yes/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nAdditionalTracks);
 617           fFakeTrackCounters->Count(Form("track:additional/run:%d/file:%s/event:%d/trig:yes/%s/acc:unknown/%s", fCurrentRunNumber, fCurrentFileName.Data(), esd->GetEventNumberInFile(), selected.Data(), centrality.Data()), nAdditionalTracks);
 618           fEventCounters->Count(Form("event:additional/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 619         } else if (!containMatchedFake) { // none matched
 620           fTrackCounters->Count(Form("track:additional/run:%d/trig:no/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nAdditionalTracks);
 621           fFakeTrackCounters->Count(Form("track:additional/run:%d/file:%s/event:%d/trig:no/%s/acc:unknown/%s", fCurrentRunNumber, fCurrentFileName.Data(), esd->GetEventNumberInFile(), selected.Data(), centrality.Data()), nAdditionalTracks);
 622         } else { // mixed
 623           fTrackCounters->Count(Form("track:additional/run:%d/trig:unknown/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nAdditionalTracks);
 624           fFakeTrackCounters->Count(Form("track:additional/run:%d/file:%s/event:%d/trig:unknown/%s/acc:unknown/%s", fCurrentRunNumber, fCurrentFileName.Data(), esd->GetEventNumberInFile(), selected.Data(), centrality.Data()), nAdditionalTracks);
 625           fEventCounters->Count(Form("event:additional/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
 626         }
 627
 628       }
 629
 630     }
 631
 632   }
 633
 634   // clean memory
 635   delete fakeTrackStore;
 636
 637   // double loop over ESD tracks, build pairs and fill histograms and counters according to their label
 638   TLorentzVector vMu1, vMu2, vDiMu;
 639   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
 640     AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
 641
 642     // skip ghosts
 643     if (!muonTrack1->ContainTrackerData()) continue;
 644
 645     // get track info
 646     Bool_t trigger1 = muonTrack1->ContainTriggerData();
 647     Double_t thetaAbs1 = TMath::ATan(muonTrack1->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 648     Double_t eta1 = muonTrack1->Eta();
 649     Bool_t acc1 = (thetaAbs1 >= 2. && thetaAbs1 <= 10. && eta1 >= -4. && eta1 <= -2.5);
 650     Short_t charge1 = muonTrack1->Charge();
 651     Int_t label1 = muonTrack1->GetLabel();
 652     muonTrack1->LorentzP(vMu1);
 653
 654     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
 655       AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
 656
 657       // skip ghosts
 658       if (!muonTrack2->ContainTrackerData()) continue;
 659
 660       // keep only opposite sign pairs
 661       Short_t charge2 = muonTrack2->Charge();
 662       if (charge1*charge2 > 0) continue;
 663
 664       // get track info
 665       Bool_t trigger2 = muonTrack2->ContainTriggerData();
 666       Double_t thetaAbs2 = TMath::ATan(muonTrack2->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 667       Double_t eta2 = muonTrack2->Eta();
 668       Bool_t acc2 = (thetaAbs2 >= 2. && thetaAbs2 <= 10. && eta2 >= -4. && eta2 <= -2.5);
 669       Int_t label2 = muonTrack2->GetLabel();
 670       muonTrack2->LorentzP(vMu2);
 671
 672       // compute kinematics of the pair
 673       vDiMu = vMu1 + vMu2;
 674       Float_t mass = vDiMu.M();
 675       Float_t p = vDiMu.P();
 676       Float_t pt = vDiMu.Pt();
 677       Float_t y = vDiMu.Rapidity();
 678       Float_t eta = vDiMu.Eta();
 679       Float_t phi = vDiMu.Phi();
 680       if (phi < 0) phi += 2.*TMath::Pi();
 681
 682       // trigger condition
 683       TString trig = "trig:";
 684       if (trigger1 && trigger2) trig += "2";
 685       else if (trigger1 || trigger2) trig += "1";
 686       else trig += "0";
 687
 688       // acceptance condition
 689       Bool_t inAcc = (acc1 && acc2 && eta >= -4. && eta <= -2.5);
 690       TString acc = inAcc ? "acc:in" : "acc:out";
 691
 692       // fill global histograms
 693       if ((!fMatchTrig || (trigger1 && trigger2)) && (!fApplyAccCut || inAcc)) {
 694         ((TH1F*)fList2->UncheckedAt(k2Mass))->Fill(mass);
 695         ((TH1F*)fList2->UncheckedAt(k2P))->Fill(p);
 696         ((TH1F*)fList2->UncheckedAt(k2Pt))->Fill(pt);
 697         ((TH1F*)fList2->UncheckedAt(k2Y))->Fill(y);
 698         ((TH1F*)fList2->UncheckedAt(k2Eta))->Fill(eta);
 699         ((TH1F*)fList2->UncheckedAt(k2Phi))->Fill(phi);
 700       }
 701
 702       TString pair = "pair:";
 703
 704       // fill histograms according to labels
 705       if (label1 >= 0 && label2 >= 0) {
 706
 707         pair += "matched";
 708
 709         if ((!fMatchTrig || (trigger1 && trigger2)) && (!fApplyAccCut || inAcc)) {
 710           ((TH1F*)fList2->UncheckedAt(k2MassM))->Fill(mass);
 711           ((TH1F*)fList2->UncheckedAt(k2PM))->Fill(p);
 712           ((TH1F*)fList2->UncheckedAt(k2PtM))->Fill(pt);
 713           ((TH1F*)fList2->UncheckedAt(k2YM))->Fill(y);
 714           ((TH1F*)fList2->UncheckedAt(k2EtaM))->Fill(eta);
 715           ((TH1F*)fList2->UncheckedAt(k2PhiM))->Fill(phi);
 716         }
 717
 718       } else if (label1 >= 0 || label2 >= 0) {
 719
 720         pair += "1fake";
 721
 722         if ((!fMatchTrig || (trigger1 && trigger2)) && (!fApplyAccCut || inAcc)) {
 723           ((TH1F*)fList2->UncheckedAt(k2MassF1))->Fill(mass);
 724           ((TH1F*)fList2->UncheckedAt(k2PF1))->Fill(p);
 725           ((TH1F*)fList2->UncheckedAt(k2PtF1))->Fill(pt);
 726           ((TH1F*)fList2->UncheckedAt(k2YF1))->Fill(y);
 727           ((TH1F*)fList2->UncheckedAt(k2EtaF1))->Fill(eta);
 728           ((TH1F*)fList2->UncheckedAt(k2PhiF1))->Fill(phi);
 729         }
 730
 731       } else {
 732
 733         pair += "2fakes";
 734
 735         if ((!fMatchTrig || (trigger1 && trigger2)) && (!fApplyAccCut || inAcc)) {
 736           ((TH1F*)fList2->UncheckedAt(k2MassF2))->Fill(mass);
 737           ((TH1F*)fList2->UncheckedAt(k2PF2))->Fill(p);
 738           ((TH1F*)fList2->UncheckedAt(k2PtF2))->Fill(pt);
 739           ((TH1F*)fList2->UncheckedAt(k2YF2))->Fill(y);
 740           ((TH1F*)fList2->UncheckedAt(k2EtaF2))->Fill(eta);
 741           ((TH1F*)fList2->UncheckedAt(k2PhiF2))->Fill(phi);
 742         }
 743
 744       }
 745
 746       // fill global counters
 747       fPairCounters->Count(Form("pair:reconstructed/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 748       fPairCounters->Count(Form("%s/run:%d/%s/%s/%s/%s", pair.Data(), fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 749
 750     }
 751
 752   }
 753
 754   // Post final data
 755   PostData(1, fList);
 756   PostData(2, fTrackCounters);
 757   PostData(3, fFakeTrackCounters);
 758   PostData(4, fMatchedTrackCounters);
 759   PostData(5, fEventCounters);
 760   PostData(6, fList2);
 761   PostData(7, fPairCounters);
 762 }
 763
 764 //________________________________________________________________________
 765 void AliAnalysisTaskMuonFakes::NotifyRun()
 766 {
 767   /// Prepare processing of new run: load corresponding OCDB objects...
 768
 769   // load necessary data from OCDB
 770   AliCDBManager::Instance()->SetDefaultStorage(fRecoParamLocation.Data());
 771   AliCDBManager::Instance()->SetRun(fCurrentRunNumber);
 772   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
 773   if (!recoParam) {
 774     fRequestedStationMask = 0;
 775     fRequest2ChInSameSt45 = kFALSE;
 776     fSigmaCut = -1.;
 777     AliError("--> skip this run");
 778     return;
 779   }
 780
 781   // compute the mask of requested stations from recoParam
 782   fRequestedStationMask = 0;
 783   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) fRequestedStationMask |= ( 1 << i );
 784
 785   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
 786   fRequest2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
 787
 788   // get sigma cut from recoParam to associate clusters with TrackRefs in case the labels are not used
 789   fSigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
 790 }
 791
 792 //________________________________________________________________________
 793 void AliAnalysisTaskMuonFakes::Terminate(Option_t *)
 794 {
 795   /// Draw results to the screen and print statistics.
 796
 797   // recover output objects
 798   fList = static_cast<TObjArray*> (GetOutputData(1));
 799   fList2 = static_cast<TObjArray*> (GetOutputData(6));
 800   if (!fList || !fList2) return;
 801   fTrackCounters = static_cast<AliCounterCollection*> (GetOutputData(2));
 802   fFakeTrackCounters = static_cast<AliCounterCollection*> (GetOutputData(3));
 803   fMatchedTrackCounters = static_cast<AliCounterCollection*> (GetOutputData(4));
 804   fEventCounters = static_cast<AliCounterCollection*> (GetOutputData(5));
 805   fPairCounters = static_cast<AliCounterCollection*> (GetOutputData(7));
 806
 807   // add canvas to compare histograms
 808   fCanvases = new TObjArray(3);
 809   fCanvases->SetOwner();
 810   TCanvas *cFakesSummary1 = new TCanvas("cFakesSummary1","cFakesSummary1",900,900);
 811   fCanvases->AddAtAndExpand(cFakesSummary1, 0);
 812   TCanvas *cFakesSummary2 = new TCanvas("cFakesSummary2","cFakesSummary2",1200,600);
 813   fCanvases->AddAtAndExpand(cFakesSummary2, 1);
 814   TCanvas *cFakesSummary3 = new TCanvas("cFakesSummary3","cFakesSummary3",900,600);
 815   fCanvases->AddAtAndExpand(cFakesSummary3, 2);
 816
 817   // display
 818   Int_t iHist1[9] = {kNumberOfClusters, kNumberOfChamberHit, kChi2PerDof, kDCA, kRAbs, kEta, kP, kPt, kPhi};
 819   cFakesSummary1->Divide(3,3);
 820   for (Int_t i=0; i<9; i++) {
 821     cFakesSummary1->cd(i+1);
 822     cFakesSummary1->GetPad(i+1)->SetLogy();
 823     ((TH1F*)fList->UncheckedAt(iHist1[i]))->SetMinimum(0.5);
 824     ((TH1F*)fList->UncheckedAt(iHist1[i]))->DrawCopy();
 825     ((TH1F*)fList->UncheckedAt(iHist1[i]+1))->SetLineColor(4);
 826     ((TH1F*)fList->UncheckedAt(iHist1[i]+1))->DrawCopy("sames");
 827     ((TH1F*)fList->UncheckedAt(iHist1[i]+2))->SetLineColor(2);
 828     ((TH1F*)fList->UncheckedAt(iHist1[i]+2))->SetFillColor(2);
 829     ((TH1F*)fList->UncheckedAt(iHist1[i]+2))->SetFillStyle(3017);
 830     ((TH1F*)fList->UncheckedAt(iHist1[i]+2))->DrawCopy("sames");
 831   }
 832
 833   Int_t iHist2[2] = {kChi2PerDofVsNClusters, kChi2PerDofVsNChamberHit};
 834   cFakesSummary2->Divide(2);
 835   for (Int_t i=0; i<2; i++) {
 836     cFakesSummary2->cd(i+1);
 837     ((TH2F*)fList->UncheckedAt(iHist2[i]+1))->SetMarkerColor(4);
 838     ((TH2F*)fList->UncheckedAt(iHist2[i]+1))->DrawCopy();
 839     ((TH2F*)fList->UncheckedAt(iHist2[i]+2))->SetMarkerColor(2);
 840     ((TH2F*)fList->UncheckedAt(iHist2[i]+2))->SetMarkerStyle(7);
 841     ((TH2F*)fList->UncheckedAt(iHist2[i]+2))->DrawCopy("sames");
 842   }
 843
 844   Int_t iHist3[6] = {k2Mass, k2P, k2Pt, k2Y, k2Eta, k2Phi};
 845   cFakesSummary3->Divide(3,2);
 846   for (Int_t i=0; i<6; i++) {
 847     cFakesSummary3->cd(i+1);
 848     cFakesSummary3->GetPad(i+1)->SetLogy();
 849     ((TH1F*)fList2->UncheckedAt(iHist3[i]))->SetMinimum(0.5);
 850     ((TH1F*)fList2->UncheckedAt(iHist3[i]))->DrawCopy();
 851     ((TH1F*)fList2->UncheckedAt(iHist3[i]+1))->SetLineColor(4);
 852     ((TH1F*)fList2->UncheckedAt(iHist3[i]+1))->DrawCopy("sames");
 853     ((TH1F*)fList2->UncheckedAt(iHist3[i]+2))->Add((TH1F*)fList2->UncheckedAt(iHist3[i]+3));
 854     ((TH1F*)fList2->UncheckedAt(iHist3[i]+2))->SetLineColor(2);
 855     ((TH1F*)fList2->UncheckedAt(iHist3[i]+2))->SetFillColor(2);
 856     ((TH1F*)fList2->UncheckedAt(iHist3[i]+2))->SetFillStyle(3017);
 857     ((TH1F*)fList2->UncheckedAt(iHist3[i]+2))->DrawCopy("sames");
 858     ((TH1F*)fList2->UncheckedAt(iHist3[i]+3))->SetLineColor(6);
 859     ((TH1F*)fList2->UncheckedAt(iHist3[i]+3))->SetFillColor(6);
 860     ((TH1F*)fList2->UncheckedAt(iHist3[i]+3))->SetFillStyle(3018);
 861     ((TH1F*)fList2->UncheckedAt(iHist3[i]+3))->DrawCopy("sames");
 862   }
 863
 864   // print
 865   if (fTrackCounters && fFakeTrackCounters && fMatchedTrackCounters && fEventCounters) {
 866     printf("\nGlobal statistics of reconstructed tracks matched or not with the trigger:\n");
 867     fTrackCounters->Print("track/trig");
 868     printf("\nGlobal statistics of pathological tracks matched or not with the trigger:\n");
 869     fFakeTrackCounters->Print("track/trig");
 870     printf("\nDetailled statistics of tracks matched per label vs position:\n");
 871     fMatchedTrackCounters->Print("label/position");
 872     printf("\nGlobal statistics of events containing pathological tracks:\n");
 873     fEventCounters->Print("event/trig");
 874   }
 875
 876   if (fPairCounters) {
 877     printf("\nGlobal statistics of reconstructed track pairs matched or not with the trigger:\n");
 878     fPairCounters->Print("pair/trig");
 879   }
 880
 881   printf("\nREMINDER: results are relevent provided that you use the same recoParams as for the reconstruction\n\n");
 882 }
 883
 884 //________________________________________________________________________
 885 Int_t AliAnalysisTaskMuonFakes::RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
 886                                                      TString &selected, TString &centrality)
 887 {
 888   /// loop over reconstructible TrackRef not associated with reconstructed track:
 889   /// for each of them, find and remove the most connected the fake track, if any,
 890   /// and fill the histograms with the fraction of connected clusters.
 891   /// Return the number of reconstructible track not connected to any fake
 892
 893   Int_t nFreeMissingTracks = 0;
 894
 895   // loop over trackRefs
 896   TIter next(trackRefStore.CreateIterator());
 897   AliMUONTrack* trackRef;
 898   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
 899
 900     // skip not reconstructible trackRefs
 901     if (!trackRef->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) continue;
 902
 903     Int_t label = trackRef->GetUniqueID();
 904
 905     // look for the most connected fake track
 906     AliMUONTrack *connectedFake = 0x0;
 907     Double_t fractionOfConnectedClusters = 0.;
 908     TIter next2(fakeTrackStore.CreateIterator());
 909     AliMUONTrack* fakeTrack;
 910     while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
 911
 912       // get the number of connected clusters
 913       Int_t nConnectedClusters = 0;
 914       if (fUseLabel) { // by using the MC label
 915         for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
 916           if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
 917             nConnectedClusters++;
 918       } else { // by comparing cluster/TrackRef positions
 919         Bool_t compTrack[10];
 920         nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, fSigmaCut, compTrack);
 921       }
 922
 923       // skip non-connected fake tracks
 924       if (nConnectedClusters == 0) continue;
 925
 926       // check if it is the most connected fake track
 927       Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters());
 928       if (f > fractionOfConnectedClusters) {
 929         connectedFake = fakeTrack;
 930         fractionOfConnectedClusters = f;
 931       }
 932
 933     }
 934
 935     if (connectedFake) {
 936
 937       // find the corresponding ESD MUON track
 938       AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 939       if (!esd) {
 940         AliError("Cannot get input event");
 941         return nFreeMissingTracks;
 942       }
 943       TIter next3(static_cast<TClonesArray*>(esd->FindListObject("MuonTracks")));
 944       AliESDMuonTrack* esdTrack = 0x0;
 945       while ((esdTrack = static_cast<AliESDMuonTrack*>(next3())) && esdTrack->GetUniqueID() != connectedFake->GetUniqueID()) {}
 946       if (!esdTrack) {
 947         AliError("unable to find the corresponding ESD track???");
 948         continue;
 949       }
 950
 951       // trigger condition
 952       TString trig = (esdTrack->ContainTriggerData()) ? "trig:yes" : "trig:no";
 953
 954       // acceptance condition
 955       Double_t thetaTrackAbsEnd = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 956       Double_t eta = esdTrack->Eta();
 957       TString acc = (thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 10. && eta >= -4. && eta <= -2.5) ? "acc:in" : "acc:out";
 958
 959       // fill histogram and counters
 960       ((TH1F*)fList->UncheckedAt(kFractionOfConnectedClusters))->Fill(fractionOfConnectedClusters);
 961       fTrackCounters->Count(Form("track:connected/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 962       fFakeTrackCounters->Count(Form("track:connected/run:%d/file:%s/event:%d/%s/%s/%s/%s", fCurrentRunNumber, fCurrentFileName.Data(),
 963                                      esd->GetEventNumberInFile(), trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
 964
 965       // remove the most connected fake track
 966       fakeTrackStore.Remove(*connectedFake);
 967
 968     } else nFreeMissingTracks++;
 969
 970   }
 971
 972   return nFreeMissingTracks;
 973
 974 }
 975