Coding conventions
[u/mrichter/AliRoot.git] / 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 }