]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGGA/GammaConv/AliAnalysisTaskdPhi.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
adding track cut filter. Only add used ME tracks to ME
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskdPhi.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Authors: Svein Lindal *
5 * Version 1.0 *
6 * *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
17////////////////////////////////////////////////
18//---------------------------------------------
19// Class doing conversion gamma dPhi correlations
20// Gamma Conversion analysis
21//---------------------------------------------
22////////////////////////////////////////////////
23
24#include "AliAnalysisTaskdPhi.h"
25
26#include <TH2I.h>
27#include <TList.h>
28#include <TChain.h>
29
30#include <AliAnalysisManager.h>
31#include <AliInputEventHandler.h>
32#include <AliESDInputHandler.h>
33#include <AliAODInputHandler.h>
34#include <AliAnalysisFilter.h>
35
36#include "AliConversionTrackCuts.h"
37#include "AliConversionCuts.h"
38#include "AliAODConversionPhoton.h"
39#include "AliAODConversionMother.h"
40#include "AliAnaConvCorrPhoton.h"
41#include "AliAnaConvCorrPion.h"
42#include "AliAnaConvIsolation.h"
43// Author Svein Lindal <slindal@fys.uio.no>
44using namespace std;
45
46ClassImp(AliAnalysisTaskdPhi)
47
48
49//________________________________________________________________________
50AliAnalysisTaskdPhi::AliAnalysisTaskdPhi(const char *name) : AliAnalysisTaskSE(name),
51 fHistograms(NULL),
52 fHistoGamma(NULL),
53 fHistoPion(NULL),
54 fV0Filter(NULL),
55 fTrackCuts(),
56 fGammas(NULL),
57 fPions(NULL),
58 hMETracks(NULL),
59 hMEPhotons(NULL),
60 hMEPions(NULL),
61 hMEvents(NULL),
62 fPhotonCorr(NULL),
63 fPionCorr(NULL),
64 fIsoAna(NULL),
65 fL1(-1),
66 fL2(-1),
67 fDeltaAODBranchName("AliAODGammaConversion_gamma"),
68 fAxistPt(),
69 fAxiscPt(),
70 fAxisEta(),
71 fAxisPhi(),
72 fAxisCent(),
73 fAxisZ(),
74 fAxisPiM()
75{
76 //constructor
77 fAxistPt.SetNameTitle("tPtAxis", "tPt");
78 fAxistPt.Set(20, 0, 100);
79
80 fAxiscPt.SetNameTitle("cPtAxis", "cPt");
81 fAxiscPt.Set(20, 0, 100);
82
83 fAxisEta.SetNameTitle("EtaAxis", "Eta");
84 fAxisEta.Set(160, -0.8, 0.8);
85
86 fAxisPhi.SetNameTitle("PhiAxis", "Phi");
87 fAxisPhi.Set(128, 0, TMath::TwoPi());
88
89 fAxisZ.SetNameTitle("ZAxis", "Z");
90 fAxisZ.Set(4, -10, 10);
91 fAxisCent.SetNameTitle("CentAxis", "Cent");
92
93 Double_t centbins[5] = {0, 10, 30, 60, 100.1};
94 fAxisCent.Set(4, centbins);
95
96 Double_t mbins[7] = {0.1, 0.11, 0.12, 0.15, 0.16, 0.18, 0.2};
97 fAxisPiM.SetNameTitle("InvMassPi0", "Invariant mass");
98 fAxisPiM.Set(6, mbins);
99
100 fGammas = new TObjArray();
101 fGammas->SetOwner(kFALSE);
102
103 fPions = new TObjArray();
104 fPions->SetOwner(kFALSE);
105
106 // Define input and output slots here
107 DefineInput(0, TChain::Class());
108 //DefineInput(1, TClonesArray::Class());
109 DefineOutput(1, TList::Class());
110 DefineOutput(2, TList::Class());
111 DefineOutput(3, TList::Class());
112}
113
114
115
116//________________________________________________________________________
117AliAnalysisTaskdPhi::~AliAnalysisTaskdPhi(){
118 //destructor
119 if(fPions)
120 delete fPions;
121 fPions = NULL;
122
123 if(fGammas)
124 delete fGammas;
125 fGammas = NULL;
126
127 if(fIsoAna)
128 delete fIsoAna;
129 fIsoAna = NULL;
130
131 if(fV0Filter)
132 delete fV0Filter;
133 fV0Filter = NULL;
134
135 if(fHistograms)
136 delete fHistograms;
137 fHistograms = NULL;
138
139if(fHistoPion)
140 delete fHistoPion;
141 fHistoPion = NULL;
142
143 if(fHistoGamma)
144 delete fHistoGamma;
145 fHistoGamma = NULL;
146
147}
148
149///________________________________________________________________________
150void AliAnalysisTaskdPhi::SetUpCorrObjects() {
151 ///Creat corr obj
152 fIsoAna = new AliAnaConvIsolation();
153
154
155 fPhotonCorr = new TObjArray();
156 fPionCorr = new TObjArray();
157
158 TList * hPhoton = new TList();
159 hPhoton->SetName("hPhotonCorr");
160 hPhoton->SetOwner(kTRUE);
161 fHistoGamma->Add(hPhoton);
162
163 TList * hPion = new TList();
164 hPion->SetName("hPionCorr");
165 hPion->SetOwner(kTRUE);
166 fHistoPion->Add(hPion);
167
168
169 for(Int_t iz = 0; iz < fAxisZ.GetNbins(); iz++) {
170 TObjArray * photonArray = new TObjArray();
171 photonArray->SetOwner(kTRUE);
172 fPhotonCorr->AddAt(photonArray, iz);
173
174 TObjArray * pionArray = new TObjArray();
175 pionArray->SetOwner(kTRUE);
176 fPionCorr->AddAt(pionArray, iz);
177
178 TList * photonList = new TList();
179 photonList->SetName(Form("photon_%d", iz));
180 photonList->SetOwner(kTRUE);
181 hPhoton->AddAt(photonList, iz);
182
183 TList * pionList = new TList();
184 pionList->SetName(Form("pion_%d", iz));
185 pionList->SetOwner(kTRUE);
186 hPion->AddAt(pionList, iz);
187
188 for(Int_t ic = 0; ic < fAxisCent.GetNbins(); ic++) {
189
190 TString nameString = Form("%d_%d", iz, ic);
191 TString titleString = Form("%f < Z < %f ... %f cent %f",
192 fAxisZ.GetBinLowEdge(iz+1), fAxisZ.GetBinUpEdge(iz+1),
193 fAxisCent.GetBinLowEdge(ic+1), fAxisCent.GetBinUpEdge(ic+1));
194
195
196
197 AliAnaConvCorrPhoton * photonCorr = new AliAnaConvCorrPhoton(Form("PhotonCorr_%s", nameString.Data()), Form("photon %s", titleString.Data()));
198 photonArray->AddAt(photonCorr, ic);
199 photonCorr->GetAxistPt().Set(fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
200 photonCorr->GetAxiscPt().Set(fAxiscPt.GetNbins(), fAxiscPt.GetXbins()->GetArray());
201 photonCorr->CreateHistograms();
202 photonList->Add(photonCorr->GetHistograms());
203
204 AliAnaConvCorrPion * pionCorr = new AliAnaConvCorrPion(Form("PionCorr_%s", nameString.Data()), Form("pion %s", titleString.Data()));
205 pionArray->AddAt(pionCorr, ic);
206 pionCorr->GetAxistPt().Set(fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
207 pionCorr->GetAxiscPt().Set(fAxiscPt.GetNbins(), fAxiscPt.GetXbins()->GetArray());
208 pionCorr->GetAxisM().Set(fAxisPiM.GetNbins(), fAxisPiM.GetXbins()->GetArray());
209 pionCorr->CreateHistograms();
210 pionList->Add(pionCorr->GetHistograms());
211 }
212 }
213}
214
215//________________________________________________________________________
216void AliAnalysisTaskdPhi::UserCreateOutputObjects() {
217 // Create histograms
218
219 fHistograms = new TList();
220 fHistograms->SetName("dPhi_histograms");
221 fHistograms->SetOwner(kTRUE);
222
223 fHistoGamma = new TList();
224 fHistoGamma->SetName("Gamma_histo");
225 fHistoGamma->SetOwner(kTRUE);
226
227 fHistoPion = new TList();
228 fHistoPion->SetName("Pion_histo");
229 fHistoPion->SetOwner(kTRUE);
230
231
232 if(fV0Filter) {
233 fV0Filter->InitCutHistograms();
234 fHistograms->Add(fV0Filter->GetCutHistograms());
235 }
236
237
238 SetUpCorrObjects();
239
240
241 ///Set up ME histograms
242 TList * MEHistograms = new TList();
243 MEHistograms->SetName("MEHistograms");
244 MEHistograms->SetOwner(kTRUE);
245 fHistograms->Add(MEHistograms);
246
247 hMETracks = new TObjArray();
248 hMETracks->SetName("TrackArray");
249 hMETracks->SetOwner(kTRUE);
250 hMEPhotons = new TObjArray();
251 hMEPhotons->SetName("PhotonArray");
252 hMEPhotons->SetOwner(kTRUE);
253 hMEPions = new TObjArray();
254 hMEPions->SetName("PionArray");
255 hMEPions->SetOwner(kTRUE);
256
257 MEHistograms->Add(hMETracks);
258 MEHistograms->Add(hMEPions);
259 MEHistograms->Add(hMEPhotons);
260
261 hMEvents = new TH2I("hMEvents", "Nevents vs centrality vertexz",
262 fAxisZ.GetNbins(), fAxisZ.GetBinLowEdge(1), fAxisZ.GetBinUpEdge(fAxisZ.GetNbins()),
263 fAxisCent.GetNbins(), fAxisCent.GetBinLowEdge(1), fAxisCent.GetBinUpEdge(fAxisCent.GetNbins()));
264 hMEvents->GetYaxis()->Set(fAxisCent.GetNbins(), fAxisCent.GetXbins()->GetArray());
265 MEHistograms->Add(hMEvents);
266
267
268 TList axesList;
269 axesList.AddAt(&GetAxisEta(), 0);
270 axesList.AddAt(&GetAxisPhi(), 1);
271 axesList.AddAt(&GetAxistPt(), 2);
272 axesList.SetOwner(kFALSE);
273
274 TList piAxesList;
275 piAxesList.AddAt(&GetAxisEta(), 0);
276 piAxesList.AddAt(&GetAxisPhi(), 1);
277 piAxesList.AddAt(&GetAxistPt(), 2);
278 piAxesList.AddAt(&GetAxisPiMass(), 3);
279 piAxesList.SetOwner(kFALSE);
280
281
282 TList * outAxesList = new TList();
283 outAxesList->Add(&fAxisCent);
284 outAxesList->Add(&fAxisZ);
285 fHistograms->Add(outAxesList);
286
287 // for(Int_t iz = 0; iz < fAxisZ.GetNbins(); iz++) {
288 // TObjArray * trackArray = new TObjArray();
289 // trackArray->SetName(Form("METracks_%d", iz));
290 // trackArray->SetOwner(kTRUE);
291 // TObjArray * photonArray = new TObjArray();
292 // photonArray->SetName(Form("MEPhotons_%d", iz));
293 // photonArray->SetOwner(kTRUE);
294 // TObjArray * pionArray = new TObjArray();
295 // pionArray->SetName(Form("MEPions_%d", iz));
296 // pionArray->SetOwner(kTRUE);
297
298
299 // hMEPions->AddAt(pionArray, iz);
300 // hMETracks->AddAt(trackArray, iz);
301 // hMEPhotons->AddAt(photonArray, iz);
302
303 // for(Int_t ic = 0; ic < fAxisCent.GetNbins(); ic++) {
304
305 // TString nameString = Form("%d_%d", iz, ic);
306 // TString titleString = Form("%f < Z < %f ... %f cent %f",
307 // fAxisZ.GetBinLowEdge(iz+1), fAxisZ.GetBinUpEdge(iz+1),
308 // fAxisCent.GetBinLowEdge(ic+1), fAxisCent.GetBinUpEdge(ic+1));
309
310
311 // THnSparseF * trackHistogram = CreateSparse(Form("tracks_%s", nameString.Data()),
312 // Form("tracks %s", titleString.Data()), &axesList );
313 // trackArray->AddAt(trackHistogram, ic);
314
315 // THnSparseF * photonHistogram = CreateSparse(Form("photons_%s", nameString.Data()),
316 // Form("photons %s", titleString.Data()), &axesList );
317 // photonArray->AddAt(photonHistogram, ic);
318
319 // THnSparseF * pionHistogram = CreateSparse(Form("pions_%s", nameString.Data()),
320 // Form("pions %s", titleString.Data()), &piAxesList );
321 // pionArray->AddAt(pionHistogram, ic);
322 // }
323 // }
324
325 PostData(1, fHistograms);
326 PostData(2, fHistoGamma);
327 PostData(3, fHistoPion);
328
329}
330
331///________________________________________________________________________
332THnSparseF * AliAnalysisTaskdPhi::CreateSparse(TString nameString, TString titleString, TList * axesList) {
333 ///Creat sparse
334 const Int_t dim = axesList->GetSize();
335
336 TAxis * axes[dim];
337 Int_t bins[dim];
338 Double_t min[dim];
339 Double_t max[dim];
340
341 for(Int_t i = 0; i<dim; i++) {
342 TAxis * axis = dynamic_cast<TAxis*>(axesList->At(i));
343 if(axis) {
344 axes[i] = axis;
345 } else {
346 cout << "AliAnalysisTaskdPhi::CreateSparse: Error error, all the axes are not present in axis list" << endl;
347 return NULL;
348 }
349 }
350
351 for(Int_t i = 0; i<dim; i++) {
352 bins[i] = axes[i]->GetNbins();
353 min[i] = axes[i]->GetBinLowEdge(1);
354 max[i] = axes[i]->GetBinUpEdge(axes[i]->GetNbins());
355 }
356
357 THnSparseF * sparse = new THnSparseF(Form("METracks_%s", nameString.Data()),
358 Form("tracks %s", titleString.Data()),
359 dim, bins, min, max);
360
361 for(Int_t i = 0; i<dim; i++) {
362 sparse->GetAxis(i)->SetNameTitle(axes[i]->GetName(), axes[i]->GetTitle() );
363 if(axes[i]->GetXbins()->GetSize() > 0) {
364 sparse->SetBinEdges(i, axes[i]->GetXbins()->GetArray() );
365 }
366 }
367
368 return sparse;
369}
370
371//________________________________________________________________________
372void AliAnalysisTaskdPhi::UserExec(Option_t *) {
373 ///User exec.
374
375 //if(! fV0Filter->EventIsSelected(fInputEvent)) return;
376
377
378 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
379 Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
380
381 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
382 if (!inputHandler) {
383 cout << "cout no input event handler"<<endl;
384 return;
385 }
386
387
388 if ( fV0Filter && !fV0Filter->GetPIDResponse() ) {
389 if ( inputHandler->GetPIDResponse() ){
390 fV0Filter->SetPIDResponse( inputHandler->GetPIDResponse() );
391 } else {
392
393 //AOD case
394 if (isAOD){
395 if (!fV0Filter->GetPIDResponse()){
396 fV0Filter->InitAODpidUtil(1);
397 }
398 }
399 }
400 }
401
402 Double_t centrality = 0.0;
403 Double_t eventPlane = 0.0;
404 Double_t vertexz = fInputEvent->GetPrimaryVertex()->GetZ();
405 if(isAOD) {
406 AliAODHeader * header = static_cast<AliAODHeader*>(fInputEvent->GetHeader());
407 centrality = header->GetCentrality();
408 eventPlane = header->GetEventplane();
409 } else {
410 centrality = static_cast<AliESDEvent*>(fInputEvent)->GetCentrality()->GetCentralityPercentile("kV0M");
411 eventPlane = fInputEvent->GetEventplane()->GetEventplane("Q");
412 }
413
414 if(DebugLevel () > 4) {
415 cout << "centrality: " << centrality << " " << GetBin(fAxisCent, centrality) << endl;
416 cout << "vertexz: " << vertexz << " " << GetBin(fAxisZ, vertexz) << endl;
417 cout << "eventPlane: " << eventPlane << " " << endl;
418 }
419
420 const Int_t centBin = GetBin(fAxisCent, centrality);
421 const Int_t vertexBin = GetBin(fAxisZ, vertexz);
422
423 if(centBin < 0 || vertexBin < 0) {
424 AliError("bin out of range");
425 return;
426 }
427
428 fGammas->Clear();
429 fPions->Clear();
430
431 TClonesArray * aodGammas = GetConversionGammas(isAOD);
432 if(!aodGammas) {
433 AliError("no aod gammas found!");
434 return;
435 }
436
437 if(aodGammas->GetEntriesFast() > 0) {
438 if( static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(0) == fL1 &&
439 static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(1) == fL2
440 ) {
441 return;
442 }
443 fL1 = static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(0);
444 fL2 = static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(1);
445 //cout << aodGammas->GetEntriesFast() << " " << fInputEvent->GetNumberOfTracks() << "c" << endl;
446 }
447
448 if(DebugLevel() > 1) printf("Number of conversion gammas %d \n", aodGammas->GetEntriesFast());
449 for(Int_t ig = 0; ig < aodGammas->GetEntriesFast(); ig++) {
450 AliAODConversionPhoton * photon = dynamic_cast<AliAODConversionPhoton*>(aodGammas->At(ig));
451
452 if(!photon) continue;
453 if(!fV0Filter || fV0Filter->PhotonIsSelected(static_cast<AliConversionPhotonBase*>(photon), fInputEvent)) {
454 fGammas->Add(static_cast<TObject*>(photon));
455 }
456 }
457
458 if(DebugLevel() > 4) printf("Number of accepted gammas %d \n", fGammas->GetEntriesFast());
459
460 // THnSparseF * trackMehist = GetMEHistogram(vertexBin, centBin, hMETracks);
461 // hMEvents->Fill(vertexz, centrality);
462
463 ///Add tracks to array
464 TObjArray tracks;
465 for(Int_t iTrack = 0; iTrack < fInputEvent->GetNumberOfTracks(); iTrack++) {
466
467 AliVTrack * track = static_cast<AliVTrack*>(fInputEvent->GetTrack(iTrack));
468 //if(track->Pt() < 0.5) continue;
469 //if(TMath::Abs(track->Eta()) > 0.8) continue;
470 //if(track->GetTPCNcls() > 70)continue;
471 if(fTrackCuts.AcceptTrack(static_cast<AliAODTrack*>(track), static_cast<AliAODEvent*>(fInputEvent))) {
472 tracks.Add(track);
473 }
474 }
475
476 Process(fGammas, &tracks, vertexBin, centBin);
477
478 PostData(1, fHistograms);
479 PostData(2, fHistoGamma);
480 PostData(3, fHistoPion);
481
482}
483
484
485//________________________________________________________________________
486void AliAnalysisTaskdPhi::Process(TObjArray * gammas, TObjArray * tracks, Int_t vertexBin, Int_t centBin) {
487 ///Process stuff
488
489 if(DebugLevel() > 4) printf("Number of accepted tracks %d \n", tracks->GetEntriesFast());
490
491 AliAnaConvCorrBase * gCorr = GetCorrObject(vertexBin, centBin, fPhotonCorr);
492 AliAnaConvCorrBase * piCorr = GetCorrObject(vertexBin, centBin, fPionCorr);
493
494 if(!gCorr || !piCorr) {
495 AliError("corr object missing");
496 return;
497 }
498
499 for(Int_t i1 = 0; i1 < gammas->GetEntriesFast(); i1++) {
500 AliAODConversionPhoton * ph1 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i1));
501 Int_t tIDs[4] = {ph1->GetLabel(0), ph1->GetLabel(1), -1, -1};
502
503 Int_t leading = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs);
504 if(ph1->Pt() > fAxistPt.GetBinLowEdge(1)) {
505 gCorr->CorrelateWithTracks( static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs, leading);
506 }
507 for(Int_t i2 = 0; i2 < i1; i2++) {
508 AliAODConversionPhoton * ph2 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i2));
509
510 if( ph2->GetTrackLabelPositive()==ph1->GetTrackLabelPositive()
511 || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelNegative()
512 || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelPositive()
513 || ph2->GetTrackLabelPositive()==ph1->GetTrackLabelNegative()) {
514 continue;
515 }
516
517 AliAODConversionMother * pion = new AliAODConversionMother(ph1, ph2);
518 pion->SetLabels(i1, i2);
519
520 if(!fV0Filter || fV0Filter->MesonIsSelected(pion, kTRUE) ) {
521 Int_t leadingpi = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(pion), tracks, tIDs);
522 tIDs[2] = ph2->GetLabel(0);
523 tIDs[3] = ph2->GetLabel(1);
524 if(pion->Pt() > fAxistPt.GetBinLowEdge(1)) {
525 piCorr->CorrelateWithTracks(pion, tracks, tIDs, leadingpi);
526 }
527 }
528 }
529 }
530}
531
532//________________________________________________________________________
533void AliAnalysisTaskdPhi::Terminate(Option_t *) {
534
535 // Draw result to the screen
536 // Called once at the end of the query
537}
538
539//________________________________________________________________________
540TClonesArray * AliAnalysisTaskdPhi::GetConversionGammas(Bool_t isAOD) {
541
542
543 if(isAOD) {
544
545 TClonesArray * gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
546 if(gammas) {
547 return gammas;
548 }
549
550 FindDeltaAODBranchName(fInputEvent);
551 gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
552 return gammas;
553
554 } else {
555 TClonesArray * gammas = dynamic_cast<TClonesArray*>(GetInputData(1));
556 return gammas;
557 }
558
559}
560
561//________________________________________________________________________
562void AliAnalysisTaskdPhi::FindDeltaAODBranchName(AliVEvent * event){
563 ///Find aod branch
564 TList *list=event->GetList();
565 for(Int_t ii=0;ii<list->GetEntries();ii++){
566 TString name((list->At(ii))->GetName());
567 if(name.BeginsWith("GammaConv")&&name.EndsWith("gamma")){
568 fDeltaAODBranchName=name;
569 AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
570 return;
571 }
572 }
573}
574
575