Fixing an error in the computation of the efficiency of a whole chamber (Matthieu)
[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskRecoCheck.cxx
CommitLineData
d71eddb1 1#include "TChain.h"
2#include "TClonesArray.h"
3#include "TLorentzVector.h"
4#include "TIterator.h"
5
6#include "AliStack.h"
d71eddb1 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"
f202486b 24#include "AliMUONESDInterface.h"
d71eddb1 25#include "AliMUONRecoCheck.h"
d71eddb1 26
27// analysis task for decoding reconstructed tracks and kinematics (AliMUONRecoCheck)
28// Authors: Bogdan Vulpescu
29
30ClassImp(AliAnalysisTaskRecoCheck)
31
32//________________________________________________________________________
33AliAnalysisTaskRecoCheck::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//________________________________________________________________________
47void 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
d71eddb1 68}
69
70//________________________________________________________________________
71void 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//________________________________________________________________________
86void 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
d71eddb1 109 fArray1Mu->Clear();
110 fArray2Mu->Clear();
111
112 TLorentzVector v;
113
114 AliMUONRecoCheck *rc = new AliMUONRecoCheck(fESDEvent, eventHandler);
115
f202486b 116 AliMUONVTrackStore* trackRefs = rc->TrackRefs(-1);
d71eddb1 117 AliStack* pstack = mcEvent->Stack();
118
f202486b 119 // loop over ESD tracks
d71eddb1 120 Int_t nreftracks = 0;
f202486b 121 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
122
123 AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
d71eddb1 124
f202486b 125 // skip ghosts
126 if (!esdTrack->ContainTrackerData()) continue;
d71eddb1 127
f202486b 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());
d71eddb1 145 v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
f202486b 146 muLight->SetPGen(v);
147 muLight->FillMuonHistory(pstack, part);
148
149 }
d71eddb1 150
f202486b 151 } // end esd track
152
d71eddb1 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//________________________________________________________________________
181void AliAnalysisTaskRecoCheck::Terminate(Option_t *)
182{
183 // the end
184
185}