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 // ROOT includes
  17 #include <TObjArray.h>
  18 #include <TClonesArray.h>
  19 #include <TH1F.h>
  20 #include <TH2F.h>
  21 #include <TCanvas.h>
  22
  23 // STEER includes
  24 #include "AliLog.h"
  25 #include "AliESDEvent.h"
  26 #include "AliESDMuonTrack.h"
  27 #include "AliESDInputHandler.h"
  28 #include "AliMCEventHandler.h"
  29
  30 // ANALYSIS includes
  31 #include "AliAnalysisDataSlot.h"
  32 #include "AliAnalysisDataContainer.h"
  33 #include "AliAnalysisManager.h"
  34 #include "AliCDBManager.h"
  35
  36 // MUON includes
  37 #include "AliMUONCDB.h"
  38 #include "AliMUONRecoParam.h"
  39 #include "AliMUONRecoCheck.h"
  40 #include "AliMUONVCluster.h"
  41 #include "AliMUONVTrackStore.h"
  42 #include "AliMUONVTriggerTrackStore.h"
  43 #include "AliMUONTrack.h"
  44 #include "AliMUONTrackParam.h"
  45 #include "AliMUONESDInterface.h"
  46
  47 #include "AliAnalysisTaskMuonFakes.h"
  48 #include "AliCounterCollection.h"
  49
  50 ClassImp(AliAnalysisTaskMuonFakes)
  51
  52 //________________________________________________________________________
  53 AliAnalysisTaskMuonFakes::AliAnalysisTaskMuonFakes() :
  54 AliAnalysisTaskSE(),
  55 fList(0x0),
  56 fCanvases(0x0),
  57 fTrackCounters(0x0),
  58 fFakeTrackCounters(0x0),
  59 fMatchedTrackCounters(0x0),
  60 fEventCounters(0x0),
  61 fCurrentFileName(""),
  62 fRequestedStationMask(0),
  63 fRequest2ChInSameSt45(kFALSE),
  64 fSigmaCut(-1.),
  65 fUseLabel(kFALSE),
  66 fRecoParamLocation("alien://folder=/alice/simulation/2008/v4-15-Release/Full")
  67 {
  68   /// Default constructor.
  69 }
  70
  71 //________________________________________________________________________
  72 AliAnalysisTaskMuonFakes::AliAnalysisTaskMuonFakes(const char *name) :
  73 AliAnalysisTaskSE(name),
  74 fList(0x0),
  75 fCanvases(0x0),
  76 fTrackCounters(0x0),
  77 fFakeTrackCounters(0x0),
  78 fMatchedTrackCounters(0x0),
  79 fEventCounters(0x0),
  80 fCurrentFileName(""),
  81 fRequestedStationMask(0),
  82 fRequest2ChInSameSt45(kFALSE),
  83 fSigmaCut(-1.),
  84 fUseLabel(kFALSE),
  85 fRecoParamLocation("alien://folder=/alice/simulation/2008/v4-15-Release/Full")
  86 {
  87   /// Constructor.
  88   // Output slot #1 writes into a TObjArray container
  89   DefineOutput(1,TObjArray::Class());
  90   // Output slot #2 writes into an AliCounterCollection container
  91   DefineOutput(2,AliCounterCollection::Class());
  92   // Output slot #3 writes into an AliCounterCollection container
  93   DefineOutput(3,AliCounterCollection::Class());
  94   // Output slot #4 writes into an AliCounterCollection container
  95   DefineOutput(4,AliCounterCollection::Class());
  96   // Output slot #5 writes into an AliCounterCollection container
  97   DefineOutput(5,AliCounterCollection::Class());
  98 }
  99
 100 //________________________________________________________________________
 101 AliAnalysisTaskMuonFakes::~AliAnalysisTaskMuonFakes()
 102 {
 103   /// Destructor.
 104   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 105     delete fList;
 106     delete fTrackCounters;
 107     delete fFakeTrackCounters;
 108     delete fMatchedTrackCounters;
 109     delete fEventCounters;
 110   }
 111   delete fCanvases;
 112 }
 113
 114 //___________________________________________________________________________
 115 void AliAnalysisTaskMuonFakes::UserCreateOutputObjects()
 116 {
 117   /// Create histograms and counters.
 118
 119   fList = new TObjArray(100);
 120   fList->SetOwner();
 121
 122   // number of tracks
 123   TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks", "nb of tracks /evt", 21, -0.5, 20.5);
 124   fList->AddAtAndExpand(hNumberOfTracks, kNumberOfTracks);
 125   TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks", "nb of fake - nb of missing track", 21, -0.5, 20.5);
 126   fList->AddAtAndExpand(hNumberOfAdditionalTracks, kNumberOfAdditionalTracks);
 127
 128   // number of clusters
 129   TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters", "nb of clusters /track", 21, -0.5, 20.5);
 130   fList->AddAtAndExpand(hNumberOfClusters, kNumberOfClusters);
 131   TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM", "nb of clusters /matched track", 21, -0.5, 20.5);
 132   fList->AddAtAndExpand(hNumberOfClustersM, kNumberOfClustersM);
 133   TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF", "nb of clusters /fake track", 21, -0.5, 20.5);
 134   fList->AddAtAndExpand(hNumberOfClustersF, kNumberOfClustersF);
 135   TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC", "nb of clusters /MC track", 21, -0.5, 20.5);
 136   fList->AddAtAndExpand(hNumberOfClustersMC, kNumberOfClustersMC);
 137   TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters", "nb of matched clusters / nb of clusters", 110, 0., 1.1);
 138   fList->AddAtAndExpand(hFractionOfMatchedClusters, kFractionOfMatchedClusters);
 139   TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters", "nb of connected clusters / nb of clusters in fake tracks", 110, 0., 1.1);
 140   fList->AddAtAndExpand(hFractionOfConnectedClusters, kFractionOfConnectedClusters);
 141
 142   // number of fired chambers
 143   TH1F *hNumberOfChamberHit = new TH1F("hNumberOfChamberHit", "nb of chambers hit /track", 16, -0.5, 15.5);
 144   fList->AddAtAndExpand(hNumberOfChamberHit, kNumberOfChamberHit);
 145   TH1F *hNumberOfChamberHitM = new TH1F("hNumberOfChamberHitM", "nb of chambers hit /matched track", 16, -0.5, 15.5);
 146   fList->AddAtAndExpand(hNumberOfChamberHitM, kNumberOfChamberHitM);
 147   TH1F *hNumberOfChamberHitF = new TH1F("hNumberOfChamberHitF", "nb of chambers hit /fake track", 16, -0.5, 15.5);
 148   fList->AddAtAndExpand(hNumberOfChamberHitF, kNumberOfChamberHitF);
 149
 150   // chi2
 151   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 100, 0., 20.);
 152   fList->AddAtAndExpand(hChi2PerDof, kChi2PerDof);
 153   TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 100, 0., 20.);
 154   fList->AddAtAndExpand(hChi2PerDofM, kChi2PerDofM);
 155   TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 100, 0., 20.);
 156   fList->AddAtAndExpand(hChi2PerDofF, kChi2PerDofF);
 157
 158   // chi2 versus number of clusters
 159   TH2F *hChi2PerDofVsNClusters = new TH2F("hChi2PerDofVsNClusters", "track chi2/d.o.f. versus nb of clusters", 21, -0.5, 20.5, 100, 0., 20.);
 160   fList->AddAtAndExpand(hChi2PerDofVsNClusters, kChi2PerDofVsNClusters);
 161   TH2F *hChi2PerDofVsNClustersM = new TH2F("hChi2PerDofVsNClustersM", "matched track chi2/d.o.f. versus nb of clusters", 21, -0.5, 20.5, 100, 0., 20.);
 162   fList->AddAtAndExpand(hChi2PerDofVsNClustersM, kChi2PerDofVsNClustersM);
 163   TH2F *hChi2PerDofVsNClustersF = new TH2F("hChi2PerDofVsNClustersF", "fake track chi2/d.o.f. versus nb of clusters", 21, -0.5, 20.5, 100, 0., 20.);
 164   fList->AddAtAndExpand(hChi2PerDofVsNClustersF, kChi2PerDofVsNClustersF);
 165
 166   // chi2 versus number of fired chambers
 167   TH2F *hChi2PerDofVsNChamberHit = new TH2F("hChi2PerDofVsNChamberHit", "track chi2/d.o.f. versus nb of fired chambers", 16, -0.5, 15.5, 100, 0., 20.);
 168   fList->AddAtAndExpand(hChi2PerDofVsNChamberHit, kChi2PerDofVsNChamberHit);
 169   TH2F *hChi2PerDofVsNChamberHitM = new TH2F("hChi2PerDofVsNChamberHitM", "matched track chi2/d.o.f. versus nb of fired chambers", 16, -0.5, 15.5, 100, 0., 20.);
 170   fList->AddAtAndExpand(hChi2PerDofVsNChamberHitM, kChi2PerDofVsNChamberHitM);
 171   TH2F *hChi2PerDofVsNChamberHitF = new TH2F("hChi2PerDofVsNChamberHitF", "fake track chi2/d.o.f. versus nb of fired chambers", 16, -0.5, 15.5, 100, 0., 20.);
 172   fList->AddAtAndExpand(hChi2PerDofVsNChamberHitF, kChi2PerDofVsNChamberHitF);
 173
 174   // physics quantities
 175   TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
 176   fList->AddAtAndExpand(hP, kP);
 177   TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
 178   fList->AddAtAndExpand(hPM, kPM);
 179   TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
 180   fList->AddAtAndExpand(hPF, kPF);
 181   TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
 182   fList->AddAtAndExpand(hPt, kPt);
 183   TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
 184   fList->AddAtAndExpand(hPtM, kPtM);
 185   TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
 186   fList->AddAtAndExpand(hPtF, kPtF);
 187   TH1F *hEta = new TH1F("hEta", "Muon pseudo-rapidity distribution", 100, -10., 0.);
 188   fList->AddAtAndExpand(hEta , kEta );
 189   TH1F *hEtaM = new TH1F("hEtaM", "matched track pseudo-rapidity distribution", 100, -10., 0.);
 190   fList->AddAtAndExpand(hEtaM, kEtaM);
 191   TH1F *hEtaF = new TH1F("hEtaF", "fake track pseudo-rapidity distribution", 100, -10., 0.);
 192   fList->AddAtAndExpand(hEtaF, kEtaF);
 193   TH1F *hPhi = new TH1F("hPhi", "Muon phi distribution", 100, -1., 9.);
 194   fList->AddAtAndExpand(hPhi, kPhi);
 195   TH1F *hPhiM = new TH1F("hPhiM", "matched track phi distribution", 100, -1., 9.);
 196   fList->AddAtAndExpand(hPhiM, kPhiM);
 197   TH1F *hPhiF = new TH1F("hPhiF", "fake track phi distribution", 100, -1., 9.);
 198   fList->AddAtAndExpand(hPhiF, kPhiF);
 199   TH1F *hDCA = new TH1F("hDCA", "Muon DCA distribution", 125, 0., 500.);
 200   fList->AddAtAndExpand(hDCA, kDCA);
 201   TH1F *hDCAM = new TH1F("hDCAM", "matched track DCA distribution", 125, 0., 500.);
 202   fList->AddAtAndExpand(hDCAM, kDCAM);
 203   TH1F *hDCAF = new TH1F("hDCAF", "fake track DCA distribution", 125, 0., 500.);
 204   fList->AddAtAndExpand(hDCAF, kDCAF);
 205
 206   // global counters of tracks:
 207   // - reconstructible = number of reconstructible tracks
 208   // - reconstructed   = number of reconstructed tracks
 209   // - matched         = number of reconstructed tracks matched with a simulated one (reconstructible or not)
 210   // - matchedyet      = number of reconstructed tracks matched with a simulated one that is not reconstructible
 211   // - fake            = number of fake tracks
 212   // - connected       = number of fake tracks connected to a reconstructible simulated track
 213   // - additional      = number of additional (fake) tracks compared to the number of reconstructible ones
 214   fTrackCounters = new AliCounterCollection(GetOutputSlot(2)->GetContainer()->GetName());
 215   fTrackCounters->AddRubric("track", "reconstructible/reconstructed/matched/matchedyet/fake/connected/additional");
 216   fTrackCounters->AddRubric("run", 1000000);
 217   fTrackCounters->AddRubric("trig", "yes/no/unknown");
 218   fTrackCounters->AddRubric("selected", "yes/no");
 219   fTrackCounters->AddRubric("acc", "in/out/unknown");
 220   fTrackCounters->Init();
 221
 222   // detailled counters of fake tracks:
 223   fFakeTrackCounters = new AliCounterCollection(GetOutputSlot(3)->GetContainer()->GetName());
 224   fFakeTrackCounters->AddRubric("track", "fake/connected/additional/matchedyet/fake?");
 225   fFakeTrackCounters->AddRubric("run", 1000000);
 226   fFakeTrackCounters->AddRubric("file", 1000000);
 227   fFakeTrackCounters->AddRubric("event", 1000000);
 228   fFakeTrackCounters->AddRubric("trig", "yes/no/unknown");
 229   fFakeTrackCounters->AddRubric("selected", "yes/no");
 230   fFakeTrackCounters->AddRubric("acc", "in/out/unknown");
 231   fFakeTrackCounters->Init();
 232
 233   // counters of tracks matched by position or by using MC labels
 234   fMatchedTrackCounters = new AliCounterCollection(GetOutputSlot(4)->GetContainer()->GetName());
 235   fMatchedTrackCounters->AddRubric("position", "match/not match");
 236   fMatchedTrackCounters->AddRubric("label", "match/not match/match other");
 237   fMatchedTrackCounters->AddRubric("run", 1000000);
 238   fMatchedTrackCounters->AddRubric("trig", "yes/no");
 239   fMatchedTrackCounters->AddRubric("selected", "yes/no");
 240   fMatchedTrackCounters->AddRubric("acc", "in/out");
 241   fMatchedTrackCounters->Init();
 242
 243   // global counters of events
 244   // - any             = total number of events with reconstructed tracks
 245   // - fake            = number of events with fake track(s)
 246   // - notconnected    = number of events with fake tracks that are not connected to a reconstructible simulated track
 247   // - additional      = number of events with additional (fake) tracks compared to the number of reconstructible ones
 248   // - matchedyet      = number of events with reconstructed tracks matched with a simulated one that is not reconstructible
 249   // if trig = yes: only the tracks matched with the trigger are considered in the above logic
 250   fEventCounters = new AliCounterCollection(GetOutputSlot(5)->GetContainer()->GetName());
 251   fEventCounters->AddRubric("event", "any/fake/notconnected/additional/matchedyet");
 252   fEventCounters->AddRubric("run", 1000000);
 253   fEventCounters->AddRubric("trig", "any/yes");
 254   fEventCounters->AddRubric("selected", "yes/no");
 255   fEventCounters->Init();
 256
 257   // Disable printout of AliMCEvent
 258   AliLog::SetClassDebugLevel("AliMCEvent",-1);
 259
 260   // Post data at least once per task to ensure data synchronisation (required for merging)
 261   PostData(1, fList);
 262   PostData(2, fTrackCounters);
 263   PostData(3, fFakeTrackCounters);
 264   PostData(4, fMatchedTrackCounters);
 265   PostData(5, fEventCounters);
 266 }
 267
 268 //________________________________________________________________________
 269 void AliAnalysisTaskMuonFakes::UserExec(Option_t *)
 270 {
 271   /// Process event: looks for fakes...
 272
 273   // check that reconstructions parameters for that run have been properly set
 274   if (fSigmaCut < 0) return;
 275
 276   // check physics selection
 277   TString selected = (fInputHandler && fInputHandler->IsEventSelected() != 0) ? "selected:yes" : "selected:no";
 278
 279   // current file name
 280   fCurrentFileName = CurrentFileName();
 281   fCurrentFileName.ReplaceAll("alien://","");
 282   fCurrentFileName.ReplaceAll("/","\\");
 283   fCurrentFileName.ReplaceAll(":",";");
 284
 285   // Load ESD event
 286   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 287
 288   // Load MC event
 289   AliMCEventHandler *mcH = 0;
 290   if(MCEvent()) mcH = static_cast<AliMCEventHandler*>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 291
 292   // get reconstructed and simulated tracks
 293   AliMUONRecoCheck rc(esd,mcH);
 294   AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(-1, kFALSE);
 295   AliMUONVTrackStore* trackRefStore = rc.TrackRefs(-1);
 296   AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(-1);
 297   if (!muonTrackStore || !trackRefStore) return;
 298
 299   // count the number of reconstructible tracks
 300   TIter next(trackRefStore->CreateIterator());
 301   AliMUONTrack* trackRef;
 302   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
 303     if (trackRef->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) {
 304       TString trig = (triggerTrackRefStore->FindObject(trackRef->GetUniqueID())) ? "trig:yes" : "trig:no";
 305       fTrackCounters->Count(Form("track:reconstructible/run:%d/%s/%s/acc:unknown", fCurrentRunNumber, trig.Data(), selected.Data()));
 306     }
 307   }
 308
 309   // loop over ESD tracks
 310   Int_t nTrackerTracks = 0;
 311   Bool_t containTrack[2] = {kFALSE, kFALSE};
 312   Bool_t containFakeTrack[2] = {kFALSE, kFALSE};
 313   Bool_t containMatchedYetTrack[2] = {kFALSE, kFALSE};
 314   AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
 315   Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
 316   for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 317
 318     AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
 319
 320     // skip ghosts
 321     if (!esdTrack->ContainTrackerData()) continue;
 322     nTrackerTracks++;
 323     containTrack[0] = kTRUE;
 324
 325     // trigger condition
 326     Bool_t trigger = esdTrack->ContainTriggerData();
 327     TString trig = "trig:";
 328     if (trigger) {
 329       trig += "yes";
 330       containTrack[1] = kTRUE;
 331     } else trig += "no";
 332
 333     // acceptance condition
 334     TString acc = "acc:";
 335     Double_t thetaTrackAbsEnd = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 336     Double_t eta = esdTrack->Eta();
 337     if (thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 9. && eta >= -4. && eta <= -2.5) acc += "in";
 338     else acc += "out";
 339
 340     // fill global counters
 341     fTrackCounters->Count(Form("track:reconstructed/run:%d/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data()));
 342
 343     // find the corresponding MUON track
 344     AliMUONTrack* muonTrack = static_cast<AliMUONTrack*>(muonTrackStore->FindObject(esdTrack->GetUniqueID()));
 345
 346     // get track info
 347     Int_t nClusters = esdTrack->GetNClusters();
 348     Int_t nChamberHit = 0;
 349     for (Int_t ich=0; ich<10; ich++) if (esdTrack->IsInMuonClusterMap(ich)) nChamberHit++;
 350     Double_t normalizedChi2 = esdTrack->GetChi2() / (2. * esdTrack->GetNHit() - 5);
 351     Double_t p = esdTrack->P();
 352     Double_t pT = esdTrack->Pt();
 353     Double_t phi = esdTrack->Phi();
 354     Double_t dca = esdTrack->GetDCA();
 355
 356     // fill global histograms
 357     ((TH1F*)fList->UncheckedAt(kNumberOfClusters))->Fill(nClusters);
 358     ((TH1F*)fList->UncheckedAt(kNumberOfChamberHit))->Fill(nChamberHit);
 359     ((TH1F*)fList->UncheckedAt(kChi2PerDof))->Fill(normalizedChi2);
 360     ((TH1F*)fList->UncheckedAt(kP))->Fill(p);
 361     ((TH1F*)fList->UncheckedAt(kPt))->Fill(pT);
 362     ((TH1F*)fList->UncheckedAt(kEta))->Fill(eta);
 363     ((TH1F*)fList->UncheckedAt(kPhi))->Fill(phi);
 364     ((TH1F*)fList->UncheckedAt(kDCA))->Fill(dca);
 365     ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNClusters))->Fill(nClusters,normalizedChi2);
 366     ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNChamberHit))->Fill(nChamberHit,normalizedChi2);
 367
 368     // try to match, by position, the reconstructed track with a simulated one
 369     Int_t nMatchClustersByPosition = 0;
 370     AliMUONTrack* matchedTrackRefByPosition = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClustersByPosition, kFALSE, fSigmaCut);
 371     Int_t MCLabelByPosition = (matchedTrackRefByPosition) ? matchedTrackRefByPosition->GetUniqueID() : -1;
 372
 373     // try to match, by using MC labels, the reconstructed track with a simulated one
 374     Int_t nMatchClustersByLabel = 0;
 375     AliMUONTrack* matchedTrackRefByLabel = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClustersByLabel, kTRUE, fSigmaCut);
 376     Int_t MCLabelByLabel = (matchedTrackRefByLabel) ? matchedTrackRefByLabel->GetUniqueID() : -1;
 377
 378     // fill global counters
 379     TString positionCase = (MCLabelByPosition >= 0) ? "position:match" : "position:not match";
 380     TString labelCase = "label:";
 381     if (MCLabelByLabel >= 0 && MCLabelByPosition >= 0 && MCLabelByLabel != MCLabelByPosition) labelCase += "match other";
 382     else if (MCLabelByLabel >= 0) labelCase += "match";
 383     else labelCase += "not match";
 384     fMatchedTrackCounters->Count(Form("%s/%s/run:%d/%s/%s/%s", positionCase.Data(), labelCase.Data(), fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data()));
 385     if ((MCLabelByLabel >= 0 && MCLabelByPosition < 0) || (MCLabelByLabel < 0 && MCLabelByPosition >= 0))
 386       fFakeTrackCounters->Count(Form("track:fake?/run:%d/file:%s/event:%d/%s/%s/%s", fCurrentRunNumber, fCurrentFileName.Data(),
 387                                      esd->GetEventNumberInFile(), trig.Data(), selected.Data(), acc.Data()));
 388
 389     // take actions according to the matching result we are interested in
 390     Int_t nMatchClusters = (fUseLabel) ? nMatchClustersByLabel : nMatchClustersByPosition;
 391     AliMUONTrack* matchedTrackRef = (fUseLabel) ? matchedTrackRefByLabel : matchedTrackRefByPosition;
 392     if (matchedTrackRef) {
 393
 394       // fill global counters
 395       fTrackCounters->Count(Form("track:matched/run:%d/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data()));
 396
 397       // track matched with a trackRef that is not reconstructible
 398       if (!matchedTrackRef->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) {
 399
 400         containMatchedYetTrack[0] = kTRUE;
 401         if (trigger) containMatchedYetTrack[1] = kTRUE;
 402
 403         // fill global counters
 404         fTrackCounters->Count(Form("track:matchedyet/run:%d/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data()));
 405         fFakeTrackCounters->Count(Form("track:matchedyet/run:%d/file:%s/event:%d/%s/%s/%s", fCurrentRunNumber, fCurrentFileName.Data(),
 406                                        esd->GetEventNumberInFile(), trig.Data(), selected.Data(), acc.Data()));
 407       }
 408
 409       // fill histograms
 410       ((TH1F*)fList->UncheckedAt(kFractionOfMatchedClusters))->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
 411       ((TH1F*)fList->UncheckedAt(kNumberOfClustersMC))->Fill(matchedTrackRef->GetNClusters());
 412       ((TH1F*)fList->UncheckedAt(kNumberOfClustersM))->Fill(nClusters);
 413       ((TH1F*)fList->UncheckedAt(kNumberOfChamberHitM))->Fill(nChamberHit);
 414       ((TH1F*)fList->UncheckedAt(kChi2PerDofM))->Fill(normalizedChi2);
 415       ((TH1F*)fList->UncheckedAt(kPM))->Fill(p);
 416       ((TH1F*)fList->UncheckedAt(kPtM))->Fill(pT);
 417       ((TH1F*)fList->UncheckedAt(kEtaM))->Fill(eta);
 418       ((TH1F*)fList->UncheckedAt(kPhiM))->Fill(phi);
 419       ((TH1F*)fList->UncheckedAt(kDCAM))->Fill(dca);
 420       ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNClustersM))->Fill(nClusters,normalizedChi2);
 421       ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNChamberHitM))->Fill(nChamberHit,normalizedChi2);
 422
 423       // remove already matched trackRefs
 424       trackRefStore->Remove(*matchedTrackRef);
 425
 426     } else {
 427
 428       containFakeTrack[0] = kTRUE;
 429       if (trigger) containFakeTrack[1] = kTRUE;
 430
 431       // fill global counters
 432       fTrackCounters->Count(Form("track:fake/run:%d/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data()));
 433       fFakeTrackCounters->Count(Form("track:fake/run:%d/file:%s/event:%d/%s/%s/%s", fCurrentRunNumber, fCurrentFileName.Data(),
 434                                      esd->GetEventNumberInFile(), trig.Data(), selected.Data(), acc.Data()));
 435
 436       // fill histograms
 437       ((TH1F*)fList->UncheckedAt(kNumberOfClustersF))->Fill(nClusters);
 438       ((TH1F*)fList->UncheckedAt(kNumberOfChamberHitF))->Fill(nChamberHit);
 439       ((TH1F*)fList->UncheckedAt(kChi2PerDofF))->Fill(normalizedChi2);
 440       ((TH1F*)fList->UncheckedAt(kPF))->Fill(p);
 441       ((TH1F*)fList->UncheckedAt(kPtF))->Fill(pT);
 442       ((TH1F*)fList->UncheckedAt(kEtaF))->Fill(eta);
 443       ((TH1F*)fList->UncheckedAt(kPhiF))->Fill(phi);
 444       ((TH1F*)fList->UncheckedAt(kDCAF))->Fill(dca);
 445       ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNClustersF))->Fill(nClusters,normalizedChi2);
 446       ((TH1F*)fList->UncheckedAt(kChi2PerDofVsNChamberHitF))->Fill(nChamberHit,normalizedChi2);
 447
 448       // store fake tracks
 449       fakeTrackStore->Add(*muonTrack);
 450
 451     }
 452
 453   } // end of loop over ESD tracks
 454
 455   // fill histogram and global counters
 456   ((TH1F*)fList->UncheckedAt(kNumberOfTracks))->Fill(nTrackerTracks);
 457   if (containTrack[0]) fEventCounters->Count(Form("event:any/run:%d/trig:any/%s", fCurrentRunNumber, selected.Data()));
 458   if (containTrack[1]) fEventCounters->Count(Form("event:any/run:%d/trig:yes/%s", fCurrentRunNumber, selected.Data()));
 459   if (containFakeTrack[0]) fEventCounters->Count(Form("event:fake/run:%d/trig:any/%s", fCurrentRunNumber, selected.Data()));
 460   if (containFakeTrack[1]) fEventCounters->Count(Form("event:fake/run:%d/trig:yes/%s", fCurrentRunNumber, selected.Data()));
 461   if (containMatchedYetTrack[0]) fEventCounters->Count(Form("event:matchedyet/run:%d/trig:any/%s", fCurrentRunNumber, selected.Data()));
 462   if (containMatchedYetTrack[1]) fEventCounters->Count(Form("event:matchedyet/run:%d/trig:yes/%s", fCurrentRunNumber, selected.Data()));
 463
 464   // count the number of not connected and additional fake tracks
 465   if (fakeTrackStore->GetSize() > 0) {
 466
 467     // remove the most connected fake tracks
 468     Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore);
 469
 470     if (fakeTrackStore->GetSize() > 0) {
 471
 472       // fill global counters
 473       fEventCounters->Count(Form("event:notconnected/run:%d/trig:any/%s", fCurrentRunNumber, selected.Data()));
 474
 475       // check status of remaining fakes with respect to the matching with trigger
 476       Bool_t containMatchedFake = kFALSE;
 477       Bool_t containUnmatchedFake = kFALSE;
 478       AliMUONTrack* fakeTrack = 0x0;
 479       TIter next3(fakeTrackStore->CreateIterator());
 480       while ( ( fakeTrack = static_cast<AliMUONTrack*>(next3()) ) ) {
 481         if (fakeTrack->GetMatchTrigger() > 0) containMatchedFake = kTRUE;
 482         else containUnmatchedFake = kTRUE;
 483       }
 484
 485       // fill global counters
 486       if (containMatchedFake) fEventCounters->Count(Form("event:notconnected/run:%d/trig:yes/%s", fCurrentRunNumber, selected.Data()));
 487
 488       // remove the remaining free reconstructible tracks
 489       Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
 490
 491       if (nAdditionalTracks > 0) {
 492
 493         // fill histogram and global counters
 494         ((TH1F*)fList->UncheckedAt(kNumberOfAdditionalTracks))->Fill(nAdditionalTracks);
 495         fEventCounters->Count(Form("event:additional/run:%d/trig:any/%s", fCurrentRunNumber, selected.Data()));
 496         if (!containUnmatchedFake) { // all matched
 497           fTrackCounters->Count(Form("track:additional/run:%d/trig:yes/%s/acc:unknown", fCurrentRunNumber, selected.Data()), nAdditionalTracks);
 498           fFakeTrackCounters->Count(Form("track:additional/run:%d/file:%s/event:%d/trig:yes/%s/acc:unknown", fCurrentRunNumber, fCurrentFileName.Data(), esd->GetEventNumberInFile(), selected.Data()), nAdditionalTracks);
 499           fEventCounters->Count(Form("event:additional/run:%d/trig:yes/%s", fCurrentRunNumber, selected.Data()));
 500         } else if (!containMatchedFake) { // none matched
 501           fTrackCounters->Count(Form("track:additional/run:%d/trig:no/%s/acc:unknown", fCurrentRunNumber, selected.Data()), nAdditionalTracks);
 502           fFakeTrackCounters->Count(Form("track:additional/run:%d/file:%s/event:%d/trig:no/%s/acc:unknown", fCurrentRunNumber, fCurrentFileName.Data(), esd->GetEventNumberInFile(), selected.Data()), nAdditionalTracks);
 503         } else { // mixed
 504           fTrackCounters->Count(Form("track:additional/run:%d/trig:unknown/%s/acc:unknown", fCurrentRunNumber, selected.Data()), nAdditionalTracks);
 505           fFakeTrackCounters->Count(Form("track:additional/run:%d/file:%s/event:%d/trig:unknown/%s/acc:unknown", fCurrentRunNumber, fCurrentFileName.Data(), esd->GetEventNumberInFile(), selected.Data()), nAdditionalTracks);
 506           fEventCounters->Count(Form("event:additional/run:%d/trig:yes/%s", fCurrentRunNumber, selected.Data()));
 507         }
 508
 509       }
 510
 511     }
 512
 513   }
 514
 515   delete fakeTrackStore;
 516
 517   // Post final data
 518   PostData(1, fList);
 519   PostData(2, fTrackCounters);
 520   PostData(3, fFakeTrackCounters);
 521   PostData(4, fMatchedTrackCounters);
 522   PostData(5, fEventCounters);
 523 }
 524
 525 //________________________________________________________________________
 526 void AliAnalysisTaskMuonFakes::NotifyRun()
 527 {
 528   /// Prepare processing of new run: load corresponding OCDB objects...
 529
 530   // load necessary data from OCDB
 531   AliCDBManager::Instance()->SetDefaultStorage(fRecoParamLocation.Data());
 532   AliCDBManager::Instance()->SetRun(fCurrentRunNumber);
 533   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
 534   if (!recoParam) {
 535     fRequestedStationMask = 0;
 536     fRequest2ChInSameSt45 = kFALSE;
 537     fSigmaCut = -1.;
 538     AliError("--> skip this run");
 539     return;
 540   }
 541
 542   // compute the mask of requested stations from recoParam
 543   fRequestedStationMask = 0;
 544   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) fRequestedStationMask |= ( 1 << i );
 545
 546   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
 547   fRequest2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
 548
 549   // get sigma cut from recoParam to associate clusters with TrackRefs in case the labels are not used
 550   fSigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
 551 }
 552
 553 //________________________________________________________________________
 554 void AliAnalysisTaskMuonFakes::Terminate(Option_t *)
 555 {
 556   /// Draw results to the screen and print statistics.
 557
 558   // recover output objects
 559   fList = static_cast<TObjArray*> (GetOutputData(1));
 560   if (!fList) return;
 561   fTrackCounters = static_cast<AliCounterCollection*> (GetOutputData(2));
 562   fFakeTrackCounters = static_cast<AliCounterCollection*> (GetOutputData(3));
 563   fMatchedTrackCounters = static_cast<AliCounterCollection*> (GetOutputData(4));
 564   fEventCounters = static_cast<AliCounterCollection*> (GetOutputData(5));
 565
 566   // add canvas to compare histograms
 567   fCanvases = new TObjArray(1000);
 568   fCanvases->SetOwner();
 569   TCanvas *cFakesSummary1 = new TCanvas("cFakesSummary1","cFakesSummary1",1200,600);
 570   fCanvases->AddAtAndExpand(cFakesSummary1, 0);
 571   TCanvas *cFakesSummary2 = new TCanvas("cFakesSummary2","cFakesSummary2",1200,600);
 572   fCanvases->AddAtAndExpand(cFakesSummary2, 1);
 573
 574   // display
 575   Int_t iHist1[8] = {kNumberOfClusters, kChi2PerDof, kP, kEta, kNumberOfChamberHit, kDCA, kPt, kPhi};
 576   cFakesSummary1->Divide(4,2);
 577   for (Int_t i=0; i<8; i++) {
 578     cFakesSummary1->cd(i+1);
 579     cFakesSummary1->GetPad(i+1)->SetLogy();
 580     ((TH1F*)fList->UncheckedAt(iHist1[i]))->SetMinimum(0.5);
 581     ((TH1F*)fList->UncheckedAt(iHist1[i]))->DrawCopy();
 582     ((TH1F*)fList->UncheckedAt(iHist1[i]+1))->SetLineColor(4);
 583     ((TH1F*)fList->UncheckedAt(iHist1[i]+1))->DrawCopy("sames");
 584     ((TH1F*)fList->UncheckedAt(iHist1[i]+2))->SetLineColor(2);
 585     ((TH1F*)fList->UncheckedAt(iHist1[i]+2))->SetFillColor(2);
 586     ((TH1F*)fList->UncheckedAt(iHist1[i]+2))->SetFillStyle(3017);
 587     ((TH1F*)fList->UncheckedAt(iHist1[i]+2))->DrawCopy("sames");
 588   }
 589
 590   Int_t iHist2[2] = {kChi2PerDofVsNClusters, kChi2PerDofVsNChamberHit};
 591   cFakesSummary2->Divide(2);
 592   for (Int_t i=0; i<2; i++) {
 593     cFakesSummary2->cd(i+1);
 594     ((TH2F*)fList->UncheckedAt(iHist2[i]+1))->SetMarkerColor(4);
 595     ((TH2F*)fList->UncheckedAt(iHist2[i]+1))->DrawCopy();
 596     ((TH2F*)fList->UncheckedAt(iHist2[i]+2))->SetMarkerColor(2);
 597     ((TH2F*)fList->UncheckedAt(iHist2[i]+2))->SetMarkerStyle(7);
 598     ((TH2F*)fList->UncheckedAt(iHist2[i]+2))->DrawCopy("sames");
 599   }
 600
 601   // print
 602   if (fTrackCounters && fFakeTrackCounters && fMatchedTrackCounters && fEventCounters) {
 603     printf("\nGlobal statistics of reconstructed tracks matched or not with the trigger:\n");
 604     fTrackCounters->Print("track/trig");
 605     printf("\nGlobal statistics of pathological tracks matched or not with the trigger:\n");
 606     fFakeTrackCounters->Print("track/trig");
 607     printf("\nDetailled statistics of tracks matched per label vs position:\n");
 608     fMatchedTrackCounters->Print("label/position");
 609     printf("\nGlobal statistics of events containing pathological tracks:\n");
 610     fEventCounters->Print("event/trig");
 611   }
 612
 613   printf("\nREMINDER: results are relevent provided that you use the same recoParams as for the reconstruction\n");
 614 }
 615
 616 //________________________________________________________________________
 617 Int_t AliAnalysisTaskMuonFakes::RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore)
 618 {
 619   /// loop over reconstructible TrackRef not associated with reconstructed track:
 620   /// for each of them, find and remove the most connected the fake track, if any,
 621   /// and fill the histograms with the fraction of connected clusters.
 622   /// Return the number of reconstructible track not connected to any fake
 623
 624   Int_t nFreeMissingTracks = 0;
 625
 626   // loop over trackRefs
 627   TIter next(trackRefStore.CreateIterator());
 628   AliMUONTrack* trackRef;
 629   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
 630
 631     // skip not reconstructible trackRefs
 632     if (!trackRef->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) continue;
 633
 634     Int_t label = trackRef->GetUniqueID();
 635
 636     // look for the most connected fake track
 637     AliMUONTrack *connectedFake = 0x0;
 638     Double_t fractionOfConnectedClusters = 0.;
 639     TIter next2(fakeTrackStore.CreateIterator());
 640     AliMUONTrack* fakeTrack;
 641     while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
 642
 643       // get the number of connected clusters
 644       Int_t nConnectedClusters = 0;
 645       if (fUseLabel) { // by using the MC label
 646         for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
 647           if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
 648             nConnectedClusters++;
 649       } else { // by comparing cluster/TrackRef positions
 650         Bool_t compTrack[10];
 651         nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, fSigmaCut, compTrack);
 652       }
 653
 654       // skip non-connected fake tracks
 655       if (nConnectedClusters == 0) continue;
 656
 657       // check if it is the most connected fake track
 658       Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters());
 659       if (f > fractionOfConnectedClusters) {
 660         connectedFake = fakeTrack;
 661         fractionOfConnectedClusters = f;
 662       }
 663
 664     }
 665
 666     if (connectedFake) {
 667
 668       // find the corresponding ESD MUON track
 669       AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 670       TIter next3(static_cast<TClonesArray*>(esd->FindListObject("MuonTracks")));
 671       AliESDMuonTrack* esdTrack = 0x0;
 672       while ((esdTrack = static_cast<AliESDMuonTrack*>(next3())) && esdTrack->GetUniqueID() != connectedFake->GetUniqueID()) {}
 673       if (!esdTrack) {
 674         AliError("unable to find the corresponding ESD track???");
 675         continue;
 676       }
 677
 678       // trigger condition
 679       TString trig = (esdTrack->ContainTriggerData()) ? "trig:yes" : "trig:no";
 680
 681       // acceptance condition
 682       Double_t thetaTrackAbsEnd = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 683       Double_t eta = esdTrack->Eta();
 684       TString acc = (thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 9. && eta >= -4. && eta <= -2.5) ? "acc:in" : "acc:out";
 685
 686       // check physics selection
 687       TString selected = (fInputHandler && fInputHandler->IsEventSelected()) ? "selected:yes" : "selected:no";
 688
 689       // fill histogram and counters
 690       ((TH1F*)fList->UncheckedAt(kFractionOfConnectedClusters))->Fill(fractionOfConnectedClusters);
 691       fTrackCounters->Count(Form("track:connected/run:%d/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data()));
 692       fFakeTrackCounters->Count(Form("track:connected/run:%d/file:%s/event:%d/%s/%s/%s", fCurrentRunNumber, fCurrentFileName.Data(),
 693                                      esd->GetEventNumberInFile(), trig.Data(), selected.Data(), acc.Data()));
 694
 695       // remove the most connected fake track
 696       fakeTrackStore.Remove(*connectedFake);
 697
 698     } else nFreeMissingTracks++;
 699
 700   }
 701
 702   return nFreeMissingTracks;
 703
 704 }
 705