[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / ClusterSelection / AliPHOSClusterSelection.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * 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.                  *
 **************************************************************************/

/* $Id$ */

#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliMultiInputEventHandler.h"
#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliESDCaloCluster.h"
#include "AliVParticle.h"
#include "AliVTrack.h"
#include "AliAODCluster.h"
#include "AliVCaloCells.h"
#include "AliOADBContainer.h"
#include "AliPHOSGeometry.h"
#include "AliAnalysisManager.h"


#include "AliPHOSClusterSelection.h"

AliPHOSClusterSelection::AliPHOSClusterSelection()
  : fMinChargedParticleTrackDistance(-1.), 
    fNotUnfolded(false),
    fMaxDispR2(-1.),
    fIsCore(-1.),
    fMaxTOF(-1.)
{
  // Defaults to the most lenient selection allowable
	return;
}

AliPHOSClusterSelection::~AliPHOSClusterSelection()
{
}

Bool_t AliPHOSClusterSelection::IsSelected(AliVCluster* cluster) const
{
  return IsSelectedCPV(cluster)
    && IsSelectedUnfolded(cluster)
    && IsSelectedDisp(cluster)
    && IsSelectedTOF(cluster);
}

Bool_t AliPHOSClusterSelection::IsSelectedCPV(AliVCluster* cluster) const
{
  if( 0 > fMinChargedParticleTrackDistance )//No selection on CPV
    return true; 

  else{
    Bool_t cpvBit = true; //Changed to false if there is a matching track in the requested radius.

    Double_t dx=cluster->GetTrackDx();//Position for Track matching
    Double_t dz=cluster->GetTrackDz();


    AliESDEvent* EventESD = dynamic_cast<AliESDEvent*> (AliPHOSClusterSelection::GetCurrentEvent());//dynamic cast to test for ESD or AOD
    AliAODEvent* EventAOD = dynamic_cast<AliAODEvent*> (AliPHOSClusterSelection::GetCurrentEvent());
    Double_t mf = 0.; //
    if(EventAOD) mf = EventAOD->GetMagneticField(); //Positive for ++ and negative for -- 
    else if(EventESD) mf = EventESD->GetMagneticField(); //Positive for ++ and negative for --


    if(EventESD){//if ESD
      AliESDCaloCluster * ESDcluster = static_cast<AliESDCaloCluster*> (cluster);//Know its ESD so static cast is fine?
      TArrayI * itracks = ESDcluster-> GetTracksMatched() ;
      if(itracks->GetSize()>0){
    	Int_t iTr = itracks->At(0);
    	if(iTr>=0 && iTr<EventESD->GetNumberOfTracks()){
    	  AliVParticle* track = EventESD->GetTrack(iTr);
    	  Double_t pt = track->Pt() ;
    	  Short_t charge = track->Charge() ;
    	  Double_t r = AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf);
    	  cpvBit=(r>fMinChargedParticleTrackDistance) ;
    	}
      }
    }


    if(EventAOD){//if AOD
      AliAODCluster * AODcluster = static_cast<AliAODCluster*> (cluster);
      int nTracksMatched = AODcluster->GetNTracksMatched();
      if(nTracksMatched > 0) {
	AliVTrack* track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
	if ( track ) {
	  Double_t pt = track->Pt();
	  Short_t charge = track->Charge();
	  Double_t r = AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf);
	  cpvBit=(r>fMinChargedParticleTrackDistance) ;
	}
      }
    }
    return cpvBit;
  }
}


Bool_t AliPHOSClusterSelection::IsSelectedUnfolded(AliVCluster* cluster) const
{
  if(!fNotUnfolded)
    return true;
  else{
    Bool_t NotUnfolded = cluster->GetNExMax()<2;//True if it was not unfolded
    return !NotUnfolded;
  }
}

Bool_t AliPHOSClusterSelection::IsSelectedDisp(AliVCluster* cluster) const
{
  if(0 > fMaxDispR2)
    return true;
  else{
    Double_t m02 = 0.,m20 = 0.;
    if(fIsCore<0){//No core calculation
      m02 = cluster->GetM02();
      m20 = cluster->GetM20();
      return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02) ;
    }
    else{//DispCore
      AliVCaloCells* cells = static_cast<AliVCaloCells*> (AliPHOSClusterSelection::GetCurrentEvent()->GetPHOSCells());//Need the cells
      AliPHOSClusterSelection::EvalCoreLambdas(cluster, cells, m02, m20);
      return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02);
    }
  }

}

Bool_t AliPHOSClusterSelection::IsSelectedTOF(AliVCluster* cluster) const
{
  if(0 > fMaxTOF)
    return true;
  else{ 
    // Time of Flight (TOF)
    Double_t tof = cluster->GetTOF();//Time of Flight for the cluster
    return  TMath::Abs(tof) < fMaxTOF;//True if the cut is passed
  }
}

AliPHOSClusterSelection* AliPHOSClusterSelection::SetMinChargedParticleTrackDistance(Float_t distance)
{
  // 'distance' set the minimal allowable distance between the cluster
  // and the nearest extrapolated track.
  // If 'distance' is negative, then all clusters are sellected, the selection
  // being "not applied" or "disabled".
  
  fMinChargedParticleTrackDistance = distance;
  return this;
}

AliPHOSClusterSelection* AliPHOSClusterSelection::SetNotUnfolded(Bool_t notUnfolded)
{
  //if notUnfolded true, it rejects Unfolded Clusters
  fNotUnfolded = notUnfolded;
  return this;
}

AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxDispR2(Float_t maxR2)
{
  // 'maxR2' sets the maximum allowed dispersion.
  // If 'maxR2' is negative, then all clusters are selected, the selection
  // being "not applied" or "disabled".
  
  fMaxDispR2 = maxR2;
  return this;
}

AliPHOSClusterSelection* AliPHOSClusterSelection::SetIsCore(Bool_t isCore)
{
  // 'isCore' sets wether core version of Disp is used. -1 gives no core.
  
  fIsCore = isCore;
  return this;
}

AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxTOF(Float_t maxTOF)
{
  // 'maxTOF' sets the maximum allowed time of flight for the cluster.
  // If 'maxTOF' is negative, all clusters are selected and the selection is "disabled".
  
  fMaxTOF = maxTOF;
  return this;
}

TString AliPHOSClusterSelection::ToString() const
{
  // returns a string an quasi-unique string for whatever selection 
  // parameters the instance contains. The uniqueness of the string
  // is limited by the precision given in the formatting of the string.
  // Take care that the precision is sufficient for your needs.

  return TString::Format("%.1f_%i_%.1f_%.1f_%.1f",
			 fMinChargedParticleTrackDistance,
			 fNotUnfolded,
			 fMaxDispR2,
			 fCoreRadius,
			 fMaxTOF
			 );
}


Float_t AliPHOSClusterSelection::GetMinChargedParticleTrackDistance(const TString& string)
{
  TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(0))->String());
  Float_t flt= s->Atof();
  //delete s; //s is defined in the stack, no delete necessary
  return flt;
}

Bool_t AliPHOSClusterSelection::GetUnfolded(const TString& string)
{
  TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(1))->String());
  Bool_t bl = s->Atoi(); 
  return bl;
}

Float_t AliPHOSClusterSelection::GetMaxDispR2(const TString& string)
{
  TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(2))->String());
  Float_t flt = s->Atof(); 
  return flt;
}

Bool_t AliPHOSClusterSelection::GetIsCore(const TString& string)
{
  TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(3))->String());
  Bool_t blt = s->Atoi(); 
  return blt;
}

Float_t AliPHOSClusterSelection::GetMaxTOF(const TString& string)
{
  TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(4))->String());
  Float_t flt = s->Atof(); 
  return flt;
}


Double_t AliPHOSClusterSelection::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, Double_t mf) const {
  //Parameterization of LHC10h period
  //_true if neutral_
  //Copied from Pi0Flow task

  Double_t meanX=0;
  Double_t meanZ=0.;
  Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
			 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
  Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
 
  if(mf<0.){ //field --
    meanZ = -0.468318 ;
    if(charge>0)
      meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
    else
      meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
  }
  else{ //Field ++
    meanZ= -0.468318;
    if(charge>0)
      meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
    else
      meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
  }

  Double_t rz=(dz-meanZ)/sz ;
  Double_t rx=(dx-meanX)/sx ;
  return TMath::Sqrt(rx*rx+rz*rz) ;
}

//____________________________________________________________________________
void  AliPHOSClusterSelection::EvalCoreLambdas(AliVCluster  * clu, AliVCaloCells * cells,Double_t &m02, Double_t &m20) const
{ 
  //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
  //Copied from pi0flowtask
  Int_t fRunNumber;
  const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;  
  // Calculates the center of gravity in the local PHOS-module coordinates
  Float_t wtot = 0;
  const Int_t mulDigit=clu->GetNCells() ;
  Double_t xc[mulDigit] ;
  Double_t zc[mulDigit] ;
  Double_t wi[mulDigit] ;
  Double_t x = 0 ;
  Double_t z = 0 ;
  const Double_t logWeight=4.5 ;
  for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
    Int_t relid[4] ;
    Float_t xi=0. ;
    Float_t zi=0. ;
    Int_t absId = clu->GetCellAbsId(iDigit) ;


    AliESDEvent* EventESD = dynamic_cast<AliESDEvent*> (AliPHOSClusterSelection::GetCurrentEvent());//dynamic cast to test for ESD or AOD

    AliAODEvent* EventAOD = dynamic_cast<AliAODEvent*> (AliPHOSClusterSelection::GetCurrentEvent());
    if(EventESD)
      fRunNumber = EventESD->GetRunNumber() ;
    if(EventAOD)
      fRunNumber = EventESD->GetRunNumber() ;

    AliOADBContainer geomContainer("phosGeo");//Initialize Geometry
    geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
    TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
    AliPHOSGeometry * fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
    for(Int_t mod=0; mod<5; mod++) {
      if(!matrixes->At(mod)) {
	AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
	continue;
      }


      fPHOSGeo->AbsToRelNumbering(absId, relid) ;
      fPHOSGeo->RelPosInModule(relid, xi, zi);
      xc[iDigit]=xi ;
      zc[iDigit]=zi ;
      Double_t ei = elist[iDigit]*cells->GetCellAmplitude(absId) ;
      wi[iDigit]=0. ;
      if (clu->E()>0 && ei>0) {
	wi[iDigit] = TMath::Max( 0., logWeight + TMath::Log( ei / clu->E() ) ) ;
	Double_t w=wi[iDigit];
	x    += xc[iDigit] * w ;
	z    += zc[iDigit] * w ;
	wtot += w ;
      }
    }
    if (wtot>0) {
      x /= wtot ;
      z /= wtot ;
    }
     
    wtot = 0. ;
    Double_t dxx  = 0.;
    Double_t dzz  = 0.;
    Double_t dxz  = 0.;
    Double_t xCut = 0. ;
    Double_t zCut = 0. ;
    for(Int_t iDigit1=0; iDigit1<mulDigit; iDigit1++) {
      Double_t w=wi[iDigit1];
      if (w>0.) {
        Double_t xi1= xc[iDigit1] ;
        Double_t zi1= zc[iDigit1] ;
	if((xi1-x)*(xi1-x)+(zi1-z)*(zi1-z) < 4.5*4.5){
          xCut += w * xi1 ;
          zCut += w * zi1 ; 
          dxx  += w * xi1 * xi1 ;
          dzz  += w * zi1 * zi1 ;
          dxz  += w * xi1 * zi1 ; 
          wtot += w ;
	}
      }
    
    }
    if (wtot>0) {
      xCut/= wtot ;
      zCut/= wtot ;
      dxx /= wtot ;
      dzz /= wtot ;
      dxz /= wtot ;
      dxx -= xCut * xCut ;
      dzz -= zCut * zCut ;
      dxz -= xCut * zCut ;

      m02 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
      m20 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
    }
    else {
      m20=m02=0.;
    }

  }
}

//_____________________________________________________________________________
Bool_t AliPHOSClusterSelection::TestLambda(Double_t pt,Double_t l1,Double_t l2) const 
{
  //Evaluates if lambdas correspond to photon cluster
  //Tuned using pp data 
  //copied from Pi0FlowTask
  Double_t l2Mean, l1Mean, l2Sigma, l1Sigma, c, R2;
  if(fIsCore<0){
    l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
    l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
    l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
    l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
    c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
    R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
      0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
      0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
  }
  else{
    //For core radius R=4.5
    l1Mean  = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ;
    l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ;
    l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ;
    l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ;
    c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ;
    R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma + 
      0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
      0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
  }
  return (R2<fMaxDispR2*fMaxDispR2);//TODO: Check if this should be 2. order
}


AliVEvent* AliPHOSClusterSelection::GetCurrentEvent() const
{
  // Hackish way of getting the current event.
  // Its probably not appropriate to call this function outside of
  // AliAnalysisTaskSE::UserExec
  
  AliAnalysisManager* analysisManager = dynamic_cast<AliAnalysisManager*>(AliAnalysisManager::GetAnalysisManager());
  AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(analysisManager->GetInputEventHandler());
  AliMultiInputEventHandler *multiInputHandler = dynamic_cast<AliMultiInputEventHandler *>(inputHandler);
  if (multiInputHandler)
    inputHandler = dynamic_cast<AliInputEventHandler *>(multiInputHandler->GetFirstInputEventHandler());
  
  AliVEvent* inputEvent = dynamic_cast<AliVEvent*>(inputHandler->GetEvent());
  if( ! inputEvent ) 
    AliError("Was not able to retrieve event!");
  
  return inputEvent;
}
Commit	Line	Data
cc909f95 HQ	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	#include "AliAnalysisManager.h"
	19	#include "AliInputEventHandler.h"
	20	#include "AliMultiInputEventHandler.h"
	21	#include "AliVEvent.h"
17ba2205	22	#include "AliESDEvent.h"
	23	#include "AliAODEvent.h"
	24	#include "AliESDCaloCluster.h"
	25	#include "AliVParticle.h"
	26	#include "AliVTrack.h"
	27	#include "AliAODCluster.h"
	28	#include "AliVCaloCells.h"
	29	#include "AliOADBContainer.h"
	30	#include "AliPHOSGeometry.h"
	31	#include "AliAnalysisManager.h"
	32
	33
cc909f95	34
9607c225 HQ	35	#include "AliPHOSClusterSelection.h"
	36
	37	AliPHOSClusterSelection::AliPHOSClusterSelection()
	38	: fMinChargedParticleTrackDistance(-1.),
	39	fNotUnfolded(false),
	40	fMaxDispR2(-1.),
ba4942ff	41	fIsCore(-1.),
9607c225 HQ	42	fMaxTOF(-1.)
	43	{
	44	// Defaults to the most lenient selection allowable
	45	return;
	46	}
	47
	48	AliPHOSClusterSelection::~AliPHOSClusterSelection()
	49	{
	50	}
	51
	52	Bool_t AliPHOSClusterSelection::IsSelected(AliVCluster* cluster) const
	53	{
	54	return IsSelectedCPV(cluster)
	55	&& IsSelectedUnfolded(cluster)
	56	&& IsSelectedDisp(cluster)
9607c225 HQ	57	&& IsSelectedTOF(cluster);
	58	}
	59
ba4942ff	60	Bool_t AliPHOSClusterSelection::IsSelectedCPV(AliVCluster* cluster) const
9607c225	61	{
17ba2205	62	if( 0 > fMinChargedParticleTrackDistance )//No selection on CPV
9607c225	63	return true;
17ba2205	64
	65	else{
	66	Bool_t cpvBit = true; //Changed to false if there is a matching track in the requested radius.
	67
	68	Double_t dx=cluster->GetTrackDx();//Position for Track matching
	69	Double_t dz=cluster->GetTrackDz();
	70
	71
	72	AliESDEvent* EventESD = dynamic_cast<AliESDEvent*> (AliPHOSClusterSelection::GetCurrentEvent());//dynamic cast to test for ESD or AOD
	73	AliAODEvent* EventAOD = dynamic_cast<AliAODEvent*> (AliPHOSClusterSelection::GetCurrentEvent());
	74	Double_t mf = 0.; //
	75	if(EventAOD) mf = EventAOD->GetMagneticField(); //Positive for ++ and negative for --
	76	else if(EventESD) mf = EventESD->GetMagneticField(); //Positive for ++ and negative for --
	77
	78
	79	if(EventESD){//if ESD
	80	AliESDCaloCluster * ESDcluster = static_cast<AliESDCaloCluster*> (cluster);//Know its ESD so static cast is fine?
	81	TArrayI * itracks = ESDcluster-> GetTracksMatched() ;
	82	if(itracks->GetSize()>0){
	83	Int_t iTr = itracks->At(0);
ba4942ff	84	if(iTr>=0 && iTr<EventESD->GetNumberOfTracks()){
ba4942ff	85	AliVParticle* track = EventESD->GetTrack(iTr);
17ba2205	86	Double_t pt = track->Pt() ;
17ba2205	87	Short_t charge = track->Charge() ;
ba4942ff	88	Double_t r = AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf);
17ba2205	89	cpvBit=(r>fMinChargedParticleTrackDistance) ;
	90	}
	91	}
	92	}
	93
	94
	95	if(EventAOD){//if AOD
	96	AliAODCluster * AODcluster = static_cast<AliAODCluster*> (cluster);
	97	int nTracksMatched = AODcluster->GetNTracksMatched();
	98	if(nTracksMatched > 0) {
	99	AliVTrack* track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
	100	if ( track ) {
	101	Double_t pt = track->Pt();
	102	Short_t charge = track->Charge();
ba4942ff	103	Double_t r = AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf);
17ba2205	104	cpvBit=(r>fMinChargedParticleTrackDistance) ;
	105	}
	106	}
	107	}
	108	return cpvBit;
	109	}
	110	}
	111
	112
ba4942ff	113	Bool_t AliPHOSClusterSelection::IsSelectedUnfolded(AliVCluster* cluster) const
17ba2205	114	{
	115	if(!fNotUnfolded)
	116	return true;
	117	else{
	118	Bool_t NotUnfolded = cluster->GetNExMax()<2;//True if it was not unfolded
	119	return !NotUnfolded;
	120	}
	121	}
	122
ba4942ff	123	Bool_t AliPHOSClusterSelection::IsSelectedDisp(AliVCluster* cluster) const
17ba2205	124	{
	125	if(0 > fMaxDispR2)
	126	return true;
	127	else{
	128	Double_t m02 = 0.,m20 = 0.;
ba4942ff	129	if(fIsCore<0){//No core calculation
17ba2205	130	m02 = cluster->GetM02();
	131	m20 = cluster->GetM20();
	132	return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02) ;
	133	}
ba4942ff	134	else{//DispCore
17ba2205	135	AliVCaloCells* cells = static_cast<AliVCaloCells*> (AliPHOSClusterSelection::GetCurrentEvent()->GetPHOSCells());//Need the cells
	136	AliPHOSClusterSelection::EvalCoreLambdas(cluster, cells, m02, m20);
	137	return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02);
	138	}
	139	}
	140
	141	}
	142
ba4942ff	143	Bool_t AliPHOSClusterSelection::IsSelectedTOF(AliVCluster* cluster) const
17ba2205	144	{
	145	if(0 > fMaxTOF)
	146	return true;
	147	else{
	148	// Time of Flight (TOF)
	149	Double_t tof = cluster->GetTOF();//Time of Flight for the cluster
	150	return TMath::Abs(tof) < fMaxTOF;//True if the cut is passed
	151	}
9607c225 HQ	152	}
	153
	154	AliPHOSClusterSelection* AliPHOSClusterSelection::SetMinChargedParticleTrackDistance(Float_t distance)
	155	{
	156	// 'distance' set the minimal allowable distance between the cluster
	157	// and the nearest extrapolated track.
	158	// If 'distance' is negative, then all clusters are sellected, the selection
	159	// being "not applied" or "disabled".
	160
	161	fMinChargedParticleTrackDistance = distance;
17ba2205	162	return this;
	163	}
	164
	165	AliPHOSClusterSelection* AliPHOSClusterSelection::SetNotUnfolded(Bool_t notUnfolded)
	166	{
	167	//if notUnfolded true, it rejects Unfolded Clusters
	168	fNotUnfolded = notUnfolded;
	169	return this;
	170	}
	171
	172	AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxDispR2(Float_t maxR2)
	173	{
	174	// 'maxR2' sets the maximum allowed dispersion.
	175	// If 'maxR2' is negative, then all clusters are selected, the selection
	176	// being "not applied" or "disabled".
	177
	178	fMaxDispR2 = maxR2;
	179	return this;
	180	}
	181
ba4942ff	182	AliPHOSClusterSelection* AliPHOSClusterSelection::SetIsCore(Bool_t isCore)
17ba2205	183	{
ba4942ff	184	// 'isCore' sets wether core version of Disp is used. -1 gives no core.
17ba2205	185
ba4942ff	186	fIsCore = isCore;
17ba2205	187	return this;
	188	}
	189
	190	AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxTOF(Float_t maxTOF)
	191	{
	192	// 'maxTOF' sets the maximum allowed time of flight for the cluster.
	193	// If 'maxTOF' is negative, all clusters are selected and the selection is "disabled".
	194
ba4942ff	195	fMaxTOF = maxTOF;
17ba2205	196	return this;
9607c225 HQ	197	}
	198
	199	TString AliPHOSClusterSelection::ToString() const
	200	{
	201	// returns a string an quasi-unique string for whatever selection
	202	// parameters the instance contains. The uniqueness of the string
	203	// is limited by the precision given in the formatting of the string.
	204	// Take care that the precision is sufficient for your needs.
	205
17ba2205	206	return TString::Format("%.1f_%i_%.1f_%.1f_%.1f",
9607c225 HQ	207	fMinChargedParticleTrackDistance,
	208	fNotUnfolded,
	209	fMaxDispR2,
17ba2205	210	fCoreRadius,
9607c225 HQ	211	fMaxTOF
	212	);
	213	}
	214
	215
17ba2205	216	Float_t AliPHOSClusterSelection::GetMinChargedParticleTrackDistance(const TString& string)
17ba2205	217	{
ba4942ff	218	TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(0))->String());
	219	Float_t flt= s->Atof();
	220	//delete s; //s is defined in the stack, no delete necessary
17ba2205	221	return flt;
	222	}
	223
	224	Bool_t AliPHOSClusterSelection::GetUnfolded(const TString& string)
	225	{
ba4942ff	226	TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(1))->String());
ba4942ff	227	Bool_t bl = s->Atoi();
17ba2205	228	return bl;
	229	}
	230
	231	Float_t AliPHOSClusterSelection::GetMaxDispR2(const TString& string)
	232	{
ba4942ff	233	TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(2))->String());
ba4942ff	234	Float_t flt = s->Atof();
17ba2205	235	return flt;
	236	}
	237
ba4942ff	238	Bool_t AliPHOSClusterSelection::GetIsCore(const TString& string)
17ba2205	239	{
ba4942ff	240	TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(3))->String());
	241	Bool_t blt = s->Atoi();
	242	return blt;
17ba2205	243	}
	244
	245	Float_t AliPHOSClusterSelection::GetMaxTOF(const TString& string)
9607c225	246	{
ba4942ff	247	TString * s =&(static_cast<TObjString*>(string.Tokenize("_")->At(4))->String());
ba4942ff	248	Float_t flt = s->Atof();
9607c225 HQ	249	return flt;
9607c225 HQ	250	}
cc909f95 HQ	251
cc909f95 HQ	252
ba4942ff	253	Double_t AliPHOSClusterSelection::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, Double_t mf) const {
17ba2205	254	//Parameterization of LHC10h period
	255	//_true if neutral_
	256	//Copied from Pi0Flow task
	257
	258	Double_t meanX=0;
	259	Double_t meanZ=0.;
	260	Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
	261	6.58365e-015.91917e-015.91917e-01/((pt-9.61306e-01)(pt-9.61306e-01)+5.91917e-015.91917e-01)+1.59219);
	262	Double_t sz=TMath::Min(2.75,4.90341e+021.91456e-021.91456e-02/(ptpt+1.91456e-021.91456e-02)+1.60) ;
	263
	264	if(mf<0.){ //field --
	265	meanZ = -0.468318 ;
	266	if(charge>0)
	267	meanX=TMath::Min(7.3, 3.899941.206791.20679/(ptpt+1.206791.20679)+0.249029+2.49088e+07TMath::Exp(-pt3.33650e+01)) ;
	268	else
	269	meanX=-TMath::Min(7.7,3.860400.9124990.912499/(ptpt+0.9124990.912499)+1.23114+4.48277e+05TMath::Exp(-pt2.57070e+01)) ;
	270	}
	271	else{ //Field ++
	272	meanZ= -0.468318;
	273	if(charge>0)
	274	meanX=-TMath::Min(8.0,3.860401.313571.31357/(ptpt+1.313571.31357)+0.880579+7.56199e+06TMath::Exp(-pt3.08451e+01)) ;
	275	else
	276	meanX= TMath::Min(6.85, 3.899941.162401.16240/(ptpt+1.162401.16240)-0.120787+2.20275e+05TMath::Exp(-pt2.40913e+01)) ;
	277	}
	278
	279	Double_t rz=(dz-meanZ)/sz ;
	280	Double_t rx=(dx-meanX)/sx ;
	281	return TMath::Sqrt(rxrx+rzrz) ;
	282	}
	283
	284	//____________________________________________________________________________
ba4942ff	285	void AliPHOSClusterSelection::EvalCoreLambdas(AliVCluster * clu, AliVCaloCells * cells,Double_t &m02, Double_t &m20) const
ba4942ff	286	{
17ba2205	287	//calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
17ba2205	288	//Copied from pi0flowtask
ba4942ff	289	Int_t fRunNumber;
17ba2205	290	const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
ba4942ff	291	// Calculates the center of gravity in the local PHOS-module coordinates
17ba2205	292	Float_t wtot = 0;
	293	const Int_t mulDigit=clu->GetNCells() ;
	294	Double_t xc[mulDigit] ;
	295	Double_t zc[mulDigit] ;
	296	Double_t wi[mulDigit] ;
	297	Double_t x = 0 ;
	298	Double_t z = 0 ;
	299	const Double_t logWeight=4.5 ;
	300	for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
	301	Int_t relid[4] ;
	302	Float_t xi=0. ;
	303	Float_t zi=0. ;
	304	Int_t absId = clu->GetCellAbsId(iDigit) ;
	305
ba4942ff	306
	307	AliESDEvent* EventESD = dynamic_cast<AliESDEvent*> (AliPHOSClusterSelection::GetCurrentEvent());//dynamic cast to test for ESD or AOD
	308
	309	AliAODEvent* EventAOD = dynamic_cast<AliAODEvent*> (AliPHOSClusterSelection::GetCurrentEvent());
	310	if(EventESD)
	311	fRunNumber = EventESD->GetRunNumber() ;
	312	if(EventAOD)
	313	fRunNumber = EventESD->GetRunNumber() ;
	314
17ba2205	315	AliOADBContainer geomContainer("phosGeo");//Initialize Geometry
	316	geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
	317	TObjArray matrixes = (TObjArray)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
	318	AliPHOSGeometry * fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
	319	for(Int_t mod=0; mod<5; mod++) {
	320	if(!matrixes->At(mod)) {
ba4942ff	321	AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
17ba2205	322	continue;
	323	}
	324
	325
ba4942ff	326	fPHOSGeo->AbsToRelNumbering(absId, relid) ;
	327	fPHOSGeo->RelPosInModule(relid, xi, zi);
	328	xc[iDigit]=xi ;
	329	zc[iDigit]=zi ;
	330	Double_t ei = elist[iDigit]*cells->GetCellAmplitude(absId) ;
	331	wi[iDigit]=0. ;
	332	if (clu->E()>0 && ei>0) {
	333	wi[iDigit] = TMath::Max( 0., logWeight + TMath::Log( ei / clu->E() ) ) ;
	334	Double_t w=wi[iDigit];
	335	x += xc[iDigit] * w ;
	336	z += zc[iDigit] * w ;
	337	wtot += w ;
	338	}
	339	}
	340	if (wtot>0) {
	341	x /= wtot ;
	342	z /= wtot ;
17ba2205	343	}
17ba2205	344
ba4942ff	345	wtot = 0. ;
	346	Double_t dxx = 0.;
	347	Double_t dzz = 0.;
	348	Double_t dxz = 0.;
	349	Double_t xCut = 0. ;
	350	Double_t zCut = 0. ;
	351	for(Int_t iDigit1=0; iDigit1<mulDigit; iDigit1++) {
	352	Double_t w=wi[iDigit1];
	353	if (w>0.) {
	354	Double_t xi1= xc[iDigit1] ;
	355	Double_t zi1= zc[iDigit1] ;
	356	if((xi1-x)(xi1-x)+(zi1-z)(zi1-z) < 4.5*4.5){
	357	xCut += w * xi1 ;
	358	zCut += w * zi1 ;
	359	dxx += w * xi1 * xi1 ;
	360	dzz += w * zi1 * zi1 ;
	361	dxz += w * xi1 * zi1 ;
17ba2205	362	wtot += w ;
17ba2205	363	}
ba4942ff	364	}
17ba2205	365
ba4942ff	366	}
	367	if (wtot>0) {
	368	xCut/= wtot ;
	369	zCut/= wtot ;
	370	dxx /= wtot ;
	371	dzz /= wtot ;
	372	dxz /= wtot ;
	373	dxx -= xCut * xCut ;
	374	dzz -= zCut * zCut ;
	375	dxz -= xCut * zCut ;
	376
	377	m02 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
	378	m20 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
	379	}
	380	else {
	381	m20=m02=0.;
	382	}
17ba2205	383
ba4942ff	384	}
17ba2205	385	}
	386
	387	//_____________________________________________________________________________
ba4942ff	388	Bool_t AliPHOSClusterSelection::TestLambda(Double_t pt,Double_t l1,Double_t l2) const
ba4942ff	389	{
17ba2205	390	//Evaluates if lambdas correspond to photon cluster
	391	//Tuned using pp data
	392	//copied from Pi0FlowTask
	393	Double_t l2Mean, l1Mean, l2Sigma, l1Sigma, c, R2;
ba4942ff	394	if(fIsCore<0){
17ba2205	395	l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07pt+1.73420e+06pt*pt) ;
	396	l1Mean = 1.12365+0.123770TMath::Exp(-pt0.246551)+5.30000e-03*pt ;
	397	l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11pt+5.01265e+05ptpt)+9.00000e-03pt;
	398	l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00pt+6.78604e-07ptpt)+9.00000e-03pt;
	399	c=-0.35-0.550TMath::Exp(-0.390730pt) ;
	400	R2=0.5(l1-l1Mean)(l1-l1Mean)/l1Sigma/l1Sigma +
	401	0.5(l2-l2Mean)(l2-l2Mean)/l2Sigma/l2Sigma +
	402	0.5c(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
	403	}
	404	else{
	405	//For core radius R=4.5
	406	l1Mean = 1.150200 + 0.097886/(1.+1.486645pt+0.000038pt*pt) ;
	407	l2Mean = 1.574706 + 0.997966exp(-0.895075pt)-0.010666*pt ;
	408	l1Sigma = 0.100255 + 0.337177exp(-0.517684pt)+0.001170*pt ;
	409	l2Sigma = 0.232580 + 0.573401exp(-0.735903pt)-0.002325*pt ;
	410	c = -0.110983 -0.017353/(1.-1.836995pt+0.934517pt*pt) ;
	411	R2=0.5(l1-l1Mean)(l1-l1Mean)/l1Sigma/l1Sigma +
	412	0.5(l2-l2Mean)(l2-l2Mean)/l2Sigma/l2Sigma +
	413	0.5c(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
	414	}
ba4942ff	415	return (R2<fMaxDispR2*fMaxDispR2);//TODO: Check if this should be 2. order
17ba2205	416	}
	417
	418
cc909f95 HQ	419	AliVEvent* AliPHOSClusterSelection::GetCurrentEvent() const
	420	{
	421	// Hackish way of getting the current event.
	422	// Its probably not appropriate to call this function outside of
	423	// AliAnalysisTaskSE::UserExec
	424
	425	AliAnalysisManager* analysisManager = dynamic_cast<AliAnalysisManager*>(AliAnalysisManager::GetAnalysisManager());
	426	AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(analysisManager->GetInputEventHandler());
ba4942ff	427	AliMultiInputEventHandler multiInputHandler = dynamic_cast<AliMultiInputEventHandler >(inputHandler);
cc909f95 HQ	428	if (multiInputHandler)
	429	inputHandler = dynamic_cast<AliInputEventHandler *>(multiInputHandler->GetFirstInputEventHandler());
	430
	431	AliVEvent* inputEvent = dynamic_cast<AliVEvent*>(inputHandler->GetEvent());
	432	if( ! inputEvent )
	433	AliError("Was not able to retrieve event!");
	434
	435	return inputEvent;
	436	}