[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtSelectorPhos.cxx

//_________________________________________________________________________
//  Utility Class for transverse energy studies
//  Selector Base class for PHOS
//  -  
// implementation file
//
//*-- Authors: Oystein Djuvsland (Bergen)
//_________________________________________________________________________
#include "AliAnalysisEtSelectorPhos.h"
#include "AliAnalysisEtCuts.h"
#include "AliESDCaloCluster.h"
#include "AliESDEvent.h"
#include "TRefArray.h"
#include "AliPHOSGeometry.h"
#include "TH2I.h"
#include "TFile.h"
#include "TMath.h"
#include "TParticle.h"
#include "AliLog.h"
#include <iostream>

ClassImp(AliAnalysisEtSelectorPhos)

AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
,fGeoUtils(0)
,fBadMapM2(0)
,fBadMapM3(0)
,fBadMapM4(0)
,fInitialized(kFALSE)
,fMatrixInitialized(kFALSE)
{
  
}

AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
{

}

TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
{ // Get clusters
  if(!fClusterArray) fClusterArray = new TRefArray;
  
  if(fClusterArray)
  {
    fEvent->GetPHOSClusters(fClusterArray);
  }
  else
  {
    Printf("Could not initialize cluster array");
  }
  
  return fClusterArray;
}

Int_t AliAnalysisEtSelectorPhos::Init(const AliESDEvent* event)
{ // Init
  
  AliAnalysisEtSelector::Init(event);
  Printf("Initializing selector for run: %d", event->GetRunNumber());
  int res = LoadGeometry();
  if(res) return -1;
  if(LoadBadMaps()) return -1;
  fInitialized = kTRUE;
  if (!fMatrixInitialized)
    {
	Printf("INITIALIZING MISALIGNMENT MATRICES");
        for (Int_t mod=0; mod<5; mod++) {
	    
            if (!event->GetPHOSMatrix(mod))
	    {
	      Printf("Could not find geo matrix for module %d", mod);
	      continue;
	    }
	    fMatrixInitialized = kTRUE;
            fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
            Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
        }
    }
  return 0;
}

Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const AliESDCaloCluster& cluster) const
{
  
//    std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
  return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
}

Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const TParticle& part) const
{
//    std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
    return part.Energy() > fCuts->GetReconstructedPhosClusterEnergyCut();
}


Bool_t AliAnalysisEtSelectorPhos::CutDistanceToBadChannel(const AliESDCaloCluster& cluster) const
{ // cut distance to bad channel
  if(!fMatrixInitialized)
  {
    Printf("Misalignment matrices are not initialized");
    return kFALSE;
  }
    Float_t gPos[3];
    cluster.GetPosition(gPos);
    Int_t relId[4];
    TVector3 glVec(gPos);
    fGeoUtils->GlobalPos2RelId(glVec, relId);

    //std::cout << "In phos distance to bad channel cut!" << std::endl;
    TVector3 locVec;
    fGeoUtils->Global2Local(locVec, glVec, relId[0]);
//    std::cout << fGeoUtils << std::endl;
    //std::cout << relId[0] << " " << cluster.IsPHOS() <<  std::endl;
    //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
    for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
    {
        for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
        {
            if (relId[0] == 3)
            {
                if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
                {
		    Int_t tmpRel[4];
		    tmpRel[0] = 3;
		    tmpRel[1] = 0;
		    tmpRel[2] = x+1;
		    tmpRel[3] = z+1;
		    
		    Float_t tmpX;
		    Float_t tmpZ;
		    fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);

                    Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
                    //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
		    
                    if (distance < fCuts->GetPhosBadDistanceCut())
                    {
//		      std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
                        return kFALSE;
                    }
                }
            }
            if (relId[0] == 2)
            {
                if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
                {
		    Int_t tmpRel[4];
		    tmpRel[0] = 2;
		    tmpRel[1] = 0;
		    tmpRel[2] = x+1;
		    tmpRel[3] = z+1;
		    
		    Float_t tmpX;
		    Float_t tmpZ;
		    fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);

                    Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));

//                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
		    //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
                    if (distance < fCuts->GetPhosBadDistanceCut())
                    {
//		      std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
                        return kFALSE;
                    }
                }
            }
            if (relId[0] == 1)
            {
                if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
                {
		    Int_t tmpRel[4];
		    tmpRel[0] = 1;
		    tmpRel[1] = 0;
		    tmpRel[2] = x+1;
		    tmpRel[3] = z+1;
		    
		    Float_t tmpX;
		    Float_t tmpZ;
		    fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);

                    Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));

//                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
		    //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
                    if (distance < fCuts->GetPhosBadDistanceCut())
                    {
//			std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
                        return kFALSE;
                    }
                }
            }

        }
    }

    return kTRUE;

}

Bool_t AliAnalysisEtSelectorPhos::CutTrackMatching(const AliESDCaloCluster& cluster) const
{ // cut track matching

  if(!fMatrixInitialized)
  {
    Printf("Misalignment matrices are not initialized");
    return kFALSE;
  }
  
  // cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
  
  Int_t nTracksMatched = cluster.GetNTracksMatched();
  if(nTracksMatched == 0) return kTRUE;
  
  Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
  if(trackMatchedIndex < 0) return kTRUE;
  
  AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
  Double_t dx = cluster.GetTrackDx();
  Double_t dz = cluster.GetTrackDz();
  Double_t pt = track->Pt();
  Int_t charge = track->Charge();
  
  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) ;
  
  Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --

  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 ;
  Double_t r = TMath::Sqrt(rx*rx+rz*rz);
  if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
  
  return kTRUE;
 
}

int AliAnalysisEtSelectorPhos::LoadGeometry()
{ // load geometry

  fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
    // ifstream f("badchannels.txt", ios::in);
  return 0;
}

int AliAnalysisEtSelectorPhos::LoadBadMaps()
{ // load bad maps
TFile *f = TFile::Open("badchannels.root", "READ");

    if(!f)
    {
      std::cout << "Could not open badchannels.root" << std::endl;
      return -1;
    }

    fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
    if(!fBadMapM2) 
    {
      std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
    }
    fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
    if(!fBadMapM2) 
    {
      std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
    }
    
    fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
    if(!fBadMapM4) 
    {
      std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
    }
    
    
    return 0;
    
}

void AliAnalysisEtSelectorPhos::SetEvent(const AliESDEvent* event)
{ // set event
    //AliAnalysisEtSelector::SetEvent(event);
    fEvent = event;
    if(!fInitialized) Init(event);
}

Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const TParticle& part) const
{
  return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut() 
	  && part.Phi() < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
	  && part.Phi() > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
}

Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track) const
{
  return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut() &&
           track.Phi() > fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180. &&
           track.Phi() < fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
}
Commit	Line	Data
4d376d01	1	//_________________________________________________________________________
	2	// Utility Class for transverse energy studies
	3	// Selector Base class for PHOS
	4	// -
	5	// implementation file
	6	//
	7	//*-- Authors: Oystein Djuvsland (Bergen)
	8	//_________________________________________________________________________
ef647350	9	#include "AliAnalysisEtSelectorPhos.h"
	10	#include "AliAnalysisEtCuts.h"
	11	#include "AliESDCaloCluster.h"
f61cec2f	12	#include "AliESDEvent.h"
ef647350	13	#include "TRefArray.h"
	14	#include "AliPHOSGeometry.h"
	15	#include "TH2I.h"
	16	#include "TFile.h"
	17	#include "TMath.h"
f61cec2f	18	#include "TParticle.h"
ef647350	19	#include "AliLog.h"
ef647350	20	#include <iostream>
f61cec2f	21
ef647350	22	ClassImp(AliAnalysisEtSelectorPhos)
f61cec2f	23
ef647350	24	AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
	25	,fGeoUtils(0)
	26	,fBadMapM2(0)
	27	,fBadMapM3(0)
	28	,fBadMapM4(0)
f61cec2f	29	,fInitialized(kFALSE)
f61cec2f	30	,fMatrixInitialized(kFALSE)
ef647350	31	{
	32
	33	}
	34
	35	AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
	36	{
	37
	38	}
	39
	40	TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
4d376d01	41	{ // Get clusters
f61cec2f	42	if(!fClusterArray) fClusterArray = new TRefArray;
	43
	44	if(fClusterArray)
	45	{
	46	fEvent->GetPHOSClusters(fClusterArray);
	47	}
	48	else
	49	{
	50	Printf("Could not initialize cluster array");
	51	}
	52
	53	return fClusterArray;
ef647350	54	}
ef647350	55
f61cec2f	56	Int_t AliAnalysisEtSelectorPhos::Init(const AliESDEvent* event)
4d376d01	57	{ // Init
ef647350	58
f61cec2f	59	AliAnalysisEtSelector::Init(event);
f61cec2f	60	Printf("Initializing selector for run: %d", event->GetRunNumber());
ef647350	61	int res = LoadGeometry();
ef647350	62	if(res) return -1;
f61cec2f	63	if(LoadBadMaps()) return -1;
	64	fInitialized = kTRUE;
	65	if (!fMatrixInitialized)
	66	{
	67	Printf("INITIALIZING MISALIGNMENT MATRICES");
	68	for (Int_t mod=0; mod<5; mod++) {
	69
	70	if (!event->GetPHOSMatrix(mod))
	71	{
	72	Printf("Could not find geo matrix for module %d", mod);
	73	continue;
	74	}
	75	fMatrixInitialized = kTRUE;
	76	fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
	77	Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
	78	}
	79	}
	80	return 0;
ef647350	81	}
	82
	83	Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const AliESDCaloCluster& cluster) const
	84	{
b2c10007	85
b2c10007	86	// std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
ef647350	87	return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
	88	}
	89
f61cec2f	90	Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const TParticle& part) const
f61cec2f	91	{
b2c10007	92	// std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
f61cec2f	93	return part.Energy() > fCuts->GetReconstructedPhosClusterEnergyCut();
	94	}
	95
	96
ef647350	97	Bool_t AliAnalysisEtSelectorPhos::CutDistanceToBadChannel(const AliESDCaloCluster& cluster) const
4d376d01	98	{ // cut distance to bad channel
f61cec2f	99	if(!fMatrixInitialized)
	100	{
	101	Printf("Misalignment matrices are not initialized");
	102	return kFALSE;
	103	}
ef647350	104	Float_t gPos[3];
	105	cluster.GetPosition(gPos);
	106	Int_t relId[4];
	107	TVector3 glVec(gPos);
	108	fGeoUtils->GlobalPos2RelId(glVec, relId);
	109
b2c10007	110	//std::cout << "In phos distance to bad channel cut!" << std::endl;
ef647350	111	TVector3 locVec;
	112	fGeoUtils->Global2Local(locVec, glVec, relId[0]);
	113	// std::cout << fGeoUtils << std::endl;
	114	//std::cout << relId[0] << " " << cluster.IsPHOS() << std::endl;
	115	//std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
	116	for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
	117	{
	118	for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
	119	{
	120	if (relId[0] == 3)
	121	{
	122	if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
	123	{
	124	Int_t tmpRel[4];
	125	tmpRel[0] = 3;
	126	tmpRel[1] = 0;
	127	tmpRel[2] = x+1;
	128	tmpRel[3] = z+1;
	129
	130	Float_t tmpX;
	131	Float_t tmpZ;
	132	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
	133
	134	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
	135	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
	136
	137	if (distance < fCuts->GetPhosBadDistanceCut())
	138	{
	139	// std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
	140	return kFALSE;
	141	}
	142	}
	143	}
	144	if (relId[0] == 2)
	145	{
	146	if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
	147	{
	148	Int_t tmpRel[4];
	149	tmpRel[0] = 2;
	150	tmpRel[1] = 0;
	151	tmpRel[2] = x+1;
	152	tmpRel[3] = z+1;
	153
	154	Float_t tmpX;
	155	Float_t tmpZ;
	156	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
	157
	158	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
	159
	160	// Float_t distance = TMath::Sqrt((x-locVec[0])(x-locVec[0]) + (z - locVec[2])(z-locVec[2]));
	161	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
	162	if (distance < fCuts->GetPhosBadDistanceCut())
	163	{
	164	// std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
	165	return kFALSE;
	166	}
	167	}
	168	}
	169	if (relId[0] == 1)
	170	{
	171	if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
	172	{
	173	Int_t tmpRel[4];
	174	tmpRel[0] = 1;
175	tmpRel[1] = 0;
176	tmpRel[2] = x+1;
177	tmpRel[3] = z+1;
178
179	Float_t tmpX;
180	Float_t tmpZ;
181	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
182
183	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
184
185	// Float_t distance = TMath::Sqrt((x-locVec[0])(x-locVec[0]) + (z - locVec[2])(z-locVec[2]));
186	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
187	if (distance < fCuts->GetPhosBadDistanceCut())
188	{
189	// std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
190	return kFALSE;
191	}
192	}
193	}
194
195	}
196	}
197
198	return kTRUE;
199
200	}
201
f61cec2f	202	Bool_t AliAnalysisEtSelectorPhos::CutTrackMatching(const AliESDCaloCluster& cluster) const
4d376d01	203	{ // cut track matching
ef647350	204
f61cec2f	205	if(!fMatrixInitialized)
	206	{
	207	Printf("Misalignment matrices are not initialized");
	208	return kFALSE;
	209	}
ef647350	210
	211	// cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
	212
	213	Int_t nTracksMatched = cluster.GetNTracksMatched();
b2c10007	214	if(nTracksMatched == 0) return kTRUE;
ef647350	215
	216	Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
	217	if(trackMatchedIndex < 0) return kTRUE;
	218
	219	AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
	220	Double_t dx = cluster.GetTrackDx();
	221	Double_t dz = cluster.GetTrackDz();
	222	Double_t pt = track->Pt();
	223	Int_t charge = track->Charge();
	224
	225	Double_t meanX=0;
	226	Double_t meanZ=0.;
	227	Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
	228	6.58365e-015.91917e-015.91917e-01/((pt-9.61306e-01)(pt-9.61306e-01)+5.91917e-015.91917e-01)+1.59219);
	229	Double_t sz=TMath::Min(2.75,4.90341e+021.91456e-021.91456e-02/(ptpt+1.91456e-021.91456e-02)+1.60) ;
	230
	231	Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --
	232
	233	if(mf<0.){ //field --
	234	meanZ = -0.468318 ;
	235	if(charge>0)
	236	meanX=TMath::Min(7.3, 3.899941.206791.20679/(ptpt+1.206791.20679)+0.249029+2.49088e+07TMath::Exp(-pt3.33650e+01)) ;
	237	else
	238	meanX=-TMath::Min(7.7,3.860400.9124990.912499/(ptpt+0.9124990.912499)+1.23114+4.48277e+05TMath::Exp(-pt2.57070e+01)) ;
	239	}
	240	else{ //Field ++
	241	meanZ= -0.468318;
	242	if(charge>0)
	243	meanX=-TMath::Min(8.0,3.860401.313571.31357/(ptpt+1.313571.31357)+0.880579+7.56199e+06TMath::Exp(-pt3.08451e+01)) ;
	244	else
	245	meanX= TMath::Min(6.85, 3.899941.162401.16240/(ptpt+1.162401.16240)-0.120787+2.20275e+05TMath::Exp(-pt2.40913e+01)) ;
	246	}
	247
	248	Double_t rz=(dz-meanZ)/sz ;
	249	Double_t rx=(dx-meanX)/sx ;
f61cec2f	250	Double_t r = TMath::Sqrt(rxrx+rzrz);
ef647350	251	if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
	252
	253	return kTRUE;
	254
	255	}
	256
	257	int AliAnalysisEtSelectorPhos::LoadGeometry()
4d376d01	258	{ // load geometry
ef647350	259
	260	fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
	261	// ifstream f("badchannels.txt", ios::in);
	262	return 0;
	263	}
	264
	265	int AliAnalysisEtSelectorPhos::LoadBadMaps()
4d376d01	266	{ // load bad maps
ef647350	267	TFile *f = TFile::Open("badchannels.root", "READ");
	268
	269	if(!f)
	270	{
	271	std::cout << "Could not open badchannels.root" << std::endl;
	272	return -1;
	273	}
f61cec2f	274
ef647350	275	fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
f61cec2f	276	if(!fBadMapM2)
ef647350	277	{
	278	std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
	279	}
	280	fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
f61cec2f	281	if(!fBadMapM2)
ef647350	282	{
	283	std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
	284	}
	285
	286	fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
f61cec2f	287	if(!fBadMapM4)
ef647350	288	{
	289	std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
	290	}
ef647350	291
ef647350	292
f61cec2f	293	return 0;
ef647350	294
ef647350	295	}
f61cec2f	296
f61cec2f	297	void AliAnalysisEtSelectorPhos::SetEvent(const AliESDEvent* event)
4d376d01	298	{ // set event
f61cec2f	299	//AliAnalysisEtSelector::SetEvent(event);
	300	fEvent = event;
	301	if(!fInitialized) Init(event);
	302	}
	303
	304	Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const TParticle& part) const
	305	{
	306	return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
	307	&& part.Phi() < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
	308	&& part.Phi() > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
	309	}
	310
	311	Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track) const
	312	{
	313	return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut() &&
	314	track.Phi() > fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180. &&
	315	track.Phi() < fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
	316	}
b2c10007	317