[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskdPhi.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *																          *
 * Authors: Svein Lindal                                                  *
 * Version 1.0			                                                  *
 *																          *
 *																          *
 * Permission to use, copy, modify and distribute this software and its	  *
 * documentation strictly for non-commercial purposes is hereby granted	  *
 * without fee, provided that the above copyright notice appears in all	  *
 * copies and that both the copyright notice and this permission notice	  *
 * appear in the supporting documentation. The authors make no claims	  *
 * about the suitability of this software for any purpose. It is		  *
 * provided "as is" without express or implied warranty.				  *
 **************************************************************************/

////////////////////////////////////////////////
//--------------------------------------------- 
// Class doing conversion gamma dPhi correlations
// Gamma Conversion analysis
//---------------------------------------------
////////////////////////////////////////////////

#include "AliAnalysisTaskdPhi.h"

#include <TH2I.h>
#include <TList.h>
#include <TChain.h>

#include <AliAnalysisManager.h>
#include <AliInputEventHandler.h>
#include <AliESDInputHandler.h>
#include <AliAODInputHandler.h>
#include <AliAnalysisFilter.h>

#include "AliConversionTrackCuts.h"
#include "AliConversionCuts.h"
#include "AliAODConversionPhoton.h"
#include "AliAODConversionMother.h"
#include "AliAnaConvCorrPhoton.h"
#include "AliAnaConvCorrPion.h"
#include "AliAnaConvIsolation.h"
// Author Svein Lindal <slindal@fys.uio.no>
using namespace std;

ClassImp(AliAnalysisTaskdPhi)


//________________________________________________________________________
AliAnalysisTaskdPhi::AliAnalysisTaskdPhi(const char *name) : AliAnalysisTaskSE(name),
  fHistograms(NULL),
  fHistoGamma(NULL),
  fHistoPion(NULL),
  fV0Filter(NULL),
  fTrackCuts(NULL),
  fGammas(NULL),
  fPions(NULL),
  hMETracks(NULL), 
  hMEPhotons(NULL), 
  hMEPions(NULL),
  hMEvents(NULL),
  fPhotonCorr(NULL),
  fPionCorr(NULL), 
  fIsoAna(NULL),
  fDeltaAODBranchName("AliAODGammaConversion_gamma"), 
  fAxistPt(),
  fAxiscPt(),
  fAxisEta(),
  fAxisPhi(),
  fAxisCent(),
  fAxisZ(), 
  fAxisPiM()
{
  //constructor
  fAxistPt.SetNameTitle("tPtAxis", "tPt");
  fAxistPt.Set(20, 0, 100);

  fAxiscPt.SetNameTitle("cPtAxis", "cPt");
  fAxiscPt.Set(20, 0, 100);

  fAxisEta.SetNameTitle("EtaAxis", "Eta");
  fAxisEta.Set(160, -0.8, 0.8);

  fAxisPhi.SetNameTitle("PhiAxis", "Phi");
  fAxisPhi.Set(128, 0, TMath::TwoPi());

  fAxisZ.SetNameTitle("ZAxis", "Z");
  fAxisZ.Set(4, -10, 10);

  fAxisCent.SetNameTitle("CentAxis", "Cent");
  Double_t centbins[5] = {0, 10, 30, 60, 100.1};
  fAxisCent.Set(4, centbins);

  Double_t mbins[8] = {0.11, 0.12, 0.13, 0.132, 0.138, 0.14, 0.15, 0.16};
  fAxisPiM.SetNameTitle("InvMassPi0", "Invariant mass");
  fAxisPiM.Set(7, mbins);

  fGammas = new TObjArray();
  fGammas->SetOwner(kFALSE);

  fPions = new TObjArray();
  fPions->SetOwner(kFALSE);

  // Define input and output slots here
  DefineInput(0, TChain::Class());
  //DefineInput(1, TClonesArray::Class());
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
  DefineOutput(3, TList::Class());
}


//________________________________________________________________________
AliAnalysisTaskdPhi::~AliAnalysisTaskdPhi(){
  //destructor
  if(fPions)
	delete fPions;
  fPions = NULL;

  if(fGammas)
	delete fGammas;
  fGammas = NULL;
  
  if(fIsoAna)
  	delete fIsoAna;
  fIsoAna = NULL;

  if(fV0Filter)
	delete fV0Filter;
  fV0Filter = NULL;

  if(fHistograms)
	delete fHistograms;
  fHistograms = NULL;

  if(fHistoPion)
	delete fHistoPion;
  fHistoPion = NULL;

  if(fHistoGamma)
	delete fHistoGamma;
  fHistoGamma = NULL;

}

///________________________________________________________________________
void AliAnalysisTaskdPhi::SetUpCorrObjects() {
  ///Creat corr obj
  fIsoAna = new AliAnaConvIsolation();


  fPhotonCorr = new TObjArray();
  fPionCorr = new TObjArray();
  
  TList * hPhoton = new TList();
  hPhoton->SetName("hPhotonCorr");
  hPhoton->SetOwner(kTRUE);
  fHistoGamma->Add(hPhoton);

  TList * hPion = new TList();
  hPion->SetName("hPionCorr");
  hPion->SetOwner(kTRUE);
  fHistoPion->Add(hPion);


	for(Int_t ic = 0; ic < fAxisCent.GetNbins(); ic++) {
	TObjArray * photonArray = new TObjArray();
	photonArray->SetOwner(kTRUE);
	fPhotonCorr->AddAt(photonArray, ic);

	TObjArray * pionArray = new TObjArray();
	pionArray->SetOwner(kTRUE);
	fPionCorr->AddAt(pionArray, ic);

	TList * photonList = new TList();
	photonList->SetName(Form("photon_%d", ic));
	photonList->SetOwner(kTRUE);
	hPhoton->AddAt(photonList, ic);

	TList * pionList = new TList();
	pionList->SetName(Form("pion_%d", ic));
	pionList->SetOwner(kTRUE);
	hPion->AddAt(pionList, ic);
	
	  
	for(Int_t iz = 0; iz < fAxisZ.GetNbins(); iz++) {
	  TString nameString = Form("%d_%d", ic, iz);
	  TString titleString = Form("%f < Z < %f ... %f cent %f", 
								 fAxisZ.GetBinLowEdge(iz+1), fAxisZ.GetBinUpEdge(iz+1), 
								 fAxisCent.GetBinLowEdge(ic+1), fAxisCent.GetBinUpEdge(ic+1));


	  AliAnaConvCorrPhoton * photonCorr = new AliAnaConvCorrPhoton(Form("PhotonCorr_%s", nameString.Data()), Form("photon %s", titleString.Data()));
	  photonArray->AddAt(photonCorr, iz);
	  photonCorr->GetAxistPt().Set(fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
	  photonCorr->GetAxiscPt().Set(fAxiscPt.GetNbins(), fAxiscPt.GetXbins()->GetArray());
	  photonCorr->CreateHistograms();
	  photonList->Add(photonCorr->GetHistograms());

	  AliAnaConvCorrPion * pionCorr = new AliAnaConvCorrPion(Form("PionCorr_%s", nameString.Data()), Form("pion %s", titleString.Data()));
	  pionArray->AddAt(pionCorr, iz);
	  pionCorr->GetAxistPt().Set(fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
	  pionCorr->GetAxiscPt().Set(fAxiscPt.GetNbins(), fAxiscPt.GetXbins()->GetArray());
	  pionCorr->GetAxisM().Set(fAxisPiM.GetNbins(), fAxisPiM.GetXbins()->GetArray());
	  pionCorr->CreateHistograms();
	  pionList->Add(pionCorr->GetHistograms());
	}
  }
}

//________________________________________________________________________
void AliAnalysisTaskdPhi::UserCreateOutputObjects() {
  // Create histograms
  
  fHistograms = new TList();
  fHistograms->SetName("dPhi_histograms");
  fHistograms->SetOwner(kTRUE);

  fHistoGamma = new TList();
  fHistoGamma->SetName("Gamma_histo");
  fHistoGamma->SetOwner(kTRUE);

  fHistoPion = new TList();
  fHistoPion->SetName("Pion_histo");
  fHistoPion->SetOwner(kTRUE);

  
  if(fV0Filter) {
	fV0Filter->InitCutHistograms();
	fHistograms->Add(fV0Filter->GetCutHistograms());
  }
  

  AliConversionTrackCuts * tc = dynamic_cast<AliConversionTrackCuts*>(fTrackCuts);
  if(tc) fHistograms->Add(tc->CreateHistograms());

  SetUpCorrObjects();


  ///Set up ME histograms
  TList * MEHistograms = new TList();
  MEHistograms->SetName("MEHistograms");
  MEHistograms->SetOwner(kTRUE);
  fHistograms->Add(MEHistograms);


  hMEvents = new TH2I("hMEvents", "Nevents vs centrality vertexz",
					  fAxisZ.GetNbins(), fAxisZ.GetBinLowEdge(1), fAxisZ.GetBinUpEdge(fAxisZ.GetNbins()),
 					  fAxisCent.GetNbins(), fAxisCent.GetBinLowEdge(1), fAxisCent.GetBinUpEdge(fAxisCent.GetNbins()));
  hMEvents->GetYaxis()->Set(fAxisCent.GetNbins(), fAxisCent.GetXbins()->GetArray());
  MEHistograms->Add(hMEvents);

  TList * outAxesList = new TList();
  outAxesList->Add(&fAxisCent);
  outAxesList->Add(&fAxisZ);
  fHistograms->Add(outAxesList);

  PostData(1, fHistograms);
  PostData(2, fHistoGamma);
  PostData(3, fHistoPion);

}

///________________________________________________________________________
THnSparseF * AliAnalysisTaskdPhi::CreateSparse(TString nameString, TString titleString, TList * axesList) {
  ///Create sparse
  const Int_t dim = axesList->GetSize();

  TAxis * axes[dim];
  Int_t bins[dim];
  Double_t min[dim];
  Double_t max[dim];

  for(Int_t i = 0; i<dim; i++) {
	TAxis * axis = dynamic_cast<TAxis*>(axesList->At(i));
	if(axis) {
	  axes[i] = axis;
  	} else {
	  cout << "AliAnalysisTaskdPhi::CreateSparse: Error error, all the axes are not present in axis list" << endl;
	  return NULL;
	}
  }

  for(Int_t i = 0; i<dim; i++) {
	bins[i] = axes[i]->GetNbins(); 
	min[i] = axes[i]->GetBinLowEdge(1);
	max[i] = axes[i]->GetBinUpEdge(axes[i]->GetNbins());
  }

  THnSparseF * sparse = new THnSparseF(Form("METracks_%s", nameString.Data()), 
											   Form("tracks %s", titleString.Data()), 
											   dim, bins, min, max);
  
  for(Int_t i = 0; i<dim; i++) {
	sparse->GetAxis(i)->SetNameTitle(axes[i]->GetName(), axes[i]->GetTitle() );
	if(axes[i]->GetXbins()->GetSize() > 0) {
	  sparse->SetBinEdges(i, axes[i]->GetXbins()->GetArray() );
	}
  }

  return sparse;
}

//________________________________________________________________________
void AliAnalysisTaskdPhi::UserExec(Option_t *) {
  ///User exec. 

  //if(! fV0Filter->EventIsSelected(fInputEvent)) return;

  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
  
  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
  if (!inputHandler) {
	cout << "cout no input event handler"<<endl;
	return;
  }

  
  if ( fV0Filter && !fV0Filter->GetPIDResponse() ) {
	if ( inputHandler->GetPIDResponse() ){
	  fV0Filter->SetPIDResponse( inputHandler->GetPIDResponse() );
	} else {
	  
	  //AOD case
	  if (isAOD){
		if (!fV0Filter->GetPIDResponse()){
		  fV0Filter->InitAODpidUtil(1);
		}
	  }
	}
  }

  Double_t centrality = 0.0;
  Double_t eventPlane = 0.0;
  Double_t vertexz = fInputEvent->GetPrimaryVertex()->GetZ();
  if(isAOD) {
    AliAODHeader * header = static_cast<AliAODHeader*>(fInputEvent->GetHeader());
	centrality = header->GetCentrality();
	eventPlane = header->GetEventplane();
  } else {
	centrality = static_cast<AliESDEvent*>(fInputEvent)->GetCentrality()->GetCentralityPercentile("V0M");
	eventPlane = fInputEvent->GetEventplane()->GetEventplane("Q");
  }


  const Int_t centBin = GetBin(fAxisCent, centrality);
  const Int_t vertexBin = GetBin(fAxisZ, vertexz);


  if(DebugLevel () > 4) {
	cout << "centrality: " << centrality <<  " " << GetBin(fAxisCent, centrality) << endl;
	cout << "vertexz: " << vertexz <<  " " << GetBin(fAxisZ, vertexz) << endl;
	cout << "eventPlane: " << eventPlane <<  " " << endl;
  }


  if(centBin < 0 || vertexBin < 0) {
	AliError("bin out of range");
	//cout << "bad bin"<<endl;
	return;
  }

  fGammas->Clear();
  fPions->Clear();

  TClonesArray * aodGammas = GetConversionGammas(isAOD);

  if(!aodGammas) {
	AliError("no aod gammas found!");
	return;
  }


  if(DebugLevel() > 1) printf("Number of conversion gammas %d \n", aodGammas->GetEntriesFast());
  for(Int_t ig = 0; ig < aodGammas->GetEntriesFast(); ig++) {
   AliAODConversionPhoton * photon = dynamic_cast<AliAODConversionPhoton*>(aodGammas->At(ig));
    
    if(!photon) {
	  cout << "can't get photon"<<endl;
	  continue;
	}
	if(!VerifyAODGamma(photon)) { 
	  continue;
	}

    if(!fV0Filter || fV0Filter->PhotonIsSelected(static_cast<AliConversionPhotonBase*>(photon), fInputEvent)) {
      fGammas->Add(static_cast<TObject*>(photon));
    }
  }
  
  if(DebugLevel() > 4) printf("Number of accepted gammas %d \n", fGammas->GetEntriesFast());
  hMEvents->Fill(vertexz, centrality);
  
  
  ///create track array
  TObjArray tracks;
  const Double_t etalim[2] = { fAxisEta.GetBinLowEdge(1), fAxisEta.GetBinUpEdge(fAxisEta.GetNbins())};
  for(Int_t iTrack = 0; iTrack < fInputEvent->GetNumberOfTracks(); iTrack++) {

	AliVTrack * track = static_cast<AliVTrack*>(fInputEvent->GetTrack(iTrack));
	if(track->Pt() < fAxiscPt.GetBinLowEdge(1) ) continue;
	if(track->Eta() < etalim[0] || track->Eta() > etalim[1]) continue;

	
	if(!fTrackCuts || fTrackCuts->IsSelected((track))) {
	  tracks.Add(track);
	}
  }
  
  Process(fGammas, &tracks, vertexBin, centBin);

  PostData(1, fHistograms);
  PostData(2, fHistoGamma);
  PostData(3, fHistoPion);
  
}


//________________________________________________________________________
void AliAnalysisTaskdPhi::Process(TObjArray * gammas, TObjArray * tracks, Int_t vertexBin, Int_t centBin) {
  ///Process stuff

  if(DebugLevel() > 4) printf("Number of accepted gammas, tracks %d  %d \n", gammas->GetEntriesFast(), tracks->GetEntriesFast());

 
  AliAnaConvCorrBase * gCorr = GetCorrObject(vertexBin, centBin, fPhotonCorr);
  AliAnaConvCorrPion * piCorr = dynamic_cast<AliAnaConvCorrPion*>(GetCorrObject(vertexBin, centBin, fPionCorr));
  
  if(!gCorr || !piCorr) {
	AliError("corr object missing");
	return;
  }

  for(Int_t i1 = 0; i1 < gammas->GetEntriesFast(); i1++) {
   	AliAODConversionPhoton * ph1 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i1));
	Int_t tIDs[4] = {ph1->GetLabel(0), ph1->GetLabel(1), -1, -1};

	Int_t leading = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs);
	if(ph1->Pt() > fAxistPt.GetBinLowEdge(1)) {
	  gCorr->CorrelateWithTracks( static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs, leading);
	}
	for(Int_t i2 = 0; i2 < i1; i2++) {
	  AliAODConversionPhoton * ph2 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i2));
	  
	  if( ph2->GetTrackLabelPositive()==ph1->GetTrackLabelPositive() 
		  || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelNegative()
		  || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelPositive()
		  || ph2->GetTrackLabelPositive()==ph1->GetTrackLabelNegative()) {
		continue;
	  }

	  AliAODConversionMother * pion = new AliAODConversionMother(ph1, ph2);
	  pion->SetLabels(i1, i2);
	  
	  if(!fV0Filter || fV0Filter->MesonIsSelected(pion, kTRUE) ) {
	
		Int_t leadingpi = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(pion), tracks, tIDs);
		piCorr->FillTriggerCounters(pion, leadingpi);
		
		tIDs[2] = ph2->GetLabel(0);
		tIDs[3] = ph2->GetLabel(1);
		if(pion->Pt() > fAxistPt.GetBinLowEdge(1) && 
		   pion->M() > fAxisPiM.GetBinLowEdge(1) && 
		   pion->M() < fAxisPiM.GetBinUpEdge(fAxisPiM.GetNbins())) {
		  piCorr->CorrelateWithTracks(pion, tracks, tIDs, leadingpi);
		}
	  }
	}
  }
}


//________________________________________________________________________
void AliAnalysisTaskdPhi::Terminate(Option_t *) {
 
  // Draw result to the screen
  // Called once at the end of the query
}

//________________________________________________________________________
TClonesArray * AliAnalysisTaskdPhi::GetConversionGammas(Bool_t isAOD) {
  

  if(isAOD) {

	TClonesArray * gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
	if(gammas) {
	  return gammas;
	}

	FindDeltaAODBranchName(fInputEvent);
	gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
	return gammas;

  } else {
	TClonesArray * gammas = dynamic_cast<TClonesArray*>(GetInputData(1));
	return gammas;
  }

}

//________________________________________________________________________
void AliAnalysisTaskdPhi::FindDeltaAODBranchName(AliVEvent * event){
  ///Find aod branch
  TList *list=event->GetList();
  for(Int_t ii=0;ii<list->GetEntries();ii++){
	TString name((list->At(ii))->GetName());
	if(name.BeginsWith("GammaConv")&&name.EndsWith("gamma")){
	  fDeltaAODBranchName=name;
	  AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
	  cout <<fDeltaAODBranchName << endl;
	  return;
	}
  }
}
  

///________________________________________________________________________
Bool_t AliAnalysisTaskdPhi::VerifyAODGamma(AliAODConversionPhoton * gamma) {

  AliAODEvent * event = static_cast<AliAODEvent*>(fInputEvent);

  //cout << "label "<< gamma->GetV0Index() << endl;

  AliAODv0 * v0 =  NULL;

  //Int_t v0idx =  gamma->GetV0Index();
  for(Int_t i = 0; i < event->GetNumberOfV0s(); i++) {
	AliAODv0 * tv0 = event->GetV0(i);

	//cout << i << " " << v0->GetID() << " " << v0->GetSecondaryVtx()->GetID() << " " << v0->GetLabel() << " " << v0->GetSecondaryVtx()->GetLabel() << endl;
   
	if(tv0->GetSecondaryVtx()->GetID() == gamma->GetV0Index() ) {
	  v0 = tv0;
	  //cout << "found it" << endl;
	  break;
	}	
  }

  if(!v0) {
	//cout << "v0 not found"<<endl;
	return kFALSE;
  } 


  AliAODTrack * d1 = dynamic_cast<AliAODTrack*>(v0->GetDaughter(0));
  AliAODTrack * d2 = dynamic_cast<AliAODTrack*>(v0->GetDaughter(1));

  Int_t t1 = -1;
  Int_t t2 = -2;

  if(d1) t1 = d1->GetID();
  if(d2) t2 = d2->GetID();

  Int_t g1 = gamma->GetLabel(0);
  Int_t g2 = gamma->GetLabel(1);

  if((t1 == g1 && 
	  t2 == g2) || 
	 (t1 == g2 && 
	  t2 == g1) ) {
	//cout <<"match"<< " " <<  gamma->Pt() << " " <<  d1->Pt() + d2->Pt() <<endl;
  }
		 
  else {

	cout << g1 << " " << g2 <<endl;
	cout << t1 << " " << t2 <<endl;
	
	for(Int_t i = 0; i < event->GetNumberOfV0s(); i++) {
	v0 = event->GetV0(i);
	cout << i << " " << v0->GetSecondaryVtx()->GetID() << " " <<dynamic_cast<AliAODTrack*>(v0->GetDaughter(0))->GetID() << " " << dynamic_cast<AliAODTrack*>(v0->GetDaughter(1))->GetID() << endl; 
	
	}
  }
  
  // Float_t sumdpt = d1->Pt() + d2->Pt();
  // cout << "pt: " << sumdpt << " " << gamma->Pt() << endl;

  return kTRUE;


}
Commit	Line	Data
3b77b2d1	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
2dd6ec26	36	#include "AliConversionTrackCuts.h"
3b77b2d1	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),
3b77b2d1	54	fV0Filter(NULL),
4267afec	55	fTrackCuts(NULL),
3b77b2d1	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),
3b77b2d1	65	fDeltaAODBranchName("AliAODGammaConversion_gamma"),
3bff49c3	66	fAxistPt(),
3bff49c3	67	fAxiscPt(),
3b77b2d1	68	fAxisEta(),
	69	fAxisPhi(),
	70	fAxisCent(),
	71	fAxisZ(),
	72	fAxisPiM()
	73	{
	74	//constructor
3bff49c3	75	fAxistPt.SetNameTitle("tPtAxis", "tPt");
3bff49c3	76	fAxistPt.Set(20, 0, 100);
3b77b2d1	77
3bff49c3	78	fAxiscPt.SetNameTitle("cPtAxis", "cPt");
3bff49c3	79	fAxiscPt.Set(20, 0, 100);
3b77b2d1	80
3b77b2d1	81	fAxisEta.SetNameTitle("EtaAxis", "Eta");
3b77b2d1	82	fAxisEta.Set(160, -0.8, 0.8);
3bff49c3	83
3b77b2d1	84	fAxisPhi.SetNameTitle("PhiAxis", "Phi");
3b77b2d1	85	fAxisPhi.Set(128, 0, TMath::TwoPi());
3bff49c3	86
3b77b2d1	87	fAxisZ.SetNameTitle("ZAxis", "Z");
3b77b2d1	88	fAxisZ.Set(4, -10, 10);
3b77b2d1	89
4267afec	90	fAxisCent.SetNameTitle("CentAxis", "Cent");
3b77b2d1	91	Double_t centbins[5] = {0, 10, 30, 60, 100.1};
	92	fAxisCent.Set(4, centbins);
	93
849102f9	94	Double_t mbins[8] = {0.11, 0.12, 0.13, 0.132, 0.138, 0.14, 0.15, 0.16};
3b77b2d1	95	fAxisPiM.SetNameTitle("InvMassPi0", "Invariant mass");
849102f9	96	fAxisPiM.Set(7, mbins);
3b77b2d1	97
	98	fGammas = new TObjArray();
	99	fGammas->SetOwner(kFALSE);
	100
	101	fPions = new TObjArray();
	102	fPions->SetOwner(kFALSE);
d64241db	103
3b77b2d1	104	// Define input and output slots here
3b77b2d1	105	DefineInput(0, TChain::Class());
43cadb8a	106	//DefineInput(1, TClonesArray::Class());
3b77b2d1	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
4267afec	137	if(fHistoPion)
3b77b2d1	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
43cadb8a	167	for(Int_t ic = 0; ic < fAxisCent.GetNbins(); ic++) {
3b77b2d1	168	TObjArray * photonArray = new TObjArray();
3b77b2d1	169	photonArray->SetOwner(kTRUE);
43cadb8a	170	fPhotonCorr->AddAt(photonArray, ic);
3b77b2d1	171
	172	TObjArray * pionArray = new TObjArray();
	173	pionArray->SetOwner(kTRUE);
43cadb8a	174	fPionCorr->AddAt(pionArray, ic);
3b77b2d1	175
3b77b2d1	176	TList * photonList = new TList();
43cadb8a	177	photonList->SetName(Form("photon_%d", ic));
3b77b2d1	178	photonList->SetOwner(kTRUE);
43cadb8a	179	hPhoton->AddAt(photonList, ic);
3b77b2d1	180
3b77b2d1	181	TList * pionList = new TList();
43cadb8a	182	pionList->SetName(Form("pion_%d", ic));
3b77b2d1	183	pionList->SetOwner(kTRUE);
43cadb8a	184	hPion->AddAt(pionList, ic);
3b77b2d1	185
3b77b2d1	186
43cadb8a	187	for(Int_t iz = 0; iz < fAxisZ.GetNbins(); iz++) {
43cadb8a	188	TString nameString = Form("%d_%d", ic, iz);
3b77b2d1	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()));
43cadb8a	196	photonArray->AddAt(photonCorr, iz);
3bff49c3	197	photonCorr->GetAxistPt().Set(fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
3bff49c3	198	photonCorr->GetAxiscPt().Set(fAxiscPt.GetNbins(), fAxiscPt.GetXbins()->GetArray());
3b77b2d1	199	photonCorr->CreateHistograms();
	200	photonList->Add(photonCorr->GetHistograms());
	201
	202	AliAnaConvCorrPion * pionCorr = new AliAnaConvCorrPion(Form("PionCorr_%s", nameString.Data()), Form("pion %s", titleString.Data()));
43cadb8a	203	pionArray->AddAt(pionCorr, iz);
3bff49c3	204	pionCorr->GetAxistPt().Set(fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
3bff49c3	205	pionCorr->GetAxiscPt().Set(fAxiscPt.GetNbins(), fAxiscPt.GetXbins()->GetArray());
3b77b2d1	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
4267afec	236	AliConversionTrackCuts * tc = dynamic_cast<AliConversionTrackCuts*>(fTrackCuts);
4267afec	237	if(tc) fHistograms->Add(tc->CreateHistograms());
5c76942e	238
3b77b2d1	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
3b77b2d1	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
3b77b2d1	255	TList * outAxesList = new TList();
	256	outAxesList->Add(&fAxisCent);
	257	outAxesList->Add(&fAxisZ);
43cadb8a	258	fHistograms->Add(outAxesList);
3b77b2d1	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) {
43cadb8a	268	///Create sparse
3b77b2d1	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;
3b77b2d1	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 eventPlane = 0.0;
	338	Double_t vertexz = fInputEvent->GetPrimaryVertex()->GetZ();
	339	if(isAOD) {
	340	AliAODHeader * header = static_cast<AliAODHeader*>(fInputEvent->GetHeader());
	341	centrality = header->GetCentrality();
	342	eventPlane = header->GetEventplane();
	343	} else {
4267afec	344	centrality = static_cast<AliESDEvent*>(fInputEvent)->GetCentrality()->GetCentralityPercentile("V0M");
3b77b2d1	345	eventPlane = fInputEvent->GetEventplane()->GetEventplane("Q");
	346	}
	347
4267afec	348
	349	const Int_t centBin = GetBin(fAxisCent, centrality);
	350	const Int_t vertexBin = GetBin(fAxisZ, vertexz);
	351
	352
3bff49c3	353	if(DebugLevel () > 4) {
3b77b2d1	354	cout << "centrality: " << centrality << " " << GetBin(fAxisCent, centrality) << endl;
	355	cout << "vertexz: " << vertexz << " " << GetBin(fAxisZ, vertexz) << endl;
	356	cout << "eventPlane: " << eventPlane << " " << endl;
	357	}
	358
3b77b2d1	359
	360	if(centBin < 0 \|\| vertexBin < 0) {
	361	AliError("bin out of range");
849102f9	362	//cout << "bad bin"<<endl;
3b77b2d1	363	return;
	364	}
	365
	366	fGammas->Clear();
	367	fPions->Clear();
	368
	369	TClonesArray * aodGammas = GetConversionGammas(isAOD);
a60bc693	370
3b77b2d1	371	if(!aodGammas) {
	372	AliError("no aod gammas found!");
	373	return;
	374	}
	375
849102f9	376
3b77b2d1	377	if(DebugLevel() > 1) printf("Number of conversion gammas %d \n", aodGammas->GetEntriesFast());
3b77b2d1	378	for(Int_t ig = 0; ig < aodGammas->GetEntriesFast(); ig++) {
849102f9	379	AliAODConversionPhoton * photon = dynamic_cast<AliAODConversionPhoton*>(aodGammas->At(ig));
3b77b2d1	380
a60bc693	381	if(!photon) {
	382	cout << "can't get photon"<<endl;
	383	continue;
	384	}
849102f9	385	if(!VerifyAODGamma(photon)) {
849102f9	386	continue;
a60bc693	387	}
a60bc693	388
3b77b2d1	389	if(!fV0Filter \|\| fV0Filter->PhotonIsSelected(static_cast<AliConversionPhotonBase*>(photon), fInputEvent)) {
	390	fGammas->Add(static_cast<TObject*>(photon));
	391	}
	392	}
	393
	394	if(DebugLevel() > 4) printf("Number of accepted gammas %d \n", fGammas->GetEntriesFast());
4267afec	395	hMEvents->Fill(vertexz, centrality);
4267afec	396
3b77b2d1	397
4267afec	398
4267afec	399	///create track array
3b77b2d1	400	TObjArray tracks;
4267afec	401	const Double_t etalim[2] = { fAxisEta.GetBinLowEdge(1), fAxisEta.GetBinUpEdge(fAxisEta.GetNbins())};
3b77b2d1	402	for(Int_t iTrack = 0; iTrack < fInputEvent->GetNumberOfTracks(); iTrack++) {
	403
	404	AliVTrack * track = static_cast<AliVTrack*>(fInputEvent->GetTrack(iTrack));
4267afec	405	if(track->Pt() < fAxiscPt.GetBinLowEdge(1) ) continue;
	406	if(track->Eta() < etalim[0] \|\| track->Eta() > etalim[1]) continue;
	407
	408
	409	if(!fTrackCuts \|\| fTrackCuts->IsSelected((track))) {
2dd6ec26	410	tracks.Add(track);
2dd6ec26	411	}
3b77b2d1	412	}
	413
	414	Process(fGammas, &tracks, vertexBin, centBin);
	415
	416	PostData(1, fHistograms);
	417	PostData(2, fHistoGamma);
	418	PostData(3, fHistoPion);
	419
	420	}
	421
	422
	423	//________________________________________________________________________
	424	void AliAnalysisTaskdPhi::Process(TObjArray * gammas, TObjArray * tracks, Int_t vertexBin, Int_t centBin) {
	425	///Process stuff
	426
4267afec	427	if(DebugLevel() > 4) printf("Number of accepted gammas, tracks %d %d \n", gammas->GetEntriesFast(), tracks->GetEntriesFast());
849102f9	428
3b77b2d1	429
3b77b2d1	430	AliAnaConvCorrBase * gCorr = GetCorrObject(vertexBin, centBin, fPhotonCorr);
4267afec	431	AliAnaConvCorrPion * piCorr = dynamic_cast<AliAnaConvCorrPion*>(GetCorrObject(vertexBin, centBin, fPionCorr));
3b77b2d1	432
	433	if(!gCorr \|\| !piCorr) {
	434	AliError("corr object missing");
	435	return;
	436	}
	437
3b77b2d1	438	for(Int_t i1 = 0; i1 < gammas->GetEntriesFast(); i1++) {
	439	AliAODConversionPhoton * ph1 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i1));
	440	Int_t tIDs[4] = {ph1->GetLabel(0), ph1->GetLabel(1), -1, -1};
	441
3bff49c3	442	Int_t leading = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs);
e3598cd0	443	if(ph1->Pt() > fAxistPt.GetBinLowEdge(1)) {
	444	gCorr->CorrelateWithTracks( static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs, leading);
	445	}
3b77b2d1	446	for(Int_t i2 = 0; i2 < i1; i2++) {
3b77b2d1	447	AliAODConversionPhoton * ph2 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i2));
3bff49c3	448
3b77b2d1	449	if( ph2->GetTrackLabelPositive()==ph1->GetTrackLabelPositive()
	450	\|\| ph2->GetTrackLabelNegative()==ph1->GetTrackLabelNegative()
	451	\|\| ph2->GetTrackLabelNegative()==ph1->GetTrackLabelPositive()
	452	\|\| ph2->GetTrackLabelPositive()==ph1->GetTrackLabelNegative()) {
	453	continue;
	454	}
	455
	456	AliAODConversionMother * pion = new AliAODConversionMother(ph1, ph2);
	457	pion->SetLabels(i1, i2);
3bff49c3	458
3b77b2d1	459	if(!fV0Filter \|\| fV0Filter->MesonIsSelected(pion, kTRUE) ) {
4267afec	460
3bff49c3	461	Int_t leadingpi = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(pion), tracks, tIDs);
d64241db	462	piCorr->FillTriggerCounters(pion, leadingpi);
849102f9	463
3bff49c3	464	tIDs[2] = ph2->GetLabel(0);
3bff49c3	465	tIDs[3] = ph2->GetLabel(1);
d64241db	466	if(pion->Pt() > fAxistPt.GetBinLowEdge(1) &&
	467	pion->M() > fAxisPiM.GetBinLowEdge(1) &&
	468	pion->M() < fAxisPiM.GetBinUpEdge(fAxisPiM.GetNbins())) {
e3598cd0	469	piCorr->CorrelateWithTracks(pion, tracks, tIDs, leadingpi);
e3598cd0	470	}
3b77b2d1	471	}
	472	}
	473	}
	474	}
	475
849102f9	476
3b77b2d1	477	//________________________________________________________________________
	478	void AliAnalysisTaskdPhi::Terminate(Option_t *) {
	479
	480	// Draw result to the screen
	481	// Called once at the end of the query
	482	}
	483
	484	//________________________________________________________________________
	485	TClonesArray * AliAnalysisTaskdPhi::GetConversionGammas(Bool_t isAOD) {
	486
	487
	488	if(isAOD) {
	489
	490	TClonesArray * gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
	491	if(gammas) {
	492	return gammas;
	493	}
e3598cd0	494
	495	FindDeltaAODBranchName(fInputEvent);
	496	gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
3b77b2d1	497	return gammas;
	498
	499	} else {
	500	TClonesArray * gammas = dynamic_cast<TClonesArray*>(GetInputData(1));
	501	return gammas;
	502	}
	503
	504	}
	505
	506	//________________________________________________________________________
e3598cd0	507	void AliAnalysisTaskdPhi::FindDeltaAODBranchName(AliVEvent * event){
3b77b2d1	508	///Find aod branch
	509	TList *list=event->GetList();
	510	for(Int_t ii=0;ii<list->GetEntries();ii++){
	511	TString name((list->At(ii))->GetName());
	512	if(name.BeginsWith("GammaConv")&&name.EndsWith("gamma")){
	513	fDeltaAODBranchName=name;
	514	AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
a60bc693	515	cout <<fDeltaAODBranchName << endl;
3b77b2d1	516	return;
	517	}
	518	}
	519	}
	520
	521
a60bc693	522
	523
	524	///________________________________________________________________________
	525	Bool_t AliAnalysisTaskdPhi::VerifyAODGamma(AliAODConversionPhoton * gamma) {
	526
	527	AliAODEvent * event = static_cast<AliAODEvent*>(fInputEvent);
	528
849102f9	529	//cout << "label "<< gamma->GetV0Index() << endl;
a60bc693	530
	531	AliAODv0 * v0 = NULL;
	532
	533	//Int_t v0idx = gamma->GetV0Index();
	534	for(Int_t i = 0; i < event->GetNumberOfV0s(); i++) {
	535	AliAODv0 * tv0 = event->GetV0(i);
	536
	537	//cout << i << " " << v0->GetID() << " " << v0->GetSecondaryVtx()->GetID() << " " << v0->GetLabel() << " " << v0->GetSecondaryVtx()->GetLabel() << endl;
	538
	539	if(tv0->GetSecondaryVtx()->GetID() == gamma->GetV0Index() ) {
	540	v0 = tv0;
	541	//cout << "found it" << endl;
	542	break;
	543	}
	544	}
	545
	546	if(!v0) {
849102f9	547	//cout << "v0 not found"<<endl;
a60bc693	548	return kFALSE;
	549	}
	550
	551
	552	AliAODTrack * d1 = dynamic_cast<AliAODTrack*>(v0->GetDaughter(0));
	553	AliAODTrack * d2 = dynamic_cast<AliAODTrack*>(v0->GetDaughter(1));
	554
	555	Int_t t1 = -1;
	556	Int_t t2 = -2;
	557
	558	if(d1) t1 = d1->GetID();
	559	if(d2) t2 = d2->GetID();
	560
	561	Int_t g1 = gamma->GetLabel(0);
	562	Int_t g2 = gamma->GetLabel(1);
	563
	564	if((t1 == g1 &&
	565	t2 == g2) \|\|
	566	(t1 == g2 &&
	567	t2 == g1) ) {
849102f9	568	//cout <<"match"<< " " << gamma->Pt() << " " << d1->Pt() + d2->Pt() <<endl;
a60bc693	569	}
	570
	571	else {
	572
	573	cout << g1 << " " << g2 <<endl;
	574	cout << t1 << " " << t2 <<endl;
	575
	576	for(Int_t i = 0; i < event->GetNumberOfV0s(); i++) {
	577	v0 = event->GetV0(i);
	578	cout << i << " " << v0->GetSecondaryVtx()->GetID() << " " <<dynamic_cast<AliAODTrack>(v0->GetDaughter(0))->GetID() << " " << dynamic_cast<AliAODTrack>(v0->GetDaughter(1))->GetID() << endl;
	579
	580	}
	581	}
	582
	583	// Float_t sumdpt = d1->Pt() + d2->Pt();
	584	// cout << "pt: " << sumdpt << " " << gamma->Pt() << endl;
	585
	586	return kTRUE;
	587
	588
	589
	590	}