Analysis Taks for RecoCheck muon code (Bogdan, Nicole)
[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"
7#include "AliMagFMaps.h"
8#include "AliTracker.h"
9
10#include "AliAnalysisTask.h"
11#include "AliAnalysisManager.h"
12
13#include "AliESDEvent.h"
14#include "AliESDInputHandler.h"
15#include "AliESDMuonTrack.h"
16
17#include "AliMCEventHandler.h"
18#include "AliMCEvent.h"
19
20#include "AliAnalysisTaskRecoCheck.h"
21
22#include "AliMUONTrackLight.h"
23#include "AliMUONPairLight.h"
24#include "AliMUONVTrackStore.h"
25#include "AliMUONTrack.h"
26#include "AliMUONRecoCheck.h"
27#include "AliMUONTrackExtrap.h"
28
29// analysis task for decoding reconstructed tracks and kinematics (AliMUONRecoCheck)
30// Authors: Bogdan Vulpescu
31
32ClassImp(AliAnalysisTaskRecoCheck)
33
34//________________________________________________________________________
35AliAnalysisTaskRecoCheck::AliAnalysisTaskRecoCheck(const char *name)
36 : AliAnalysisTask(name, ""), fESDEvent(0), fTree(0), fArray1Mu(0), fArray2Mu(0), fL3Current(30000.0)
37{
38 // Constructor
39
40 // Define input and output slots here
41 // Input slot #0 works with a TChain
42 DefineInput(0, TChain::Class());
43 // Output slot #0 writes into a TTree container
44 DefineOutput(0, TTree::Class());
45
46}
47
48//________________________________________________________________________
49void AliAnalysisTaskRecoCheck::ConnectInputData(Option_t *)
50{
51 // Connect ESD or AOD here
52 // Called once
53
54 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
55 if (!tree) {
56 Printf("ERROR: Could not read chain from input slot 0");
57 } else {
58 tree->SetBranchStatus("*", kFALSE);
59 tree->SetBranchStatus("MuonTracks.*", kTRUE);
60 tree->SetBranchStatus("AliESDHeader.*", kTRUE);
61
62 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
63
64 if (!esdH) {
65 Printf("ERROR: Could not get ESDInputHandler");
66 } else
67 fESDEvent = esdH->GetEvent();
68 }
69
70 // calculate the filed map in the L3 magnet using the current value
71 AliMagFMaps* field = 0x0;
72 if (fL3Current == 30000.0) {
73 field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
74 }
75
76 // set the tracker field map
77 AliTracker::SetFieldMap(field, kFALSE);
78 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
79
80}
81
82//________________________________________________________________________
83void AliAnalysisTaskRecoCheck::CreateOutputObjects()
84{
85 // Create output tree with branches of arrays
86 // Called once
87
88 fArray1Mu = new TClonesArray("AliMUONTrackLight",0);
89 fArray2Mu = new TClonesArray("AliMUONPairLight",0);
90 fTree = new TTree("tree","tree",0);
91 fTree->Branch("muons", &fArray1Mu);
92 fTree->Branch("dimuons",&fArray2Mu);
93 fTree->SetDirectory(0);
94
95}
96
97//________________________________________________________________________
98void AliAnalysisTaskRecoCheck::Exec(Option_t *)
99{
100 // Main loop
101 // Called for each event
102
103 // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
104 // This handler can return the current MC event
105
106 Int_t nTracks = fESDEvent->GetNumberOfMuonTracks();
107 if (nTracks == 0) return;
108
109 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
110 if (!eventHandler) {
111 Printf("ERROR: Could not retrieve MC event handler");
112 return;
113 }
114
115 AliMCEvent* mcEvent = eventHandler->MCEvent();
116 if (!mcEvent) {
117 Printf("ERROR: Could not retrieve MC event");
118 return;
119 }
120
121 Int_t eventNumber = fESDEvent->GetEventNumberInFile();
122
123 fArray1Mu->Clear();
124 fArray2Mu->Clear();
125
126 TLorentzVector v;
127
128 AliMUONRecoCheck *rc = new AliMUONRecoCheck(fESDEvent, eventHandler);
129
130 AliMUONVTrackStore* trackRefs = rc->TrackRefs(eventNumber);
131 AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(eventNumber);
132 AliStack* pstack = mcEvent->Stack();
133
134 TIter next(recoTracks->CreateIterator());
135 AliMUONTrack* trackReco = NULL;
136
137 Int_t nTrackReco = recoTracks->GetSize();
138 Int_t nTracksESD = fESDEvent->GetNumberOfMuonTracks();
139 //printf("Tracks in recoTracks (%d) and in ESD (%d).\n", nTrackReco, nTracksESD);
140
141 if (nTrackReco != nTracksESD) printf ("Tracks in recoTracks (%d) and in ESD (%d) do not match!\n", nTrackReco, nTracksESD);
142
143 Int_t nreftracks = 0;
144 Int_t itrRec = 0;
145 while ( (trackReco = static_cast<AliMUONTrack*>(next())) != NULL ) {
146 // assign parameters concerning the reconstructed tracks
147 AliMUONTrackLight muLight;
148
149 muLight.FillFromESD(fESDEvent->GetMuonTrack(itrRec));
150
151 // find the reference track and store further information
152 TParticle *part = muLight.FindRefTrack(trackReco, trackRefs, pstack);
153 if (part) {
154 v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
155 muLight.SetPGen(v);
156 muLight.FillMuonHistory(pstack, part);
157 //muLight.PrintInfo("A");
158 //store the referenced track in the muonArray:
159 TClonesArray &muons = *fArray1Mu;
160 new (muons[nreftracks++]) AliMUONTrackLight(muLight);
161 }
162
163 itrRec++;
164 } // end reco track
165
166 // now loop over muon pairs to build dimuons
167 Int_t nmuons = fArray1Mu->GetEntriesFast();
168 Int_t ndimuons = 0;
169 for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) {
170 AliMUONTrackLight* mu1 = (AliMUONTrackLight*) fArray1Mu->At(itrRec1);
171 for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
172 AliMUONTrackLight* mu2 = (AliMUONTrackLight*) fArray1Mu->At(itrRec2);
173 AliMUONPairLight dimuLight;
174 dimuLight.SetMuons(*mu1, *mu2);
175 //dimuLight.PrintInfo("A");
176 TClonesArray &dimuons = *fArray2Mu;
177 new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
178 }
179 }
180
181 delete rc;
182
183 // Post output data.
184 if (nreftracks != 0) {
185 fTree->Fill();
186 }
187
188 PostData(0, fTree);
189
190
191}
192
193//________________________________________________________________________
194void AliAnalysisTaskRecoCheck::Terminate(Option_t *)
195{
196 // the end
197
198}