changes from gsi. Using mult if no centrality. testfilterbit 128
[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(NULL),
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   fDeltaAODBranchName("AliAODGammaConversion_gamma"), 
66   fAxistPt(),
67   fAxiscPt(),
68   fAxisEta(),
69   fAxisPhi(),
70   fAxisCent(),
71   fAxisZ(), 
72   fAxisPiM()
73 {
74   //constructor
75   fAxistPt.SetNameTitle("tPtAxis", "tPt");
76   fAxistPt.Set(20, 0, 100);
77
78   fAxiscPt.SetNameTitle("cPtAxis", "cPt");
79   fAxiscPt.Set(20, 0, 100);
80
81   fAxisEta.SetNameTitle("EtaAxis", "Eta");
82   fAxisEta.Set(160, -0.8, 0.8);
83
84   fAxisPhi.SetNameTitle("PhiAxis", "Phi");
85   fAxisPhi.Set(128, 0, TMath::TwoPi());
86
87   fAxisZ.SetNameTitle("ZAxis", "Z");
88   fAxisZ.Set(4, -10, 10);
89
90   fAxisCent.SetNameTitle("CentAxis", "Cent");
91   Double_t centbins[5] = {0, 10, 30, 60, 100.1};
92   fAxisCent.Set(4, centbins);
93
94   Double_t mbins[8] = {0.11, 0.12, 0.13, 0.132, 0.138, 0.14, 0.15, 0.16};
95   fAxisPiM.SetNameTitle("InvMassPi0", "Invariant mass");
96   fAxisPiM.Set(7, mbins);
97
98   fGammas = new TObjArray();
99   fGammas->SetOwner(kFALSE);
100
101   fPions = new TObjArray();
102   fPions->SetOwner(kFALSE);
103
104   // Define input and output slots here
105   DefineInput(0, TChain::Class());
106   //DefineInput(1, TClonesArray::Class());
107   DefineOutput(1, TList::Class());
108   DefineOutput(2, TList::Class());
109   DefineOutput(3, TList::Class());
110 }
111
112
113
114 //________________________________________________________________________
115 AliAnalysisTaskdPhi::~AliAnalysisTaskdPhi(){
116   //destructor
117   if(fPions)
118         delete fPions;
119   fPions = NULL;
120
121   if(fGammas)
122         delete fGammas;
123   fGammas = NULL;
124   
125   if(fIsoAna)
126         delete fIsoAna;
127   fIsoAna = NULL;
128
129   if(fV0Filter)
130         delete fV0Filter;
131   fV0Filter = NULL;
132
133   if(fHistograms)
134         delete fHistograms;
135   fHistograms = NULL;
136
137   if(fHistoPion)
138         delete fHistoPion;
139   fHistoPion = NULL;
140
141   if(fHistoGamma)
142         delete fHistoGamma;
143   fHistoGamma = NULL;
144
145 }
146
147 ///________________________________________________________________________
148 void AliAnalysisTaskdPhi::SetUpCorrObjects() {
149   ///Creat corr obj
150   fIsoAna = new AliAnaConvIsolation();
151
152
153   fPhotonCorr = new TObjArray();
154   fPionCorr = new TObjArray();
155   
156   TList * hPhoton = new TList();
157   hPhoton->SetName("hPhotonCorr");
158   hPhoton->SetOwner(kTRUE);
159   fHistoGamma->Add(hPhoton);
160
161   TList * hPion = new TList();
162   hPion->SetName("hPionCorr");
163   hPion->SetOwner(kTRUE);
164   fHistoPion->Add(hPion);
165
166
167         for(Int_t ic = 0; ic < fAxisCent.GetNbins(); ic++) {
168         TObjArray * photonArray = new TObjArray();
169         photonArray->SetOwner(kTRUE);
170         fPhotonCorr->AddAt(photonArray, ic);
171
172         TObjArray * pionArray = new TObjArray();
173         pionArray->SetOwner(kTRUE);
174         fPionCorr->AddAt(pionArray, ic);
175
176         TList * photonList = new TList();
177         photonList->SetName(Form("photon_%d", ic));
178         photonList->SetOwner(kTRUE);
179         hPhoton->AddAt(photonList, ic);
180
181         TList * pionList = new TList();
182         pionList->SetName(Form("pion_%d", ic));
183         pionList->SetOwner(kTRUE);
184         hPion->AddAt(pionList, ic);
185         
186           
187         for(Int_t iz = 0; iz < fAxisZ.GetNbins(); iz++) {
188           TString nameString = Form("%d_%d", ic, iz);
189           TString titleString = Form("%f < Z < %f ... %f cent %f", 
190                                                                  fAxisZ.GetBinLowEdge(iz+1), fAxisZ.GetBinUpEdge(iz+1), 
191                                                                  fAxisCent.GetBinLowEdge(ic+1), fAxisCent.GetBinUpEdge(ic+1));
192
193
194
195           AliAnaConvCorrPhoton * photonCorr = new AliAnaConvCorrPhoton(Form("PhotonCorr_%s", nameString.Data()), Form("photon %s", titleString.Data()));
196           photonArray->AddAt(photonCorr, iz);
197           photonCorr->GetAxistPt().Set(fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
198           photonCorr->GetAxiscPt().Set(fAxiscPt.GetNbins(), fAxiscPt.GetXbins()->GetArray());
199           photonCorr->CreateHistograms();
200           photonList->Add(photonCorr->GetHistograms());
201
202           AliAnaConvCorrPion * pionCorr = new AliAnaConvCorrPion(Form("PionCorr_%s", nameString.Data()), Form("pion %s", titleString.Data()));
203           pionArray->AddAt(pionCorr, iz);
204           pionCorr->GetAxistPt().Set(fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
205           pionCorr->GetAxiscPt().Set(fAxiscPt.GetNbins(), fAxiscPt.GetXbins()->GetArray());
206           pionCorr->GetAxisM().Set(fAxisPiM.GetNbins(), fAxisPiM.GetXbins()->GetArray());
207           pionCorr->CreateHistograms();
208           pionList->Add(pionCorr->GetHistograms());
209         }
210   }
211 }
212
213 //________________________________________________________________________
214 void AliAnalysisTaskdPhi::UserCreateOutputObjects() {
215   // Create histograms
216   
217   fHistograms = new TList();
218   fHistograms->SetName("dPhi_histograms");
219   fHistograms->SetOwner(kTRUE);
220
221   fHistoGamma = new TList();
222   fHistoGamma->SetName("Gamma_histo");
223   fHistoGamma->SetOwner(kTRUE);
224
225   fHistoPion = new TList();
226   fHistoPion->SetName("Pion_histo");
227   fHistoPion->SetOwner(kTRUE);
228
229   
230   if(fV0Filter) {
231         fV0Filter->InitCutHistograms();
232         fHistograms->Add(fV0Filter->GetCutHistograms());
233   }
234   
235
236   AliConversionTrackCuts * tc = dynamic_cast<AliConversionTrackCuts*>(fTrackCuts);
237   if(tc) fHistograms->Add(tc->CreateHistograms());
238
239   SetUpCorrObjects();
240
241
242   ///Set up ME histograms
243   TList * MEHistograms = new TList();
244   MEHistograms->SetName("MEHistograms");
245   MEHistograms->SetOwner(kTRUE);
246   fHistograms->Add(MEHistograms);
247
248
249   hMEvents = new TH2I("hMEvents", "Nevents vs centrality vertexz",
250                                           fAxisZ.GetNbins(), fAxisZ.GetBinLowEdge(1), fAxisZ.GetBinUpEdge(fAxisZ.GetNbins()),
251                                           fAxisCent.GetNbins(), fAxisCent.GetBinLowEdge(1), fAxisCent.GetBinUpEdge(fAxisCent.GetNbins()));
252   hMEvents->GetYaxis()->Set(fAxisCent.GetNbins(), fAxisCent.GetXbins()->GetArray());
253   MEHistograms->Add(hMEvents);
254
255   TList * outAxesList = new TList();
256   outAxesList->Add(&fAxisCent);
257   outAxesList->Add(&fAxisZ);
258   fHistograms->Add(outAxesList);
259
260   PostData(1, fHistograms);
261   PostData(2, fHistoGamma);
262   PostData(3, fHistoPion);
263
264 }
265
266 ///________________________________________________________________________
267 THnSparseF * AliAnalysisTaskdPhi::CreateSparse(TString nameString, TString titleString, TList * axesList) {
268   ///Create sparse
269   const Int_t dim = axesList->GetSize();
270
271   TAxis * axes[dim];
272   Int_t bins[dim];
273   Double_t min[dim];
274   Double_t max[dim];
275
276   for(Int_t i = 0; i<dim; i++) {
277         TAxis * axis = dynamic_cast<TAxis*>(axesList->At(i));
278         if(axis) {
279           axes[i] = axis;
280         } else {
281           cout << "AliAnalysisTaskdPhi::CreateSparse: Error error, all the axes are not present in axis list" << endl;
282           return NULL;
283         }
284   }
285
286   for(Int_t i = 0; i<dim; i++) {
287         bins[i] = axes[i]->GetNbins(); 
288         min[i] = axes[i]->GetBinLowEdge(1);
289         max[i] = axes[i]->GetBinUpEdge(axes[i]->GetNbins());
290   }
291
292   THnSparseF * sparse = new THnSparseF(Form("METracks_%s", nameString.Data()), 
293                                                                                            Form("tracks %s", titleString.Data()), 
294                                                                                            dim, bins, min, max);
295   
296   for(Int_t i = 0; i<dim; i++) {
297         sparse->GetAxis(i)->SetNameTitle(axes[i]->GetName(), axes[i]->GetTitle() );
298         if(axes[i]->GetXbins()->GetSize() > 0) {
299           sparse->SetBinEdges(i, axes[i]->GetXbins()->GetArray() );
300         }
301   }
302
303   return sparse;
304 }
305
306 //________________________________________________________________________
307 void AliAnalysisTaskdPhi::UserExec(Option_t *) {
308   ///User exec. 
309
310   //if(! fV0Filter->EventIsSelected(fInputEvent)) return;
311
312   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
313   Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
314   
315   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
316   if (!inputHandler) {
317         cout << "cout no input event handler"<<endl;
318         return;
319   }
320
321   
322   if ( fV0Filter && !fV0Filter->GetPIDResponse() ) {
323         if ( inputHandler->GetPIDResponse() ){
324           fV0Filter->SetPIDResponse( inputHandler->GetPIDResponse() );
325         } else {
326           
327           //AOD case
328           if (isAOD){
329                 if (!fV0Filter->GetPIDResponse()){
330                   fV0Filter->InitAODpidUtil(1);
331                 }
332           }
333         }
334   }
335
336   Double_t centrality = 0.0;
337   Double_t multiplicity = 0.0;
338   Double_t eventPlane = 0.0;
339   Double_t vertexz = fInputEvent->GetPrimaryVertex()->GetZ();
340   
341   if(isAOD) {
342     AliAODHeader * header = static_cast<AliAODHeader*>(fInputEvent->GetHeader());
343         centrality = header->GetCentrality();
344         eventPlane = header->GetEventplane();
345         multiplicity = header->GetRefMultiplicity();
346         if(centrality < 0) centrality = multiplicity;
347   } else {
348         centrality = static_cast<AliESDEvent*>(fInputEvent)->GetCentrality()->GetCentralityPercentile("V0M");
349         eventPlane = fInputEvent->GetEventplane()->GetEventplane("Q");
350   }
351
352
353   const Int_t centBin = GetBin(fAxisCent, centrality);
354   const Int_t vertexBin = GetBin(fAxisZ, vertexz);
355
356
357   if(DebugLevel () > 4) {
358         cout << "centrality: " << centrality <<  " " << GetBin(fAxisCent, centrality) << endl;
359         cout << "vertexz: " << vertexz <<  " " << GetBin(fAxisZ, vertexz) << endl;
360         cout << "eventPlane: " << eventPlane <<  " " << endl;
361         cout << "multiplicity: "<<  multiplicity << endl;
362
363   }
364
365
366   if(centBin < 0 || vertexBin < 0) {
367         AliError("bin out of range");
368         //cout << "bad bin"<<endl;
369         return;
370   }
371
372   fGammas->Clear();
373   fPions->Clear();
374
375   TClonesArray * aodGammas = GetConversionGammas(isAOD);
376
377   if(!aodGammas) {
378         AliError("no aod gammas found!");
379         return;
380   }
381
382
383   if(DebugLevel() > 1) printf("Number of conversion gammas %d \n", aodGammas->GetEntriesFast());
384   for(Int_t ig = 0; ig < aodGammas->GetEntriesFast(); ig++) {
385    AliAODConversionPhoton * photon = dynamic_cast<AliAODConversionPhoton*>(aodGammas->At(ig));
386     
387     if(!photon) {
388           cout << "can't get photon"<<endl;
389           continue;
390         }
391         if(!VerifyAODGamma(photon)) { 
392           continue;
393         }
394
395     if(!fV0Filter || fV0Filter->PhotonIsSelected(static_cast<AliConversionPhotonBase*>(photon), fInputEvent)) {
396       fGammas->Add(static_cast<TObject*>(photon));
397     }
398   }
399   
400   if(DebugLevel() > 4) printf("Number of accepted gammas %d \n", fGammas->GetEntriesFast());
401   hMEvents->Fill(vertexz, centrality);
402   
403   
404   
405   ///create track array
406   TObjArray tracks;
407   const Double_t etalim[2] = { fAxisEta.GetBinLowEdge(1), fAxisEta.GetBinUpEdge(fAxisEta.GetNbins())};
408   for(Int_t iTrack = 0; iTrack < fInputEvent->GetNumberOfTracks(); iTrack++) {
409
410     AliVTrack * track = static_cast<AliVTrack*>(fInputEvent->GetTrack(iTrack));
411     if(track->Pt() < fAxiscPt.GetBinLowEdge(1) ) continue;
412     if(track->Eta() < etalim[0] || track->Eta() > etalim[1]) continue;
413     
414     
415     if(!fTrackCuts || fTrackCuts->IsSelected((track))) {
416       tracks.Add(track);
417       //cout <<"a"<<endl;
418     } // else {
419     //       cout <<"b"<<endl;
420     //     }
421   }
422
423   Process(fGammas, &tracks, vertexBin, centBin);
424
425   PostData(1, fHistograms);
426   PostData(2, fHistoGamma);
427   PostData(3, fHistoPion);
428   
429 }
430
431
432 //________________________________________________________________________
433 void AliAnalysisTaskdPhi::Process(TObjArray * gammas, TObjArray * tracks, Int_t vertexBin, Int_t centBin) {
434   ///Process stuff
435
436   if(DebugLevel() > 4) printf("Number of accepted gammas, tracks %d  %d \n", gammas->GetEntriesFast(), tracks->GetEntriesFast());
437
438  
439   AliAnaConvCorrBase * gCorr = GetCorrObject(vertexBin, centBin, fPhotonCorr);
440   AliAnaConvCorrPion * piCorr = dynamic_cast<AliAnaConvCorrPion*>(GetCorrObject(vertexBin, centBin, fPionCorr));
441   
442   if(!gCorr || !piCorr) {
443         AliError("corr object missing");
444         return;
445   }
446
447   for(Int_t i1 = 0; i1 < gammas->GetEntriesFast(); i1++) {
448         AliAODConversionPhoton * ph1 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i1));
449         Int_t tIDs[4] = {ph1->GetLabel(0), ph1->GetLabel(1), -1, -1};
450
451         Int_t leading = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs);
452         if(ph1->Pt() > fAxistPt.GetBinLowEdge(1)) {
453           gCorr->CorrelateWithTracks( static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs, leading);
454         }
455         for(Int_t i2 = 0; i2 < i1; i2++) {
456           AliAODConversionPhoton * ph2 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i2));
457           
458           if( ph2->GetTrackLabelPositive()==ph1->GetTrackLabelPositive() 
459                   || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelNegative()
460                   || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelPositive()
461                   || ph2->GetTrackLabelPositive()==ph1->GetTrackLabelNegative()) {
462                 continue;
463           }
464
465           AliAODConversionMother * pion = new AliAODConversionMother(ph1, ph2);
466           pion->SetLabels(i1, i2);
467           
468           if(!fV0Filter || fV0Filter->MesonIsSelected(pion, kTRUE) ) {
469         
470                 Int_t leadingpi = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(pion), tracks, tIDs);
471                 piCorr->FillTriggerCounters(pion, leadingpi);
472                 
473                 tIDs[2] = ph2->GetLabel(0);
474                 tIDs[3] = ph2->GetLabel(1);
475                 if(pion->Pt() > fAxistPt.GetBinLowEdge(1) && 
476                    pion->M() > fAxisPiM.GetBinLowEdge(1) && 
477                    pion->M() < fAxisPiM.GetBinUpEdge(fAxisPiM.GetNbins())) {
478                   piCorr->CorrelateWithTracks(pion, tracks, tIDs, leadingpi);
479                 }
480           }
481         }
482   }
483 }
484
485
486 //________________________________________________________________________
487 void AliAnalysisTaskdPhi::Terminate(Option_t *) {
488  
489   // Draw result to the screen
490   // Called once at the end of the query
491 }
492
493 //________________________________________________________________________
494 TClonesArray * AliAnalysisTaskdPhi::GetConversionGammas(Bool_t isAOD) {
495   
496
497   if(isAOD) {
498
499         TClonesArray * gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
500         if(gammas) {
501           return gammas;
502         }
503
504         FindDeltaAODBranchName(fInputEvent);
505         gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
506         return gammas;
507
508   } else {
509         TClonesArray * gammas = dynamic_cast<TClonesArray*>(GetInputData(1));
510         return gammas;
511   }
512
513 }
514
515 //________________________________________________________________________
516 void AliAnalysisTaskdPhi::FindDeltaAODBranchName(AliVEvent * event){
517   ///Find aod branch
518   TList *list=event->GetList();
519   for(Int_t ii=0;ii<list->GetEntries();ii++){
520         TString name((list->At(ii))->GetName());
521         if(name.BeginsWith("GammaConv")&&name.EndsWith("gamma")){
522           fDeltaAODBranchName=name;
523           AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
524           cout <<fDeltaAODBranchName << endl;
525           return;
526         }
527   }
528 }
529   
530
531
532
533 ///________________________________________________________________________
534 Bool_t AliAnalysisTaskdPhi::VerifyAODGamma(AliAODConversionPhoton * gamma) {
535
536   AliAODEvent * event = static_cast<AliAODEvent*>(fInputEvent);
537
538   //cout << "label "<< gamma->GetV0Index() << endl;
539
540   AliAODv0 * v0 =  NULL;
541
542   //Int_t v0idx =  gamma->GetV0Index();
543   for(Int_t i = 0; i < event->GetNumberOfV0s(); i++) {
544         AliAODv0 * tv0 = event->GetV0(i);
545
546         //cout << i << " " << v0->GetID() << " " << v0->GetSecondaryVtx()->GetID() << " " << v0->GetLabel() << " " << v0->GetSecondaryVtx()->GetLabel() << endl;
547    
548         if(tv0->GetSecondaryVtx()->GetID() == gamma->GetV0Index() ) {
549           v0 = tv0;
550           //cout << "found it" << endl;
551           break;
552         }       
553   }
554
555   if(!v0) {
556         //cout << "v0 not found"<<endl;
557         return kFALSE;
558   } 
559
560
561   AliAODTrack * d1 = dynamic_cast<AliAODTrack*>(v0->GetDaughter(0));
562   AliAODTrack * d2 = dynamic_cast<AliAODTrack*>(v0->GetDaughter(1));
563
564   Int_t t1 = -1;
565   Int_t t2 = -2;
566
567   if(d1) t1 = d1->GetID();
568   if(d2) t2 = d2->GetID();
569
570   Int_t g1 = gamma->GetLabel(0);
571   Int_t g2 = gamma->GetLabel(1);
572
573   if((t1 == g1 && 
574           t2 == g2) || 
575          (t1 == g2 && 
576           t2 == g1) ) {
577         //cout <<"match"<< " " <<  gamma->Pt() << " " <<  d1->Pt() + d2->Pt() <<endl;
578   }
579                  
580   else {
581
582         cout << g1 << " " << g2 <<endl;
583         cout << t1 << " " << t2 <<endl;
584         
585         for(Int_t i = 0; i < event->GetNumberOfV0s(); i++) {
586         v0 = event->GetV0(i);
587         cout << i << " " << v0->GetSecondaryVtx()->GetID() << " " <<dynamic_cast<AliAODTrack*>(v0->GetDaughter(0))->GetID() << " " << dynamic_cast<AliAODTrack*>(v0->GetDaughter(1))->GetID() << endl; 
588         
589         }
590   }
591   
592   // Float_t sumdpt = d1->Pt() + d2->Pt();
593   // cout << "pt: " << sumdpt << " " << gamma->Pt() << endl;
594
595   return kTRUE;
596
597
598
599 }