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