-Move cent, z axis into sparse
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskMaterial.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 "AliAnalysisTaskMaterial.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(AliAnalysisTaskMaterial)
36
37 AliAnalysisTaskMaterial::AliAnalysisTaskMaterial() : AliAnalysisTaskSE(),
38    fV0Reader(NULL),
39    fConversionGammas(NULL),
40    fConversionCuts(NULL),
41    fOutputList(NULL),
42    fEventList(NULL),
43    fRecGammaList(NULL),
44    fAllMCGammaList(NULL),
45    fAllMCConvGammaList(NULL),
46    fTreeEvent(NULL),
47    fTreeMaterialRec(NULL),
48    fTreeMaterialAllGamma(NULL),
49    fTreeMaterialConvGamma(NULL),
50    fPrimVtxZ(0.),
51    fNContrVtx(0),
52    fNESDtracksEta09(0),
53    fNESDtracksEta0914(0),
54    fNESDtracksEta14(0),
55    fGammaMCPt(0.),
56    fGammaMCTheta(0.),
57    fGammaMCConvPt(0.),
58    fGammaMCConvTheta(0.),
59    fMCConvCords(5),
60    fMCConvDaughterProp(4),
61    fGammaPt(0.),
62    fGammaTheta(0.),
63    fGammaChi2NDF(0.),
64    fRecCords(5),
65    fDaughterProp(4),
66    fKind(0),      
67    fIsHeavyIon(kFALSE),
68    fIsMC(kFALSE),
69    fESDEvent(NULL),
70    fMCEvent(NULL)
71 {
72
73    
74 }
75
76
77 //________________________________________________________________________
78 AliAnalysisTaskMaterial::AliAnalysisTaskMaterial(const char *name) : AliAnalysisTaskSE(name),
79    fV0Reader(NULL),
80    fConversionGammas(NULL),
81    fConversionCuts(NULL),
82    fOutputList(NULL),
83    fEventList(NULL),
84    fRecGammaList(NULL),
85    fAllMCGammaList(NULL),
86    fAllMCConvGammaList(NULL),
87    fTreeEvent(NULL),
88    fTreeMaterialRec(NULL),
89    fTreeMaterialAllGamma(NULL),
90    fTreeMaterialConvGamma(NULL),
91    fPrimVtxZ(0.),
92    fNContrVtx(0),
93    fNESDtracksEta09(0),
94    fNESDtracksEta0914(0),
95    fNESDtracksEta14(0),
96    fGammaMCPt(0.),
97    fGammaMCTheta(0.),
98    fGammaMCConvPt(0.),
99    fGammaMCConvTheta(0.),
100    fMCConvCords(5),
101    fMCConvDaughterProp(4),
102    fGammaPt(0.),
103    fGammaTheta(0.),
104    fGammaChi2NDF(0.),
105    fRecCords(5),
106    fDaughterProp(4),
107    fKind(0),      
108    fIsHeavyIon(kFALSE),
109    fIsMC(kFALSE),
110    fESDEvent(NULL),
111    fMCEvent(NULL)
112 {
113    // Default constructor
114
115    
116    DefineInput(0, TChain::Class());
117    DefineOutput(1, TList::Class());
118 }
119
120 //________________________________________________________________________
121 AliAnalysisTaskMaterial::~AliAnalysisTaskMaterial()
122 {
123    // default deconstructor
124 }
125 //________________________________________________________________________
126 void AliAnalysisTaskMaterial::UserCreateOutputObjects()
127 {
128    // Create User Output Objects
129
130    if(fOutputList != NULL){
131       delete fOutputList;
132       fOutputList = NULL;
133    }
134    if(fOutputList == NULL){
135       fOutputList = new TList();
136       fOutputList->SetOwner(kTRUE);
137    }
138    
139    fEventList = new TList();
140    fEventList->SetName("EventList");
141    fEventList->SetOwner(kTRUE);
142    fOutputList->Add(fEventList);
143    
144    fTreeEvent = new TTree("Event","Event");   
145    fTreeEvent->Branch("primVtxZ",&fPrimVtxZ,"fPrimVtxZ/F");
146    fTreeEvent->Branch("nContrVtx",&fNContrVtx,"fNContrVtx/I");
147    fTreeEvent->Branch("nGoodTracksEta09",&fNESDtracksEta09,"fNESDtracksEta09/I");
148    fTreeEvent->Branch("nGoodTracksEta0914",&fNESDtracksEta0914,"fNESDtracksEta0914/I");
149    fTreeEvent->Branch("nGoodTracksEta14",&fNESDtracksEta14,"fNESDtracksEta14/I");
150    fEventList->Add(fTreeEvent);
151    
152    fRecGammaList= new TList();
153    fRecGammaList->SetName("RecGammaList");
154    fRecGammaList->SetOwner(kTRUE);
155    fOutputList->Add(fRecGammaList);
156    
157    
158    fTreeMaterialRec = new TTree("ConvPointRec","ConvPointRec");   
159    fTreeMaterialRec->Branch("recCords",&fRecCords);
160    fTreeMaterialRec->Branch("daughterProp",&fDaughterProp);
161    fTreeMaterialRec->Branch("pt",&fGammaPt,"fGammaPt/F");
162    fTreeMaterialRec->Branch("theta",&fGammaTheta,"fGammaTheta/F");
163    fTreeMaterialRec->Branch("chi2ndf",&fGammaChi2NDF,"fGammaChi2NDF/F");
164    if (fIsMC) {
165       fTreeMaterialRec->Branch("kind",&fKind,"fKind/b");
166    }   
167    fRecGammaList->Add(fTreeMaterialRec);
168    
169    if (fIsMC) {
170       fAllMCGammaList = new TList();
171       fAllMCGammaList->SetName("AllMCGammaList");
172       fAllMCGammaList->SetOwner(kTRUE);
173       fOutputList->Add(fAllMCGammaList);
174       
175       fTreeMaterialAllGamma = new TTree("AllGamma","AllGamma");   
176       fTreeMaterialAllGamma->Branch("pt",&fGammaMCPt,"fGammaMCPt/F");
177       fTreeMaterialAllGamma->Branch("theta",&fGammaMCTheta,"fGammaMCTheta/I");
178       fAllMCGammaList->Add(fTreeMaterialAllGamma);
179       
180       fAllMCConvGammaList = new TList();
181       fAllMCConvGammaList->SetName("AllMCGammaConvList");
182       fAllMCConvGammaList->SetOwner(kTRUE);
183       fOutputList->Add(fAllMCConvGammaList);
184
185 //       fMCConvCords = new Float_t[5];
186 //       fMCConvDaughterProp = new Float_t[4];
187
188       
189       fTreeMaterialConvGamma = new TTree("ConvGammaMC","ConvGammaMC");   
190       fTreeMaterialConvGamma->Branch("Cords",&fMCConvCords);
191       fTreeMaterialConvGamma->Branch("daughterProp",&fMCConvDaughterProp);
192       fTreeMaterialConvGamma->Branch("Pt",&fGammaMCConvPt,"fGammaMCConvPt/F");
193       fTreeMaterialConvGamma->Branch("Theta",&fGammaMCConvTheta,"fGammaMCConvTheta/F");   
194       fAllMCConvGammaList->Add(fTreeMaterialConvGamma);
195    }
196    
197    // V0 Reader Cuts
198    TString cutnumber = fConversionCuts->GetCutNumber();
199    PostData(1, fOutputList);
200    
201 }
202
203 //________________________________________________________________________
204 void AliAnalysisTaskMaterial::UserExec(Option_t *){
205
206    fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
207
208    Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
209         if(eventQuality != 0){// Event Not Accepted
210       return;
211    }
212    fESDEvent = (AliESDEvent*) InputEvent();
213    if (fESDEvent==NULL) return;
214    if(fIsHeavyIon && !fConversionCuts->IsCentralitySelected(fESDEvent)) return;
215         fNESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
216         fNESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
217         fNESDtracksEta14 = fNESDtracksEta09 + fNESDtracksEta0914;
218         if(fESDEvent){
219                 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
220                         fNContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
221                 } else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
222                         if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
223                                 fNContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
224                         } else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
225                                 fNContrVtx = 0;
226                         }
227                 }
228    }
229         fPrimVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
230         
231         if (fTreeEvent){
232       fTreeEvent->Fill();
233         }
234         
235    fConversionGammas=fV0Reader->GetReconstructedGammas();
236         if(MCEvent()){
237                 fMCEvent = MCEvent();
238         }
239    ProcessPhotons();
240         if(MCEvent()){
241                 ProcessMCPhotons();
242         }
243    PostData(1, fOutputList);
244 }
245
246 ///________________________________________________________________________
247 void AliAnalysisTaskMaterial::FillMCTree(Int_t stackPos){
248         AliStack *MCStack = fMCEvent->Stack();
249         TParticle* candidate = (TParticle *)MCStack->Particle(stackPos);
250         if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kFALSE)){
251                 fGammaMCPt = candidate->Pt();
252                 fGammaMCTheta = candidate->Theta();
253                 if (fTreeMaterialAllGamma){
254                         fTreeMaterialAllGamma->Fill();
255                 }       
256         }
257         if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kTRUE)){
258                 fGammaMCConvPt = candidate->Pt();
259                 fGammaMCConvTheta = candidate->Theta();
260                 TParticle* daughter1 = (TParticle *)MCStack->Particle(candidate->GetFirstDaughter()); 
261                 TParticle* daughter2 = (TParticle *)MCStack->Particle(candidate->GetLastDaughter()); 
262       fMCConvCords(0) = (Float_t)daughter1->Vx();
263       fMCConvCords(1) = (Float_t)daughter1->Vy();
264       fMCConvCords(2) = (Float_t)daughter1->Vz();
265       fMCConvCords(3) = (Float_t)daughter1->R();
266       fMCConvCords(4) = (Float_t)daughter1->Phi();
267       
268       fMCConvDaughterProp(0) = (Float_t)daughter1->Pt();
269       fMCConvDaughterProp(1) = (Float_t)daughter1->Theta();
270       fMCConvDaughterProp(2) = (Float_t)daughter2->Pt();
271       fMCConvDaughterProp(3) = (Float_t)daughter2->Theta();      
272                 
273                 if (fTreeMaterialConvGamma){
274                         fTreeMaterialConvGamma->Fill();
275                 }
276         } // Converted MC Gamma         
277 }
278
279 ///________________________________________________________________________
280 void AliAnalysisTaskMaterial::ProcessMCPhotons(){
281         // Loop over all primary MC particle
282         AliStack *ffMCStack = fMCEvent->Stack();
283         for(Int_t i = 0; i < ffMCStack->GetNprimary(); i++) {
284                 TParticle* particle = (TParticle *)ffMCStack->Particle(i);
285                 if (!particle) continue;                
286                 if (particle->GetPdgCode() == 111 && particle->GetFirstDaughter() >= ffMCStack->GetNprimary()){
287 //                      cout << "Undecayed pi0 found with mother: " << particle->GetMother(0) << endl;
288                         for (Int_t j = 0; j < 2 ; j++){
289                                 FillMCTree(particle->GetDaughter(j));
290                         }
291                 } else {
292                         FillMCTree(i);
293                 }       
294         }       
295 }
296
297 ///________________________________________________________________________
298 void AliAnalysisTaskMaterial::ProcessPhotons(){
299
300         // Fill Histograms for QA and MC
301    for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast();firstGammaIndex++){
302       AliAODConversionPhoton *gamma=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(firstGammaIndex));
303       if (gamma ==NULL) continue;
304       if(!fConversionCuts->PhotonIsSelected(gamma,fESDEvent)) continue;
305
306       fGammaPt = gamma->GetPhotonPt();
307                 fGammaTheta = gamma->GetPhotonTheta();
308                 fGammaChi2NDF = gamma->GetChi2perNDF();
309       fRecCords(0) = (Float_t)gamma->GetConversionX();
310       fRecCords(1) = (Float_t)gamma->GetConversionY();
311       fRecCords(2) = (Float_t)gamma->GetConversionZ();
312       fRecCords(3) = (Float_t)gamma->GetConversionRadius();
313       fRecCords(4) = (Float_t)gamma->GetPhotonPhi();
314       
315                 AliESDtrack * negTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelNegative());
316       AliESDtrack * posTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelPositive());
317       fDaughterProp(0) = (Float_t)posTrack->Pt();
318       fDaughterProp(1) = (Float_t)posTrack->Theta();
319       fDaughterProp(2) = (Float_t)negTrack->Pt();
320       fDaughterProp(3) = (Float_t)negTrack->Theta();
321
322                 fKind = 9;      
323                 
324                 if(fMCEvent){
325 //                      cout << "generating MC stack"<< endl;
326                         AliStack *fMCStack = fMCEvent->Stack();
327                         if (!fMCStack) continue;
328                         TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
329                         TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
330 //                      cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
331                         
332                         if(posDaughter == NULL || negDaughter == NULL){ 
333                                 fKind = 9; // garbage
334 //                              cout << "one of the daughters not available" << endl;
335                         } else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){ 
336                                 // Not Same Mother == Combinatorial Bck
337                                 fKind = 1;
338 //                              cout << "not the same mother" << endl;
339                                 Int_t pdgCodePos; 
340                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
341                                 Int_t pdgCodeNeg; 
342                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
343 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
344                                 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11)
345                                         fKind = 10; //Electron Combinatorial
346                                 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11 && (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1))
347                                         fKind = 15; //direct Electron Combinatorial
348                                 
349                                 if(TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==211)
350                                         fKind = 11; //Pion Combinatorial
351                                 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==2212) ||
352                                         (TMath::Abs(pdgCodePos)==2212 && TMath::Abs(pdgCodeNeg)==211))
353                                         fKind = 12; //Pion, Proton Combinatorics
354                                 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==11) ||
355                                         (TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==211))
356                                         fKind = 13; //Pion, Electron Combinatorics
357                                 if (TMath::Abs(pdgCodePos)==321 || TMath::Abs(pdgCodeNeg)==321) 
358                                         fKind = 14; //Kaon combinatorics
359
360                         } else {
361 //                              cout << "same mother" << endl;
362                                 Int_t pdgCodePos; 
363                                 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
364                                 Int_t pdgCodeNeg; 
365                                 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
366 //                              cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
367                                 Int_t pdgCode; 
368                                 if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
369 //                              cout << "PDG code: " << pdgCode << endl;
370                                 if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
371                                         fKind = 2; // combinatorics from hadronic decays
372                                 else if ( !(pdgCodeNeg==pdgCodePos)){
373                                         TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
374                                         Int_t motherLabelPhoton = truePhotonCanditate->GetMother(0);
375                                         if(pdgCode == 111) 
376                                                 fKind = 3; // pi0 Dalitz
377                                         else if (pdgCode == 221) 
378                                                 fKind = 4; // eta Dalitz
379                                         else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
380                                                 if(pdgCode == 22 && motherLabelPhoton < fMCStack->GetNprimary()){
381                                                         fKind = 0; // primary photons
382                                                 } else if (pdgCode == 22){
383                                                         fKind = 5; //secondary photons
384                                                 }               
385                                         } else  fKind = 9; //garbage
386                                 } else fKind = 9; //garbage
387                         }                                       
388                 }
389                 if (fTreeMaterialRec){
390          fTreeMaterialRec->Fill();
391       }
392         }
393 }
394
395 //________________________________________________________________________
396 Int_t AliAnalysisTaskMaterial::CountTracks09(){
397    Int_t fNumberOfESDTracks = 0;
398    if(fInputEvent->IsA()==AliESDEvent::Class()){
399    // Using standard function for setting Cuts
400       Bool_t selectPrimaries=kTRUE;
401       AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
402       EsdTrackCuts->SetMaxDCAToVertexZ(2);
403       EsdTrackCuts->SetEtaRange(-0.9, 0.9);
404       EsdTrackCuts->SetPtRange(0.15);
405       
406       for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
407          AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
408          if(!curTrack) continue;
409          if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
410       }
411       delete EsdTrackCuts;
412       EsdTrackCuts=0x0;
413    }
414    else if(fInputEvent->IsA()==AliAODEvent::Class()){
415       for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
416          AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
417          if(!curTrack->IsPrimaryCandidate()) continue;
418          if(abs(curTrack->Eta())>0.9) continue;
419          if(curTrack->Pt()<0.15) continue;
420          if(abs(curTrack->ZAtDCA())>2) continue;
421          fNumberOfESDTracks++;
422       }
423    }
424
425    return fNumberOfESDTracks;
426 }
427
428 Int_t AliAnalysisTaskMaterial::CountTracks0914(){
429
430     Int_t fNumberOfESDTracks = 0;
431     if(fInputEvent->IsA()==AliESDEvent::Class()){
432       // Using standard function for setting Cuts
433       AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
434       EsdTrackCuts->SetMaxDCAToVertexZ(5);
435       EsdTrackCuts->SetEtaRange(0.9, 1.4);
436       EsdTrackCuts->SetPtRange(0.15);
437       
438       for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
439          AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
440          if(!curTrack) continue;
441          if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
442       }
443       EsdTrackCuts->SetEtaRange(-1.4, -0.9);
444       for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
445          AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
446          if(!curTrack) continue;
447          if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
448       }
449       delete EsdTrackCuts;
450       EsdTrackCuts=0x0;
451    }
452    else if(fInputEvent->IsA()==AliAODEvent::Class()){
453       for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
454          AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
455          if(abs(curTrack->Eta())<0.9 || abs(curTrack->Eta())>1.4 ) continue;
456          if(curTrack->Pt()<0.15) continue;
457          if(abs(curTrack->ZAtDCA())>5) continue;
458          fNumberOfESDTracks++;
459       }
460    }
461    
462    return fNumberOfESDTracks;
463 }
464
465 //________________________________________________________________________
466 void AliAnalysisTaskMaterial::Terminate(Option_t *)
467 {
468 //    if (fStreamMaterial){
469 //       fStreamMaterial->GetFile()->Write();
470 //    }
471 }