1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Authors: Friederike Bock *
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 **************************************************************************/
16 ////////////////////////////////////////////////
17 //---------------------------------------------
18 // QA Task for V0 Reader V1
19 //---------------------------------------------
20 ////////////////////////////////////////////////
22 #include "AliAnalysisTaskMaterial.h"
24 #include "AliAnalysisManager.h"
25 #include "TParticle.h"
27 #include "AliPIDResponse.h"
28 #include "AliESDtrackCuts.h"
35 ClassImp(AliAnalysisTaskMaterial)
37 //________________________________________________________________________
38 AliAnalysisTaskMaterial::AliAnalysisTaskMaterial(const char *name) : AliAnalysisTaskSE(name),
40 fConversionGammas(NULL),
41 fConversionCuts(NULL),
42 fStreamMaterial(NULL),
43 fStreamResolution(NULL),
49 // Default constructor
51 DefineInput(0, TChain::Class());
52 DefineOutput(1, TList::Class());
55 //________________________________________________________________________
56 AliAnalysisTaskMaterial::~AliAnalysisTaskMaterial()
58 // default deconstructor
60 delete fStreamMaterial;
61 fStreamMaterial = 0x0;
63 if(fStreamResolution){
64 delete fStreamResolution;
65 fStreamResolution = 0x0;
68 //________________________________________________________________________
69 void AliAnalysisTaskMaterial::UserCreateOutputObjects()
71 // Create User Output Objects
73 if(fOutputList != NULL){
77 if(fOutputList == NULL){
78 fOutputList = new TList();
79 fOutputList->SetOwner(kTRUE);
83 TString cutnumber = fConversionCuts->GetCutNumber();
85 fStreamMaterial = new TTreeSRedirector(Form("GammaConvV1_Material_%s.root",cutnumber.Data()),"recreate");
86 fStreamResolution = new TTreeSRedirector(Form("GammaConvV1_Resolution_%s.root",cutnumber.Data()),"recreate");
87 PostData(1, fOutputList);
90 //________________________________________________________________________
91 void AliAnalysisTaskMaterial::UserExec(Option_t *){
93 fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
95 Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
96 if(eventQuality != 0){// Event Not Accepted
99 fESDEvent = (AliESDEvent*) InputEvent();
100 if (fESDEvent==NULL) return;
101 if(fIsHeavyIon && !fConversionCuts->IsCentralitySelected(fESDEvent)) return;
102 Int_t nESDtracksEta09 = CountTracks09(); // Estimate Event Multiplicity
103 Int_t nESDtracksEta0914 = CountTracks0914(); // Estimate Event Multiplicity
104 Int_t nESDtracksEta14; // Estimate Event Multiplicity
105 nESDtracksEta14= nESDtracksEta09 + nESDtracksEta0914;
108 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
109 nContrVtx = fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
110 } else if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
111 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
112 nContrVtx = fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
114 } else if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
119 Float_t primVtxZ = fESDEvent->GetPrimaryVertex()->GetZ();
121 if (fStreamMaterial){
122 (*fStreamMaterial)<<"Event"
123 << "primVtxZ=" << primVtxZ
124 << "nContrVtx=" << nContrVtx
125 << "nGoodTracksEta09=" << nESDtracksEta09
126 << "nGoodTracksEta0914=" << nESDtracksEta0914
127 << "nGoodTracksEta14=" << nESDtracksEta14
131 fConversionGammas=fV0Reader->GetReconstructedGammas();
133 fMCEvent = MCEvent();
139 PostData(1, fOutputList);
142 ///________________________________________________________________________
143 void AliAnalysisTaskMaterial::FillMCTree(Int_t stackPos){
144 AliStack *MCStack = fMCEvent->Stack();
145 TParticle* candidate = (TParticle *)MCStack->Particle(stackPos);
146 if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kFALSE)){
147 Float_t gammaPt = candidate->Pt();
148 Float_t gammaTheta = candidate->Theta();
149 if (fStreamMaterial){
150 (*fStreamMaterial)<<"AllGamma"
153 << "theta=" << gammaTheta
157 if(fConversionCuts->PhotonIsSelectedMC(candidate,MCStack,kTRUE)){
158 Float_t gammaPt = candidate->Pt();
159 Float_t gammaTheta = candidate->Theta();
160 TParticle* daughter1 = (TParticle *)MCStack->Particle(candidate->GetFirstDaughter());
161 TParticle* daughter2 = (TParticle *)MCStack->Particle(candidate->GetLastDaughter());
163 coord(0) = daughter1->Vx();
164 coord(1) = daughter1->Vy();
165 coord(2) = daughter1->Vz();
166 coord(3) = daughter1->R();
167 coord(4) = candidate->Phi();
168 TVectorF daughterProp(4);
169 if (daughter1-> GetPdgCode() < 0){
170 daughterProp(0) = daughter2->Pt();
171 daughterProp(1) = daughter2->Theta();
172 daughterProp(2) = daughter1->Pt();
173 daughterProp(3) = daughter1->Theta();
175 daughterProp(0) = daughter1->Pt();
176 daughterProp(1) = daughter1->Theta();
177 daughterProp(2) = daughter2->Pt();
178 daughterProp(3) = daughter2->Theta();
180 if (fStreamMaterial){
181 (*fStreamMaterial)<<"ConvGammaMC"
184 << "theta=" << gammaTheta
185 << "coord.=" << &coord
186 << "daughterProp.=" << &daughterProp
189 } // Converted MC Gamma
192 ///________________________________________________________________________
193 void AliAnalysisTaskMaterial::ProcessMCPhotons(){
194 // Loop over all primary MC particle
195 AliStack *ffMCStack = fMCEvent->Stack();
196 for(Int_t i = 0; i < ffMCStack->GetNprimary(); i++) {
197 TParticle* particle = (TParticle *)ffMCStack->Particle(i);
198 if (!particle) continue;
199 if (particle->GetPdgCode() == 111 && particle->GetFirstDaughter() >= ffMCStack->GetNprimary()){
200 // cout << "Undecayed pi0 found with mother: " << particle->GetMother(0) << endl;
201 for (Int_t j = 0; j < 2 ; j++){
202 FillMCTree(particle->GetDaughter(j));
210 ///________________________________________________________________________
211 void AliAnalysisTaskMaterial::ProcessPhotons(){
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 // cout << "i= " <<firstGammaIndex << " of "<< fConversionGammas->GetEntriesFast() << endl;
219 Float_t gammaPt = gamma->GetPhotonPt();
220 // Float_t gammaP = gamma->GetPhotonP();
221 Float_t gammaTheta = gamma->GetPhotonTheta();
222 Float_t gammaEta = gamma->GetPhotonEta();
223 Float_t gammaChi2NDF = gamma->GetChi2perNDF();
224 // Float_t gammaX = gamma->GetConversionX();
225 // Float_t gammaY = gamma->GetConversionY();
226 Float_t gammaZ = gamma->GetConversionZ();
227 Float_t gammaR = gamma->GetConversionRadius();
228 Float_t gammaPhi = gamma->GetPhotonPhi();
231 coord(0) = gamma->GetConversionX();
232 coord(1) = gamma->GetConversionY();
233 coord(2) = gamma->GetConversionZ();
234 coord(3) = gamma->GetConversionRadius();
235 coord(4) = gamma->GetPhotonPhi();
236 TVectorF daughterProp(4);
237 AliESDtrack * negTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelNegative());
238 AliESDtrack * posTrack = fConversionCuts->GetESDTrack(fESDEvent, gamma->GetTrackLabelPositive());
239 daughterProp(0) = posTrack->Pt();
240 daughterProp(1) = posTrack->Theta();
241 daughterProp(2) = negTrack->Pt();
242 daughterProp(3) = negTrack->Theta();
246 // cout << "generating MC stack"<< endl;
247 AliStack *fMCStack = fMCEvent->Stack();
248 if (!fMCStack) continue;
249 TParticle *posDaughter = gamma->GetPositiveMCDaughter(fMCStack);
250 TParticle *negDaughter = gamma->GetNegativeMCDaughter(fMCStack);
251 // cout << "generate Daughters: "<<posDaughter << "\t" << negDaughter << endl;
253 if(posDaughter == NULL || negDaughter == NULL){
255 // cout << "one of the daughters not available" << endl;
256 } else if(posDaughter->GetMother(0) != negDaughter->GetMother(0) || (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1)){
257 // Not Same Mother == Combinatorial Bck
259 // cout << "not the same mother" << endl;
261 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
263 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
264 // cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
265 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11)
266 kind = 10; //Electron Combinatorial
267 if(TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==11 && (posDaughter->GetMother(0) == negDaughter->GetMother(0) && posDaughter->GetMother(0) ==-1))
268 kind = 15; //direct Electron Combinatorial
270 if(TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==211)
271 kind = 11; //Pion Combinatorial
272 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==2212) ||
273 (TMath::Abs(pdgCodePos)==2212 && TMath::Abs(pdgCodeNeg)==211))
274 kind = 12; //Pion, Proton Combinatorics
275 if((TMath::Abs(pdgCodePos)==211 && TMath::Abs(pdgCodeNeg)==11) ||
276 (TMath::Abs(pdgCodePos)==11 && TMath::Abs(pdgCodeNeg)==211))
277 kind = 13; //Pion, Electron Combinatorics
278 if (TMath::Abs(pdgCodePos)==321 || TMath::Abs(pdgCodeNeg)==321)
279 kind = 14; //Kaon combinatorics
282 // cout << "same mother" << endl;
284 if (posDaughter->GetPdgCode()) pdgCodePos = posDaughter->GetPdgCode(); else continue;
286 if (negDaughter->GetPdgCode()) pdgCodeNeg = negDaughter->GetPdgCode(); else continue;
287 // cout << "PDG codes daughters: " << pdgCodePos << "\t" << pdgCodeNeg << endl;
289 if (gamma->GetMCParticle(fMCStack)->GetPdgCode()) pdgCode = gamma->GetMCParticle(fMCStack)->GetPdgCode(); else continue;
290 // cout << "PDG code: " << pdgCode << endl;
291 if(TMath::Abs(pdgCodePos)!=11 || TMath::Abs(pdgCodeNeg)!=11)
292 kind = 2; // combinatorics from hadronic decays
293 else if ( !(pdgCodeNeg==pdgCodePos)){
294 TParticle *truePhotonCanditate = gamma->GetMCParticle(fMCStack);
295 Int_t motherLabelPhoton = truePhotonCanditate->GetMother(0);
297 kind = 3; // pi0 Dalitz
298 else if (pdgCode == 221)
299 kind = 4; // eta Dalitz
300 else if (!(negDaughter->GetUniqueID() != 5 || posDaughter->GetUniqueID() !=5)){
301 if(pdgCode == 22 && motherLabelPhoton < fMCStack->GetNprimary()){
302 kind = 0; // primary photons
303 } else if (pdgCode == 22){
304 kind = 5; //secondary photons
307 Float_t mcPt = truePhotonCanditate->Pt();
308 Float_t mcR = gamma->GetNegativeMCDaughter(fMCStack)->R();
309 Float_t mcZ = gamma->GetNegativeMCDaughter(fMCStack)->Vz();
310 Float_t mcPhi = gamma->GetNegativeMCDaughter(fMCStack)->Phi();
311 Float_t mcEta = gamma->GetNegativeMCDaughter(fMCStack)->Eta();
312 if (fStreamResolution){
313 (*fStreamResolution)<<"Resolution"
314 << "ESDpt=" << gammaPt
315 << "ESDphi=" << gammaPhi
316 << "ESDeta=" << gammaEta
324 << "chi2ndf=" << gammaChi2NDF
328 } else kind = 9; //garbage
329 } else kind = 9; //garbage
331 // cout << gammaPt << "\t" << gammaP<< "\t" << gammaEta<< "\t" <<gammaChi2NDF << "\t" << gammaX<< "\t" <<gammaY << "\t" << gammaZ<< "\t" << gammaR<< "\t" << gammaPhi<< "\t" <<kind << endl;
333 if (fStreamMaterial){
334 (*fStreamMaterial)<<"ConvPointRec"
336 << "theta=" << gammaTheta
337 << "chi2ndf=" << gammaChi2NDF
339 << "coord.=" << &coord
340 << "daughterProp.=" << &daughterProp
344 if (fStreamMaterial){
345 (*fStreamMaterial)<<"ConvPointRec"
347 << "theta=" << gammaTheta
348 << "chi2ndf=" << gammaChi2NDF
349 << "coord.=" << &coord
350 << "daughterProp.=" << &daughterProp
357 //________________________________________________________________________
358 Int_t AliAnalysisTaskMaterial::CountTracks09(){
359 Int_t fNumberOfESDTracks = 0;
360 if(fInputEvent->IsA()==AliESDEvent::Class()){
361 // Using standard function for setting Cuts
362 Bool_t selectPrimaries=kTRUE;
363 AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
364 EsdTrackCuts->SetMaxDCAToVertexZ(2);
365 EsdTrackCuts->SetEtaRange(-0.9, 0.9);
366 EsdTrackCuts->SetPtRange(0.15);
368 for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
369 AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
370 if(!curTrack) continue;
371 // if(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
372 // if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) continue;
374 if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
379 else if(fInputEvent->IsA()==AliAODEvent::Class()){
380 for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
381 AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
382 if(!curTrack->IsPrimaryCandidate()) continue;
383 if(abs(curTrack->Eta())>0.9) continue;
384 if(curTrack->Pt()<0.15) continue;
385 if(abs(curTrack->ZAtDCA())>2) continue;
386 fNumberOfESDTracks++;
390 return fNumberOfESDTracks;
393 Int_t AliAnalysisTaskMaterial::CountTracks0914(){
395 // Using standard function for setting Cuts ; We use TPCOnlyTracks for outer eta region
396 //Bool_t selectPrimaries=kTRUE;
397 // EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
398 Int_t fNumberOfESDTracks = 0;
399 if(fInputEvent->IsA()==AliESDEvent::Class()){
400 // Using standard function for setting Cuts
401 AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
402 EsdTrackCuts->SetMaxDCAToVertexZ(5);
403 EsdTrackCuts->SetEtaRange(0.9, 1.4);
404 EsdTrackCuts->SetPtRange(0.15);
406 for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
407 AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
408 if(!curTrack) continue;
409 // if(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
410 // if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(abs(curTrack->GetLabel()), fMCStack)) continue;
412 if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
414 EsdTrackCuts->SetEtaRange(-1.4, -0.9);
415 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
416 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
417 if(!curTrack) continue;
418 if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
423 else if(fInputEvent->IsA()==AliAODEvent::Class()){
424 for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
425 AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
426 // if(!curTrack->IsPrimaryCandidate()) continue;
427 if(abs(curTrack->Eta())<0.9 || abs(curTrack->Eta())>1.4 ) continue;
428 if(curTrack->Pt()<0.15) continue;
429 if(abs(curTrack->ZAtDCA())>5) continue;
430 fNumberOfESDTracks++;
434 return fNumberOfESDTracks;
437 //________________________________________________________________________
438 void AliAnalysisTaskMaterial::Terminate(Option_t *)
440 // if (fStreamMaterial){
441 // fStreamMaterial->GetFile()->Write();
443 // if (fStreamResolution){
444 // fStreamResolution->GetFile()->Write();