Update from prabhat: Dpt-Dpt Corr
[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(),
39 fV0Reader(NULL),
40 fConversionGammas(NULL),
41 fConversionCuts(NULL),
42 fTreeEvent(NULL),
43 fTreeResolution(NULL),
44 fPrimVtxZ(0.),
45 fNContrVtx(0),
46 fNESDtracksEta09(0),
47 fNESDtracksEta0914(0),
48 fNESDtracksEta14(0),
31216736 49 fGammaRecCoords(5),
50 fGammaMCCoords(5),
72395bd9 51 fChi2ndf(0),
52 fIsHeavyIon(kFALSE),
53 fOutputList(NULL),
54 fEventList(NULL),
55 fResolutionList(NULL),
56 fESDEvent(NULL),
57 fMCEvent(NULL)
58{
59
60
61}
62
63
64//________________________________________________________________________
65AliAnalysisTaskResolution::AliAnalysisTaskResolution(const char *name) : AliAnalysisTaskSE(name),
66 fV0Reader(NULL),
67 fConversionGammas(NULL),
68 fConversionCuts(NULL),
69 fTreeEvent(NULL),
70 fTreeResolution(NULL),
71 fPrimVtxZ(0.),
72 fNContrVtx(0),
73 fNESDtracksEta09(0),
74 fNESDtracksEta0914(0),
75 fNESDtracksEta14(0),
31216736 76 fGammaRecCoords(5),
77 fGammaMCCoords(5),
72395bd9 78 fChi2ndf(0),
79 fIsHeavyIon(kFALSE),
80 fOutputList(NULL),
81 fEventList(NULL),
82 fResolutionList(NULL),
83 fESDEvent(NULL),
84 fMCEvent(NULL)
85{
86 // Default constructor
87
88 DefineInput(0, TChain::Class());
89 DefineOutput(1, TList::Class());
90}
91
92//________________________________________________________________________
93AliAnalysisTaskResolution::~AliAnalysisTaskResolution()
94{
95 // default deconstructor
96
97}
98//________________________________________________________________________
99void AliAnalysisTaskResolution::UserCreateOutputObjects()
100{
101 // Create User Output Objects
102
103 if(fOutputList != NULL){
104 delete fOutputList;
105 fOutputList = NULL;
106 }
107 if(fOutputList == NULL){
108 fOutputList = new TList();
109 fOutputList->SetOwner(kTRUE);
110 }
111
112 fEventList = new TList();
113 fEventList->SetName("EventList");
114 fEventList->SetOwner(kTRUE);
115 fOutputList->Add(fEventList);
116
117 fTreeEvent = new TTree("Event","Event");
118 fTreeEvent->Branch("primVtxZ",&fPrimVtxZ,"fPrimVtxZ/F");
119 fTreeEvent->Branch("nContrVtx",&fNContrVtx,"fNContrVtx/I");
120 fTreeEvent->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
121 fTreeEvent->Branch("nGoodTracksEta0914",&fNESDtracksEta0914,"fNESDtracksEta0914/I");
122 fTreeEvent->Branch("nGoodTracksEta14",&fNESDtracksEta14,"fNESDtracksEta14/I");
123 fEventList->Add(fTreeEvent);
124
125 fResolutionList = new TList();
126 fResolutionList->SetName("ResolutionList");
127 fResolutionList->SetOwner(kTRUE);
128 fOutputList->Add(fResolutionList);
72395bd9 129
130 fTreeResolution = new TTree("Resolution","Resolution");
31216736 131 fTreeResolution->Branch("RecCoords",&fGammaRecCoords);
132 fTreeResolution->Branch("MCCoords",&fGammaMCCoords);
72395bd9 133 fTreeResolution->Branch("chi2ndf",&fChi2ndf,"fChi2ndf/F");
134 fResolutionList->Add(fTreeResolution);
135 // V0 Reader Cuts
136 TString cutnumber = fConversionCuts->GetCutNumber();
137
138 PostData(1, fOutputList);
139}
140
141//________________________________________________________________________
142void AliAnalysisTaskResolution::UserExec(Option_t *){
143
51ead2dc 144 fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
72395bd9 145
51ead2dc 146 Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
147 if(eventQuality != 0){// Event Not Accepted
148 return;
149 }
150 fESDEvent = (AliESDEvent*) InputEvent();
151 if (fESDEvent==NULL) return;
72395bd9 152 if(MCEvent()){
153 fMCEvent = MCEvent();
154 }
51ead2dc 155
156 if(fMCEvent){
157 // Process MC Particle
158 if(fConversionCuts->GetSignalRejection() != 0){
159// if(fESDEvent->IsA()==AliESDEvent::Class()){
160 fConversionCuts->GetNotRejectedParticles( fConversionCuts->GetSignalRejection(),
161 fConversionCuts->GetAcceptedHeader(),
162 fMCEvent);
163// }
164// else if(fInputEvent->IsA()==AliAODEvent::Class()){
165// fConversionCuts->GetNotRejectedParticles( fConversionCuts->GetSignalRejection(),
166// fConversionCuts->GetAcceptedHeader(),
167// fInputEvent);
168// }
169 }
170 }
171
172
173 if(fIsHeavyIon && !fConversionCuts->IsCentralitySelected(fESDEvent)) return;
174 fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
175 fNESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
176 fNESDtracksEta14 = fNESDtracksEta09 + fNESDtracksEta0914;
177 if(fESDEvent){
178 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
179 fNContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
180 } else {
181 fNContrVtx = 0;
182 }
183// else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
184// if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
185// fNContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
186// } else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
187// fNContrVtx = 0;
188// }
189// }
190 }
191 fPrimVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
192
193 if (fIsHeavyIon){
194 if (!(fNESDtracksEta09 > 20 && fNESDtracksEta09 < 80)) return;
195 }
196
197
198 if (fTreeEvent){
199 fTreeEvent->Fill();
200 }
201
202 fConversionGammas=fV0Reader->GetReconstructedGammas();
203 ProcessPhotons();
204 PostData(1, fOutputList);
72395bd9 205}
206
207
208///________________________________________________________________________
209void AliAnalysisTaskResolution::ProcessPhotons(){
210
211 // Fill Histograms for QA and MC
212 for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast();firstGammaIndex++){
213 AliAODConversionPhoton *gamma=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(firstGammaIndex));
214 if (gamma ==NULL) continue;
215 if(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;
31216736 216 fGammaRecCoords(0) = gamma->GetPhotonPt();
217 fGammaRecCoords(1) = gamma->GetPhotonPhi();
218 fGammaRecCoords(2) = gamma->GetPhotonEta();
219 fGammaRecCoords(3) = gamma->GetConversionRadius();
220 fGammaRecCoords(4) = gamma->GetConversionZ();
72395bd9 221 fChi2ndf = gamma->GetChi2perNDF();
222 if(fMCEvent){
223// cout << "generating MC stack"<< endl;
224 AliStack *fMCStack = fMCEvent->Stack();
225 if (!fMCStack) continue;
226 TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
227 TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
228// cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
51ead2dc 229 if(fMCStack && fConversionCuts->GetSignalRejection() != 0){
230 Int_t isPosFromMBHeader
231 = fConversionCuts->IsParticleFromBGEvent(gamma->GetMCLabelPositive(), fMCStack, fESDEvent);
232 Int_t isNegFromMBHeader
233 = fConversionCuts->IsParticleFromBGEvent(gamma->GetMCLabelNegative(), fMCStack, fESDEvent);
234 if( (isNegFromMBHeader < 1) || (isPosFromMBHeader < 1)) continue;
235 }
72395bd9 236
237 if(posDaughter == NULL || negDaughter == NULL){
238 continue;
239 } else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
240 continue;
241 } else {
242// cout << "same mother" << endl;
243 Int_t pdgCodePos;
244 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
245 Int_t pdgCodeNeg;
246 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
247// cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
248 Int_t pdgCode;
249 if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
250// cout << "PDG code: " << pdgCode << endl;
251 if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
252 continue;
253 else if ( !(pdgCodeNeg==pdgCodePos)){
254 TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
255 if(pdgCode == 111)
256 continue;
257 else if (pdgCode == 221)
258 continue;
259 else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
260 if(pdgCode == 22){
31216736 261 fGammaMCCoords(0) = truePhotonCanditate->Pt();
262 fGammaMCCoords(1) = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
263 fGammaMCCoords(2) = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
264 fGammaMCCoords(3) = gamma->GetNegativeMCDaughter(fMCStack)->R();
265 fGammaMCCoords(4) = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
72395bd9 266
267 if (fTreeResolution){
268 fTreeResolution->Fill();
269 }
270 }
271 } else continue; //garbage
272 } else continue; //garbage
273 }
274 }
275 }
276}
277
278//________________________________________________________________________
279Int_t AliAnalysisTaskResolution::CountTracks09(){
51ead2dc 280 Int_t fNumberOfESDTracks = 0;
281 if(fInputEvent->IsA()==AliESDEvent::Class()){
282 // Using standard function for setting Cuts
283
284 AliStack *fMCStack = NULL;
285 if (fMCEvent){
286 fMCStack= fMCEvent->Stack();
287 if (!fMCStack) return 0;
288 }
289
290 Bool_t selectPrimaries=kTRUE;
291 AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
292 EsdTrackCuts->SetMaxDCAToVertexZ(2);
293 EsdTrackCuts->SetEtaRange(-0.9, 0.9);
294 EsdTrackCuts->SetPtRange(0.15);
295
296 for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
297 AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
298 if(!curTrack) continue;
299 if(EsdTrackCuts->AcceptTrack(curTrack) ){
300 if (fMCEvent){
301 if(fMCStack && fConversionCuts->GetSignalRejection() != 0){
302 Int_t isFromMBHeader
303 = fConversionCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
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;
333 if (fMCEvent){
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){
348 if(fMCStack && fConversionCuts->GetSignalRejection() != 0){
349 Int_t isFromMBHeader
350 = fConversionCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
351 if( (isFromMBHeader < 1) ) continue;
352 }
353 }
354 fNumberOfESDTracks++;
355 }
356 }
357 EsdTrackCuts->SetEtaRange(-1.4, -0.9);
358 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
359 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
360 if(!curTrack) continue;
361 if(EsdTrackCuts->AcceptTrack(curTrack) ){
362 if (fMCEvent){
363 if(fMCStack && fConversionCuts->GetSignalRejection() != 0){
364 Int_t isFromMBHeader
365 = fConversionCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
366 if( (isFromMBHeader < 1) ) continue;
367 }
368 }
369 fNumberOfESDTracks++;
370 }
371 }
372 delete EsdTrackCuts;
373 EsdTrackCuts=0x0;
374 } else if(fInputEvent->IsA()==AliAODEvent::Class()){
375 for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
376 AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
377 if(abs(curTrack->Eta())<0.9 || abs(curTrack->Eta())>1.4 ) continue;
378 if(curTrack->Pt()<0.15) continue;
379 if(abs(curTrack->ZAtDCA())>5) continue;
380 fNumberOfESDTracks++;
381 }
382 }
383
384 return fNumberOfESDTracks;
72395bd9 385}
386
387//________________________________________________________________________
388void AliAnalysisTaskResolution::Terminate(Option_t *)
389{
390// if (fStreamMaterial){
391// fStreamMaterial->GetFile()->Write();
392// }
393// if (fStreamResolution){
394// fStreamResolution->GetFile()->Write();
395// }
396}