Coding conventions
[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskRecoCheck.cxx
CommitLineData
27de2dfb 1
2/* $Id$ */
3
d71eddb1 4#include "TChain.h"
5#include "TClonesArray.h"
6#include "TLorentzVector.h"
7#include "TIterator.h"
8
9#include "AliStack.h"
d71eddb1 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"
f202486b 27#include "AliMUONESDInterface.h"
d71eddb1 28#include "AliMUONRecoCheck.h"
d71eddb1 29
30// analysis task for decoding reconstructed tracks and kinematics (AliMUONRecoCheck)
31// Authors: Bogdan Vulpescu
32
33ClassImp(AliAnalysisTaskRecoCheck)
34
35//________________________________________________________________________
36AliAnalysisTaskRecoCheck::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//________________________________________________________________________
50void 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
d71eddb1 71}
72
73//________________________________________________________________________
74void 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//________________________________________________________________________
89void 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
d71eddb1 112 fArray1Mu->Clear();
113 fArray2Mu->Clear();
114
115 TLorentzVector v;
116
117 AliMUONRecoCheck *rc = new AliMUONRecoCheck(fESDEvent, eventHandler);
118
f202486b 119 AliMUONVTrackStore* trackRefs = rc->TrackRefs(-1);
d71eddb1 120 AliStack* pstack = mcEvent->Stack();
121
f202486b 122 // loop over ESD tracks
d71eddb1 123 Int_t nreftracks = 0;
f202486b 124 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
125
126 AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
d71eddb1 127
f202486b 128 // skip ghosts
129 if (!esdTrack->ContainTrackerData()) continue;
d71eddb1 130
f202486b 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());
d71eddb1 148 v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
f202486b 149 muLight->SetPGen(v);
150 muLight->FillMuonHistory(pstack, part);
151
152 }
d71eddb1 153
f202486b 154 } // end esd track
155
d71eddb1 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//________________________________________________________________________
184void AliAnalysisTaskRecoCheck::Terminate(Option_t *)
185{
186 // the end
187
188}