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