PWG3/muondep/AliAnalysisTaskRecoCheck.cxx

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