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