The present commit corresponds to an important change in the way the
[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 #include "AliMagF.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
32 ClassImp(AliAnalysisTaskRecoCheck)
33
34 //________________________________________________________________________
35 AliAnalysisTaskRecoCheck::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 //________________________________________________________________________
49 void 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   AliMUONTrackExtrap::SetField();
72
73 }
74
75 //________________________________________________________________________
76 void AliAnalysisTaskRecoCheck::CreateOutputObjects() 
77 {
78   // Create output tree with branches of arrays
79   // Called once
80
81   fArray1Mu = new TClonesArray("AliMUONTrackLight",0); 
82   fArray2Mu = new TClonesArray("AliMUONPairLight",0); 
83   fTree = new TTree("tree","tree",0); 
84   fTree->Branch("muons",  &fArray1Mu); 
85   fTree->Branch("dimuons",&fArray2Mu); 
86   fTree->SetDirectory(0);
87   
88 }
89
90 //________________________________________________________________________
91 void AliAnalysisTaskRecoCheck::Exec(Option_t *) 
92 {
93   // Main loop
94   // Called for each event
95   
96   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
97   // This handler can return the current MC event
98   
99   Int_t nTracks = fESDEvent->GetNumberOfMuonTracks();
100   if (nTracks == 0) return;
101   
102   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
103   if (!eventHandler) {
104     Printf("ERROR: Could not retrieve MC event handler");
105     return;
106   }
107
108   AliMCEvent* mcEvent = eventHandler->MCEvent();
109   if (!mcEvent) {
110      Printf("ERROR: Could not retrieve MC event");
111      return;
112   }
113
114   Int_t eventNumber = fESDEvent->GetEventNumberInFile();
115
116   fArray1Mu->Clear();
117   fArray2Mu->Clear();
118
119   TLorentzVector v; 
120
121   AliMUONRecoCheck *rc = new AliMUONRecoCheck(fESDEvent, eventHandler);
122
123   AliMUONVTrackStore* trackRefs = rc->TrackRefs(eventNumber);
124   AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(eventNumber);
125   AliStack* pstack = mcEvent->Stack();
126   
127   TIter next(recoTracks->CreateIterator());
128   AliMUONTrack* trackReco = NULL;
129   
130   Int_t nTrackReco = recoTracks->GetSize();
131   Int_t nTracksESD = fESDEvent->GetNumberOfMuonTracks();
132   //printf("Tracks in recoTracks (%d) and in ESD (%d).\n", nTrackReco, nTracksESD); 
133
134   if (nTrackReco != nTracksESD) printf ("Tracks in recoTracks (%d) and in ESD (%d) do not match!\n", nTrackReco, nTracksESD);
135   
136   Int_t nreftracks = 0;
137   Int_t itrRec = 0;
138   while ( (trackReco = static_cast<AliMUONTrack*>(next())) != NULL ) {
139     // assign parameters concerning the reconstructed tracks
140     AliMUONTrackLight muLight;
141     
142     muLight.FillFromESD(fESDEvent->GetMuonTrack(itrRec));
143     
144     // find the reference track and store further information   
145     TParticle *part = muLight.FindRefTrack(trackReco, trackRefs, pstack);
146     if (part) { 
147       v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
148       muLight.SetPGen(v); 
149       muLight.FillMuonHistory(pstack, part);
150       //muLight.PrintInfo("A");
151       //store the referenced track in the muonArray:
152       TClonesArray &muons = *fArray1Mu;
153       new (muons[nreftracks++]) AliMUONTrackLight(muLight);
154     } 
155     
156     itrRec++;
157   }  // end reco track
158
159   // now loop over muon pairs to build dimuons
160   Int_t nmuons = fArray1Mu->GetEntriesFast(); 
161   Int_t ndimuons = 0; 
162   for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) { 
163     AliMUONTrackLight* mu1 = (AliMUONTrackLight*) fArray1Mu->At(itrRec1); 
164     for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
165       AliMUONTrackLight* mu2 = (AliMUONTrackLight*) fArray1Mu->At(itrRec2); 
166       AliMUONPairLight dimuLight;
167       dimuLight.SetMuons(*mu1, *mu2);
168       //dimuLight.PrintInfo("A");
169       TClonesArray &dimuons = *fArray2Mu;
170       new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
171     }
172   }
173     
174   delete rc;
175
176   // Post output data.
177   if (nreftracks != 0) {
178     fTree->Fill();
179   }
180   
181   PostData(0, fTree);
182   
183   
184 }      
185
186 //________________________________________________________________________
187 void AliAnalysisTaskRecoCheck::Terminate(Option_t *) 
188 {
189   // the end
190
191 }