TENDER becomes Tender, removing .so
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskResolution.cxx
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
31 class iostream;
32
33 using namespace std;
34
35 ClassImp(AliAnalysisTaskResolution)
36
37 //________________________________________________________________________
38 AliAnalysisTaskResolution::AliAnalysisTaskResolution() : AliAnalysisTaskSE(),
39         fV0Reader(NULL),
40         fConversionGammas(NULL),
41         fEventCuts(NULL),
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)
60 {
61
62    
63 }
64
65 //________________________________________________________________________
66 AliAnalysisTaskResolution::AliAnalysisTaskResolution(const char *name) : AliAnalysisTaskSE(name),
67         fV0Reader(NULL),
68         fConversionGammas(NULL),
69         fEventCuts(NULL),
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)
88 {
89         // Default constructor
90
91         DefineInput(0, TChain::Class());
92         DefineOutput(1, TList::Class());
93 }
94
95 //________________________________________________________________________
96 AliAnalysisTaskResolution::~AliAnalysisTaskResolution()
97 {
98         // default deconstructor
99    
100 }
101 //________________________________________________________________________
102 void AliAnalysisTaskResolution::UserCreateOutputObjects()
103 {
104         // Create User Output Objects
105
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");
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                                                                 
133                 fTreeResolution = new TTree("Resolution","Resolution");
134                 fTreeResolution->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
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);
141 }
142
143 //________________________________________________________________________
144 void AliAnalysisTaskResolution::UserExec(Option_t *){
145
146         fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
147
148         Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
149         if(eventQuality != 0){// Event Not Accepted
150                 return;
151         }
152         fESDEvent = (AliESDEvent*) InputEvent();
153         if (fESDEvent==NULL) return;
154         if(MCEvent()){
155                 fMCEvent = MCEvent();
156         }
157  
158         if(MCEvent()){
159                 // Process MC Particle
160                 if(fEventCuts->GetSignalRejection() != 0){
161 //              if(fESDEvent->IsA()==AliESDEvent::Class()){
162                         fEventCuts->GetNotRejectedParticles(    fEventCuts->GetSignalRejection(),
163                                                                                                         fEventCuts->GetAcceptedHeader(),
164                                                                                                         fMCEvent);
165 //              }
166 //              else if(fInputEvent->IsA()==AliAODEvent::Class()){
167 //                      fEventCuts->GetNotRejectedParticles(    fEventCuts->GetSignalRejection(),
168 //                                                                                                      fEventCuts->GetAcceptedHeader(),
169 //                                                                                                      fInputEvent);
170 //              }
171                 }
172         }
173
174    
175         if(fIsHeavyIon > 0 && !fEventCuts->IsCentralitySelected(fESDEvent)) return;
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         
195 //      if (fIsHeavyIon == 2){
196 //              if (!(fNESDtracksEta09 > 20 && fNESDtracksEta09 < 80)) return;
197 //      }       
198
199         
200         if (fTreeEvent){
201                 fTreeEvent->Fill();
202         }
203
204         fConversionGammas=fV0Reader->GetReconstructedGammas();
205         ProcessPhotons();
206         PostData(1, fOutputList);
207 }
208
209
210 ///________________________________________________________________________
211 void AliAnalysisTaskResolution::ProcessPhotons(){
212
213         // Fill Histograms for QA and MC
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;
218                 fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
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()){
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;
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);
235                                 if( (isNegFromMBHeader < 1) || (isPosFromMBHeader < 1)) continue;
236                         }
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){
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();
270                                                         }
271                                                 }               
272                                         } else continue; //garbage
273                                 } else continue; //garbage
274                         }
275                 }
276         }
277 }
278
279 //________________________________________________________________________
280 Int_t AliAnalysisTaskResolution::CountTracks09(){
281         Int_t fNumberOfESDTracks = 0;
282         if(fInputEvent->IsA()==AliESDEvent::Class()){
283         // Using standard function for setting Cuts
284                 
285                 AliStack *fMCStack = NULL;
286                 if (MCEvent()){
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){
302                                         if(fMCStack && fEventCuts->GetSignalRejection() != 0){
303                                                 Int_t isFromMBHeader = fEventCuts->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         }
322
323         return fNumberOfESDTracks;
324 }
325
326 //________________________________________________________________________
327 Int_t AliAnalysisTaskResolution::CountTracks0914(){
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 (MCEvent()){
334                         fMCStack= fMCEvent->Stack();
335                         if (!fMCStack) return 0;
336                 }       
337
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 && fEventCuts->GetSignalRejection() != 0){
349                                                 Int_t isFromMBHeader = fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
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){
362                                         if(fMCStack && fEventCuts->GetSignalRejection() != 0){
363                                                 Int_t isFromMBHeader = fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack, fESDEvent);
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;
383 }
384
385 //________________________________________________________________________
386 void AliAnalysisTaskResolution::Terminate(Option_t *)
387 {
388         
389 }