[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)
,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::PassMinEnergyCut(const AliESDCaloCluster& cluster) const
{
  
//    std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
  return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
}

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


Bool_t AliAnalysisEtSelectorPhos::PassDistanceToBadChannelCut(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::PassTrackMatchingCut(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;
    
}


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