PWG3/muondep/AliAnalysisTaskRecoCheck.cxx

   1
   2 /* $Id$ */
   3
   4 #include "TChain.h"
   5 #include "TClonesArray.h"
   6 #include "TLorentzVector.h"
   7 #include "TIterator.h"
   8
   9 #include "AliStack.h"
  10
  11 #include "AliAnalysisTask.h"
  12 #include "AliAnalysisManager.h"
  13
  14 #include "AliESDEvent.h"
  15 #include "AliESDInputHandler.h"
  16 #include "AliESDMuonTrack.h"
  17
  18 #include "AliMCEventHandler.h"
  19 #include "AliMCEvent.h"
  20
  21 #include "AliAnalysisTaskRecoCheck.h"
  22
  23 #include "AliMUONTrackLight.h"
  24 #include "AliMUONPairLight.h"
  25 #include "AliMUONVTrackStore.h"
  26 #include "AliMUONTrack.h"
  27 #include "AliMUONESDInterface.h"
  28 #include "AliMUONRecoCheck.h"
  29
  30 // analysis task for decoding reconstructed tracks and kinematics (AliMUONRecoCheck)
  31 // Authors: Bogdan Vulpescu
  32
  33 ClassImp(AliAnalysisTaskRecoCheck)
  34
  35 //________________________________________________________________________
  36 AliAnalysisTaskRecoCheck::AliAnalysisTaskRecoCheck(const char *name)
  37   : AliAnalysisTask(name, ""), fESDEvent(0), fTree(0), fArray1Mu(0), fArray2Mu(0), fL3Current(30000.0)
  38 {
  39   // Constructor
  40
  41   // Define input and output slots here
  42   // Input slot #0 works with a TChain
  43   DefineInput(0, TChain::Class());
  44   // Output slot #0 writes into a TTree container
  45   DefineOutput(0, TTree::Class());
  46
  47 }
  48
  49 //________________________________________________________________________
  50 void AliAnalysisTaskRecoCheck::ConnectInputData(Option_t *)
  51 {
  52   // Connect ESD or AOD here
  53   // Called once
  54
  55   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  56   if (!tree) {
  57     Printf("ERROR: Could not read chain from input slot 0");
  58   } else {
  59     tree->SetBranchStatus("*", kFALSE);
  60     tree->SetBranchStatus("MuonTracks.*", kTRUE);
  61     tree->SetBranchStatus("AliESDHeader.*", kTRUE);
  62
  63     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  64
  65     if (!esdH) {
  66       Printf("ERROR: Could not get ESDInputHandler");
  67     } else
  68       fESDEvent = esdH->GetEvent();
  69   }
  70
  71 }
  72
  73 //________________________________________________________________________
  74 void AliAnalysisTaskRecoCheck::CreateOutputObjects()
  75 {
  76   // Create output tree with branches of arrays
  77   // Called once
  78
  79   fArray1Mu = new TClonesArray("AliMUONTrackLight",0);
  80   fArray2Mu = new TClonesArray("AliMUONPairLight",0);
  81   fTree = new TTree("tree","tree",0);
  82   fTree->Branch("muons",  &fArray1Mu);
  83   fTree->Branch("dimuons",&fArray2Mu);
  84   fTree->SetDirectory(0);
  85
  86 }
  87
  88 //________________________________________________________________________
  89 void AliAnalysisTaskRecoCheck::Exec(Option_t *)
  90 {
  91   // Main loop
  92   // Called for each event
  93
  94   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
  95   // This handler can return the current MC event
  96
  97   Int_t nTracks = fESDEvent->GetNumberOfMuonTracks();
  98   if (nTracks == 0) return;
  99
 100   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 101   if (!eventHandler) {
 102     Printf("ERROR: Could not retrieve MC event handler");
 103     return;
 104   }
 105
 106   AliMCEvent* mcEvent = eventHandler->MCEvent();
 107   if (!mcEvent) {
 108      Printf("ERROR: Could not retrieve MC event");
 109      return;
 110   }
 111
 112   fArray1Mu->Clear();
 113   fArray2Mu->Clear();
 114
 115   TLorentzVector v;
 116
 117   AliMUONRecoCheck *rc = new AliMUONRecoCheck(fESDEvent, eventHandler);
 118
 119   AliMUONVTrackStore* trackRefs = rc->TrackRefs(-1);
 120   AliStack* pstack = mcEvent->Stack();
 121
 122   // loop over ESD tracks
 123   Int_t nreftracks = 0;
 124   for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 125
 126     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
 127
 128     // skip ghosts
 129     if (!esdTrack->ContainTrackerData()) continue;
 130
 131     // convert ESD track to MUON track (without recomputing track parameters at each clusters)
 132     AliMUONTrack muonTrack;
 133     AliMUONESDInterface::ESDToMUON(*esdTrack, muonTrack, kFALSE);
 134
 135     // try to match the reconstructed track with a simulated one
 136     Int_t nMatchClusters = 0;
 137     AliMUONTrack* matchedTrackRef = rc->FindCompatibleTrack(muonTrack, *trackRefs, nMatchClusters, kFALSE, 10.);
 138
 139     if (matchedTrackRef) {
 140
 141       //store new referenced track in the muonArray with parameters from the reconstructed tracks
 142       AliMUONTrackLight* muLight = new ((*fArray1Mu)[nreftracks++]) AliMUONTrackLight(esdTrack);
 143
 144       // store further information related to the simulated track
 145       muLight->SetTrackPythiaLine(matchedTrackRef->GetUniqueID());
 146       TParticle *part = pstack->Particle(matchedTrackRef->GetUniqueID());
 147       muLight->SetTrackPDGCode(part->GetPdgCode());
 148       v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
 149       muLight->SetPGen(v);
 150       muLight->FillMuonHistory(pstack, part);
 151
 152     }
 153
 154   } // end esd track
 155
 156   // now loop over muon pairs to build dimuons
 157   Int_t nmuons = fArray1Mu->GetEntriesFast();
 158   Int_t ndimuons = 0;
 159   for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) {
 160     AliMUONTrackLight* mu1 = (AliMUONTrackLight*) fArray1Mu->At(itrRec1);
 161     for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
 162       AliMUONTrackLight* mu2 = (AliMUONTrackLight*) fArray1Mu->At(itrRec2);
 163       AliMUONPairLight dimuLight;
 164       dimuLight.SetMuons(*mu1, *mu2);
 165       //dimuLight.PrintInfo("A");
 166       TClonesArray &dimuons = *fArray2Mu;
 167       new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
 168     }
 169   }
 170
 171   delete rc;
 172
 173   // Post output data.
 174   if (nreftracks != 0) {
 175     fTree->Fill();
 176   }
 177
 178   PostData(0, fTree);
 179
 180
 181 }
 182
 183 //________________________________________________________________________
 184 void AliAnalysisTaskRecoCheck::Terminate(Option_t *)
 185 {
 186   // the end
 187
 188 }