adding track cut filter. Only add used ME tracks to ME
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskdPhi.cxx
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>
44 using namespace std;
45
46 ClassImp(AliAnalysisTaskdPhi)
47
48
49 //________________________________________________________________________
50 AliAnalysisTaskdPhi::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 //________________________________________________________________________
117 AliAnalysisTaskdPhi::~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
139 if(fHistoPion)
140         delete fHistoPion;
141   fHistoPion = NULL;
142
143   if(fHistoGamma)
144         delete fHistoGamma;
145   fHistoGamma = NULL;
146
147 }
148
149 ///________________________________________________________________________
150 void 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 //________________________________________________________________________
216 void 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 ///________________________________________________________________________
332 THnSparseF * 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 //________________________________________________________________________
372 void 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 //________________________________________________________________________
486 void 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 //________________________________________________________________________
533 void AliAnalysisTaskdPhi::Terminate(Option_t *) {
534  
535   // Draw result to the screen
536   // Called once at the end of the query
537 }
538
539 //________________________________________________________________________
540 TClonesArray * 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 //________________________________________________________________________
562 void 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