TENDER becomes Tender, removing .so
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskResolution.cxx
CommitLineData
72395bd9 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Authors: Friederike Bock *
5 * Version 1.0 *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16////////////////////////////////////////////////
17//---------------------------------------------
18// QA Task for V0 Reader V1
19//---------------------------------------------
20////////////////////////////////////////////////
21
22#include "AliAnalysisTaskResolution.h"
23#include "TChain.h"
24#include "AliAnalysisManager.h"
25#include "TParticle.h"
26#include "TVectorF.h"
27#include "AliPIDResponse.h"
28#include "AliESDtrackCuts.h"
29#include "TFile.h"
30
31class iostream;
32
33using namespace std;
34
35ClassImp(AliAnalysisTaskResolution)
36
37//________________________________________________________________________
38AliAnalysisTaskResolution::AliAnalysisTaskResolution() : AliAnalysisTaskSE(),
63674d67 39 fV0Reader(NULL),
40 fConversionGammas(NULL),
344100c4 41 fEventCuts(NULL),
63674d67 42 fConversionCuts(NULL),
43 fTreeEvent(NULL),
44 fTreeResolution(NULL),
45 fPrimVtxZ(0.),
46 fNContrVtx(0),
47 fNESDtracksEta09(0),
48 fNESDtracksEta0914(0),
49 fNESDtracksEta14(0),
50 fGammaRecCoords(5),
51 fGammaMCCoords(5),
52 fChi2ndf(0),
53 fIsHeavyIon(0),
54 fIsMC(kFALSE),
55 fOutputList(NULL),
56 fEventList(NULL),
57 fResolutionList(NULL),
58 fESDEvent(NULL),
59 fMCEvent(NULL)
72395bd9 60{
61
62
63}
64
72395bd9 65//________________________________________________________________________
66AliAnalysisTaskResolution::AliAnalysisTaskResolution(const char *name) : AliAnalysisTaskSE(name),
63674d67 67 fV0Reader(NULL),
68 fConversionGammas(NULL),
344100c4 69 fEventCuts(NULL),
63674d67 70 fConversionCuts(NULL),
71 fTreeEvent(NULL),
72 fTreeResolution(NULL),
73 fPrimVtxZ(0.),
74 fNContrVtx(0),
75 fNESDtracksEta09(0),
76 fNESDtracksEta0914(0),
77 fNESDtracksEta14(0),
78 fGammaRecCoords(5),
79 fGammaMCCoords(5),
80 fChi2ndf(0),
81 fIsHeavyIon(0),
82 fIsMC(kFALSE),
83 fOutputList(NULL),
84 fEventList(NULL),
85 fResolutionList(NULL),
86 fESDEvent(NULL),
87 fMCEvent(NULL)
72395bd9 88{
63674d67 89 // Default constructor
72395bd9 90
63674d67 91 DefineInput(0, TChain::Class());
92 DefineOutput(1, TList::Class());
72395bd9 93}
94
95//________________________________________________________________________
96AliAnalysisTaskResolution::~AliAnalysisTaskResolution()
97{
63674d67 98 // default deconstructor
72395bd9 99
100}
101//________________________________________________________________________
102void AliAnalysisTaskResolution::UserCreateOutputObjects()
103{
63674d67 104 // Create User Output Objects
72395bd9 105
63674d67 106 if(fOutputList != NULL){
107 delete fOutputList;
108 fOutputList = NULL;
109 }
110 if(fOutputList == NULL){
111 fOutputList = new TList();
112 fOutputList->SetOwner(kTRUE);
113 }
114
115 fEventList = new TList();
116 fEventList->SetName("EventList");
117 fEventList->SetOwner(kTRUE);
118 fOutputList->Add(fEventList);
119
120 fTreeEvent = new TTree("Event","Event");
121 fTreeEvent->Branch("primVtxZ",&fPrimVtxZ,"fPrimVtxZ/F");
122 fTreeEvent->Branch("nContrVtx",&fNContrVtx,"fNContrVtx/I");
123 fTreeEvent->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
63674d67 124 fTreeEvent->Branch("nGoodTracksEta14",&fNESDtracksEta14,"fNESDtracksEta14/I");
125 fEventList->Add(fTreeEvent);
126
127 if (fIsMC){
128 fResolutionList = new TList();
129 fResolutionList->SetName("ResolutionList");
130 fResolutionList->SetOwner(kTRUE);
131 fOutputList->Add(fResolutionList);
132
97146ba4 133 fTreeResolution = new TTree("Resolution","Resolution");
134 fTreeResolution->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
63674d67 135 fTreeResolution->Branch("RecCoords",&fGammaRecCoords);
136 fTreeResolution->Branch("MCCoords",&fGammaMCCoords);
137 fTreeResolution->Branch("chi2ndf",&fChi2ndf,"fChi2ndf/F");
138 fResolutionList->Add(fTreeResolution);
139 }
140 PostData(1, fOutputList);
72395bd9 141}
142
143//________________________________________________________________________
144void AliAnalysisTaskResolution::UserExec(Option_t *){
145
51ead2dc 146 fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
72395bd9 147
344100c4 148 Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
51ead2dc 149 if(eventQuality != 0){// Event Not Accepted
150 return;
151 }
152 fESDEvent = (AliESDEvent*) InputEvent();
153 if (fESDEvent==NULL) return;
72395bd9 154 if(MCEvent()){
155 fMCEvent = MCEvent();
156 }
51ead2dc 157
63674d67 158 if(MCEvent()){
51ead2dc 159 // Process MC Particle
344100c4 160 if(fEventCuts->GetSignalRejection() != 0){
51ead2dc 161// if(fESDEvent->IsA()==AliESDEvent::Class()){
344100c4 162 fEventCuts->GetNotRejectedParticles( fEventCuts->GetSignalRejection(),
163 fEventCuts->GetAcceptedHeader(),
164 fMCEvent);
51ead2dc 165// }
166// else if(fInputEvent->IsA()==AliAODEvent::Class()){
344100c4 167// fEventCuts->GetNotRejectedParticles( fEventCuts->GetSignalRejection(),
168// fEventCuts->GetAcceptedHeader(),
169// fInputEvent);
51ead2dc 170// }
171 }
172 }
173
174
344100c4 175 if(fIsHeavyIon > 0 && !fEventCuts->IsCentralitySelected(fESDEvent)) return;
51ead2dc 176 fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
177 fNESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
178 fNESDtracksEta14 = fNESDtracksEta09 + fNESDtracksEta0914;
179 if(fESDEvent){
180 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
181 fNContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
182 } else {
183 fNContrVtx = 0;
184 }
185// else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
186// if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
187// fNContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
188// } else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
189// fNContrVtx = 0;
190// }
191// }
192 }
193 fPrimVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
194
97146ba4 195// if (fIsHeavyIon == 2){
196// if (!(fNESDtracksEta09 > 20 && fNESDtracksEta09 < 80)) return;
197// }
51ead2dc 198
199
200 if (fTreeEvent){
201 fTreeEvent->Fill();
202 }
203
204 fConversionGammas=fV0Reader->GetReconstructedGammas();
205 ProcessPhotons();
206 PostData(1, fOutputList);
72395bd9 207}
208
209
210///________________________________________________________________________
211void AliAnalysisTaskResolution::ProcessPhotons(){
212
213 // Fill Histograms for QA and MC
63674d67 214 for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast();firstGammaIndex++){
215 AliAODConversionPhoton *gamma=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(firstGammaIndex));
216 if (gamma ==NULL) continue;
217 if(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;
97146ba4 218 fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
63674d67 219 fGammaRecCoords(0) = gamma->GetPhotonPt();
220 fGammaRecCoords(1) = gamma->GetPhotonPhi();
221 fGammaRecCoords(2) = gamma->GetPhotonEta();
222 fGammaRecCoords(3) = gamma->GetConversionRadius();
223 fGammaRecCoords(4) = gamma->GetConversionZ();
224 fChi2ndf = gamma->GetChi2perNDF();
225 if(MCEvent()){
72395bd9 226// cout << "generating MC stack"<< endl;
227 AliStack *fMCStack = fMCEvent->Stack();
228 if (!fMCStack) continue;
229 TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
230 TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
231// cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
344100c4 232 if(fMCStack && fEventCuts->GetSignalRejection() != 0){
233 Int_t isPosFromMBHeader = fEventCuts->IsParticleFromBGEvent(gamma->GetMCLabelPositive(), fMCStack, fESDEvent);
234 Int_t isNegFromMBHeader = fEventCuts->IsParticleFromBGEvent(gamma->GetMCLabelNegative(), fMCStack, fESDEvent);
51ead2dc 235 if( (isNegFromMBHeader < 1) || (isPosFromMBHeader < 1)) continue;
236 }
72395bd9 237
238 if(posDaughter == NULL || negDaughter == NULL){
239 continue;
240 } else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
241 continue;
242 } else {
243// cout << "same mother" << endl;
244 Int_t pdgCodePos;
245 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
246 Int_t pdgCodeNeg;
247 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
248// cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
249 Int_t pdgCode;
250 if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
251// cout << "PDG code: " << pdgCode << endl;
252 if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
253 continue;
254 else if ( !(pdgCodeNeg==pdgCodePos)){
255 TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
256 if(pdgCode == 111)
257 continue;
258 else if (pdgCode == 221)
259 continue;
260 else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
261 if(pdgCode == 22){
63674d67 262 fGammaMCCoords(0) = truePhotonCanditate->Pt();
263 fGammaMCCoords(1) = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
264 fGammaMCCoords(2) = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
265 fGammaMCCoords(3) = gamma->GetNegativeMCDaughter(fMCStack)->R();
266 fGammaMCCoords(4) = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
267
268 if (fTreeResolution){
269 fTreeResolution->Fill();
72395bd9 270 }
271 }
272 } else continue; //garbage
273 } else continue; //garbage
274 }
63674d67 275 }
72395bd9 276 }
277}
278
279//________________________________________________________________________
280Int_t AliAnalysisTaskResolution::CountTracks09(){
51ead2dc 281 Int_t fNumberOfESDTracks = 0;
282 if(fInputEvent->IsA()==AliESDEvent::Class()){
283 // Using standard function for setting Cuts
284
285 AliStack *fMCStack = NULL;
63674d67 286 if (MCEvent()){
51ead2dc 287 fMCStack= fMCEvent->Stack();
288 if (!fMCStack) return 0;
289 }
290
291 Bool_t selectPrimaries=kTRUE;
292 AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
293 EsdTrackCuts->SetMaxDCAToVertexZ(2);
294 EsdTrackCuts->SetEtaRange(-0.9, 0.9);
295 EsdTrackCuts->SetPtRange(0.15);
296
297 for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
298 AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
299 if(!curTrack) continue;
300 if(EsdTrackCuts->AcceptTrack(curTrack) ){
301 if (fMCEvent){
344100c4 302 if(fMCStack && fEventCuts->GetSignalRejection() != 0){
303 Int_t isFromMBHeader = fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
51ead2dc 304 if( (isFromMBHeader < 1) ) continue;
305 }
306 }
307 fNumberOfESDTracks++;
308 }
309 }
310 delete EsdTrackCuts;
311 EsdTrackCuts=0x0;
312 } else if(fInputEvent->IsA()==AliAODEvent::Class()){
313 for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
314 AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
315 if(!curTrack->IsPrimaryCandidate()) continue;
316 if(abs(curTrack->Eta())>0.9) continue;
317 if(curTrack->Pt()<0.15) continue;
318 if(abs(curTrack->ZAtDCA())>2) continue;
319 fNumberOfESDTracks++;
320 }
321 }
72395bd9 322
51ead2dc 323 return fNumberOfESDTracks;
72395bd9 324}
325
51ead2dc 326//________________________________________________________________________
72395bd9 327Int_t AliAnalysisTaskResolution::CountTracks0914(){
51ead2dc 328 Int_t fNumberOfESDTracks = 0;
329 if(fInputEvent->IsA()==AliESDEvent::Class()){
330 // Using standard function for setting Cuts
331
332 AliStack *fMCStack = NULL;
63674d67 333 if (MCEvent()){
51ead2dc 334 fMCStack= fMCEvent->Stack();
335 if (!fMCStack) return 0;
336 }
72395bd9 337
51ead2dc 338 AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
339 EsdTrackCuts->SetMaxDCAToVertexZ(5);
340 EsdTrackCuts->SetEtaRange(0.9, 1.4);
341 EsdTrackCuts->SetPtRange(0.15);
342
343 for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
344 AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
345 if(!curTrack) continue;
346 if(EsdTrackCuts->AcceptTrack(curTrack) ){
347 if (fMCEvent){
344100c4 348 if(fMCStack && fEventCuts->GetSignalRejection() != 0){
349 Int_t isFromMBHeader = fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
51ead2dc 350 if( (isFromMBHeader < 1) ) continue;
351 }
352 }
353 fNumberOfESDTracks++;
354 }
355 }
356 EsdTrackCuts->SetEtaRange(-1.4, -0.9);
357 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
358 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
359 if(!curTrack) continue;
360 if(EsdTrackCuts->AcceptTrack(curTrack) ){
361 if (fMCEvent){
344100c4 362 if(fMCStack && fEventCuts->GetSignalRejection() != 0){
363 Int_t isFromMBHeader = fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
51ead2dc 364 if( (isFromMBHeader < 1) ) continue;
365 }
366 }
367 fNumberOfESDTracks++;
368 }
369 }
370 delete EsdTrackCuts;
371 EsdTrackCuts=0x0;
372 } else if(fInputEvent->IsA()==AliAODEvent::Class()){
373 for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
374 AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
375 if(abs(curTrack->Eta())<0.9 || abs(curTrack->Eta())>1.4 ) continue;
376 if(curTrack->Pt()<0.15) continue;
377 if(abs(curTrack->ZAtDCA())>5) continue;
378 fNumberOfESDTracks++;
379 }
380 }
381
382 return fNumberOfESDTracks;
72395bd9 383}
384
385//________________________________________________________________________
386void AliAnalysisTaskResolution::Terminate(Option_t *)
387{
63674d67 388
72395bd9 389}