[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>
#include "AliStack.h"

ClassImp(AliAnalysisEtSelectorPhos)

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

AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(): AliAnalysisEtSelector()
,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
{
  
  Float_t pos[3];
  cluster.GetPosition(pos);
  TVector3 cp(pos);
//    std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
  return TMath::Sin(cp.Theta())*cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
}

Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(const TParticle& part) const
{
//    std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
    return TMath::Sin(part.Theta())*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);
  if(track->Pt()<0.5) return kTRUE;//Track matches below about 500 MeV are mostly random.  It takes ~460 MeV to reach the EMCal
  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)
{
  float myphi = part.Phi();
  myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
  return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut() 
    && myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
    && myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
}

Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track)
{
  float myphi = track.Phi();
  myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
  return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
    && myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
    && myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
}


Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliESDCaloCluster& cluster)
{
  Float_t pos[3];
  cluster.GetPosition(pos);
  TVector3 cp(pos);
  float myphi = cp.Phi();
  myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
  return TMath::Abs(cp.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
    && myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
    && myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
}
UInt_t AliAnalysisEtSelectorPhos::GetLabel(const AliESDCaloCluster *cluster, AliStack *stack){
  //Finds primary and estimates if it unique one?
  Int_t n=cluster->GetNLabels() ;
  Int_t iMax=0;

  if (n > 0) {
    Double_t*  Ekin=  new  Double_t[n] ;
    for(Int_t i=0;  i<n;  i++){
      TParticle*  p=  stack->Particle(cluster->GetLabelAt(i)) ;
      Ekin[i]=p->P() ;  // estimate of kinetic energy
      if(p->GetPdgCode()==-2212  ||  p->GetPdgCode()==-2112){
	Ekin[i]+=1.8  ;  //due to annihilation
      }
    }
    Double_t eMax=0.;//eSubMax=0. ;
    for(Int_t i=0;  i<n;  i++){
      if(Ekin[i]>eMax){
	//      eSubMax=eMax;
	eMax=Ekin[i];
	iMax=i;
      }
    }
//   if(eSubMax>0.8*eMax)//not obvious primary
//     sure=kFALSE;
//   else
//     sure=kTRUE;
    delete [] Ekin;

  } // n>0
  return  cluster->GetLabelAt(iMax) ; // DS: should this line be inside n>0 check, and return another value if n<=0 ? 

}
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>
32503dac	21	#include "AliStack.h"
f61cec2f	22
ef647350	23	ClassImp(AliAnalysisEtSelectorPhos)
f61cec2f	24
ef647350	25	AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
	26	,fGeoUtils(0)
	27	,fBadMapM2(0)
	28	,fBadMapM3(0)
	29	,fBadMapM4(0)
f61cec2f	30	,fMatrixInitialized(kFALSE)
c31562f7	31	{
	32
	33	}
	34
	35	AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(): AliAnalysisEtSelector()
	36	,fGeoUtils(0)
	37	,fBadMapM2(0)
	38	,fBadMapM3(0)
	39	,fBadMapM4(0)
	40	,fMatrixInitialized(kFALSE)
ef647350	41	{
	42
	43	}
	44
	45	AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
	46	{
	47
	48	}
	49
	50	TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
4d376d01	51	{ // Get clusters
f61cec2f	52	if(!fClusterArray) fClusterArray = new TRefArray;
	53
	54	if(fClusterArray)
	55	{
	56	fEvent->GetPHOSClusters(fClusterArray);
	57	}
	58	else
	59	{
	60	Printf("Could not initialize cluster array");
	61	}
	62
	63	return fClusterArray;
ef647350	64	}
ef647350	65
f61cec2f	66	Int_t AliAnalysisEtSelectorPhos::Init(const AliESDEvent* event)
4d376d01	67	{ // Init
ef647350	68
f61cec2f	69	AliAnalysisEtSelector::Init(event);
f61cec2f	70	Printf("Initializing selector for run: %d", event->GetRunNumber());
ef647350	71	int res = LoadGeometry();
ef647350	72	if(res) return -1;
f61cec2f	73	if(LoadBadMaps()) return -1;
	74	fInitialized = kTRUE;
	75	if (!fMatrixInitialized)
	76	{
	77	Printf("INITIALIZING MISALIGNMENT MATRICES");
	78	for (Int_t mod=0; mod<5; mod++) {
	79
	80	if (!event->GetPHOSMatrix(mod))
	81	{
	82	Printf("Could not find geo matrix for module %d", mod);
	83	continue;
	84	}
	85	fMatrixInitialized = kTRUE;
	86	fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
	87	Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
	88	}
	89	}
	90	return 0;
ef647350	91	}
ef647350	92
86e7d5db	93	Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(const AliESDCaloCluster& cluster) const
ef647350	94	{
b2c10007	95
cfe23ff0	96	Float_t pos[3];
	97	cluster.GetPosition(pos);
	98	TVector3 cp(pos);
b2c10007	99	// std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
cfe23ff0	100	return TMath::Sin(cp.Theta())*cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
ef647350	101	}
ef647350	102
86e7d5db	103	Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(const TParticle& part) const
f61cec2f	104	{
b2c10007	105	// std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
cfe23ff0	106	return TMath::Sin(part.Theta())*part.Energy() > fCuts->GetReconstructedPhosClusterEnergyCut();
f61cec2f	107	}
	108
	109
86e7d5db	110	Bool_t AliAnalysisEtSelectorPhos::PassDistanceToBadChannelCut(const AliESDCaloCluster& cluster) const
4d376d01	111	{ // cut distance to bad channel
f61cec2f	112	if(!fMatrixInitialized)
	113	{
	114	Printf("Misalignment matrices are not initialized");
	115	return kFALSE;
	116	}
ef647350	117	Float_t gPos[3];
	118	cluster.GetPosition(gPos);
	119	Int_t relId[4];
	120	TVector3 glVec(gPos);
	121	fGeoUtils->GlobalPos2RelId(glVec, relId);
	122
b2c10007	123	//std::cout << "In phos distance to bad channel cut!" << std::endl;
ef647350	124	TVector3 locVec;
	125	fGeoUtils->Global2Local(locVec, glVec, relId[0]);
	126	// std::cout << fGeoUtils << std::endl;
	127	//std::cout << relId[0] << " " << cluster.IsPHOS() << std::endl;
	128	//std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
	129	for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
	130	{
	131	for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
	132	{
	133	if (relId[0] == 3)
	134	{
	135	if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
	136	{
	137	Int_t tmpRel[4];
	138	tmpRel[0] = 3;
	139	tmpRel[1] = 0;
	140	tmpRel[2] = x+1;
	141	tmpRel[3] = z+1;
	142
	143	Float_t tmpX;
	144	Float_t tmpZ;
	145	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
	146
	147	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
	148	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
	149
	150	if (distance < fCuts->GetPhosBadDistanceCut())
	151	{
	152	// std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
	153	return kFALSE;
	154	}
	155	}
	156	}
	157	if (relId[0] == 2)
	158	{
	159	if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
	160	{
	161	Int_t tmpRel[4];
	162	tmpRel[0] = 2;
	163	tmpRel[1] = 0;
	164	tmpRel[2] = x+1;
	165	tmpRel[3] = z+1;
	166
	167	Float_t tmpX;
	168	Float_t tmpZ;
	169	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
	170
	171	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
	172
	173	// Float_t distance = TMath::Sqrt((x-locVec[0])(x-locVec[0]) + (z - locVec[2])(z-locVec[2]));
	174	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
	175	if (distance < fCuts->GetPhosBadDistanceCut())
	176	{
	177	// std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
	178	return kFALSE;
	179	}
	180	}
	181	}
	182	if (relId[0] == 1)
	183	{
	184	if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
	185	{
	186	Int_t tmpRel[4];
	187	tmpRel[0] = 1;
188	tmpRel[1] = 0;
189	tmpRel[2] = x+1;
190	tmpRel[3] = z+1;
191
192	Float_t tmpX;
193	Float_t tmpZ;
194	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
195
196	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
197
198	// Float_t distance = TMath::Sqrt((x-locVec[0])(x-locVec[0]) + (z - locVec[2])(z-locVec[2]));
199	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
200	if (distance < fCuts->GetPhosBadDistanceCut())
201	{
202	// std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
203	return kFALSE;
204	}
205	}
206	}
207
208	}
209	}
210
211	return kTRUE;
212
213	}
214
86e7d5db	215	Bool_t AliAnalysisEtSelectorPhos::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
4d376d01	216	{ // cut track matching
ef647350	217
f61cec2f	218	if(!fMatrixInitialized)
	219	{
	220	Printf("Misalignment matrices are not initialized");
	221	return kFALSE;
	222	}
ef647350	223
ef647350	224	// cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
6a152780	225
ef647350	226	Int_t nTracksMatched = cluster.GetNTracksMatched();
b2c10007	227	if(nTracksMatched == 0) return kTRUE;
ef647350	228
	229	Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
	230	if(trackMatchedIndex < 0) return kTRUE;
	231
	232	AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
ea3ce170	233	if(track->Pt()<0.5) return kTRUE;//Track matches below about 500 MeV are mostly random. It takes ~460 MeV to reach the EMCal
ef647350	234	Double_t dx = cluster.GetTrackDx();
	235	Double_t dz = cluster.GetTrackDz();
	236	Double_t pt = track->Pt();
	237	Int_t charge = track->Charge();
	238
	239	Double_t meanX=0;
	240	Double_t meanZ=0.;
	241	Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
	242	6.58365e-015.91917e-015.91917e-01/((pt-9.61306e-01)(pt-9.61306e-01)+5.91917e-015.91917e-01)+1.59219);
	243	Double_t sz=TMath::Min(2.75,4.90341e+021.91456e-021.91456e-02/(ptpt+1.91456e-021.91456e-02)+1.60) ;
	244
	245	Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --
	246
	247	if(mf<0.){ //field --
	248	meanZ = -0.468318 ;
	249	if(charge>0)
	250	meanX=TMath::Min(7.3, 3.899941.206791.20679/(ptpt+1.206791.20679)+0.249029+2.49088e+07TMath::Exp(-pt3.33650e+01)) ;
	251	else
	252	meanX=-TMath::Min(7.7,3.860400.9124990.912499/(ptpt+0.9124990.912499)+1.23114+4.48277e+05TMath::Exp(-pt2.57070e+01)) ;
	253	}
	254	else{ //Field ++
	255	meanZ= -0.468318;
	256	if(charge>0)
	257	meanX=-TMath::Min(8.0,3.860401.313571.31357/(ptpt+1.313571.31357)+0.880579+7.56199e+06TMath::Exp(-pt3.08451e+01)) ;
	258	else
	259	meanX= TMath::Min(6.85, 3.899941.162401.16240/(ptpt+1.162401.16240)-0.120787+2.20275e+05TMath::Exp(-pt2.40913e+01)) ;
	260	}
	261
	262	Double_t rz=(dz-meanZ)/sz ;
	263	Double_t rx=(dx-meanX)/sx ;
f61cec2f	264	Double_t r = TMath::Sqrt(rxrx+rzrz);
ef647350	265	if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
	266
	267	return kTRUE;
	268
	269	}
	270
	271	int AliAnalysisEtSelectorPhos::LoadGeometry()
4d376d01	272	{ // load geometry
ef647350	273
	274	fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
	275	// ifstream f("badchannels.txt", ios::in);
	276	return 0;
	277	}
	278
	279	int AliAnalysisEtSelectorPhos::LoadBadMaps()
4d376d01	280	{ // load bad maps
ef647350	281	TFile *f = TFile::Open("badchannels.root", "READ");
	282
	283	if(!f)
	284	{
	285	std::cout << "Could not open badchannels.root" << std::endl;
	286	return -1;
	287	}
f61cec2f	288
ef647350	289	fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
f61cec2f	290	if(!fBadMapM2)
ef647350	291	{
	292	std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
	293	}
	294	fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
f61cec2f	295	if(!fBadMapM2)
ef647350	296	{
	297	std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
	298	}
	299
	300	fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
f61cec2f	301	if(!fBadMapM4)
ef647350	302	{
	303	std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
	304	}
ef647350	305
ef647350	306
f61cec2f	307	return 0;
ef647350	308
ef647350	309	}
f61cec2f	310
f61cec2f	311
43dd5a38	312	Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const TParticle& part)
f61cec2f	313	{
43dd5a38	314	float myphi = part.Phi();
43dd5a38	315	myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
f61cec2f	316	return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
43dd5a38	317	&& myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
43dd5a38	318	&& myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
f61cec2f	319	}
f61cec2f	320
43dd5a38	321	Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track)
f61cec2f	322	{
43dd5a38	323	float myphi = track.Phi();
	324	myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
	325	return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
	326	&& myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
	327	&& myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
f61cec2f	328	}
b2c10007	329
31c813d5	330
43dd5a38	331	Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliESDCaloCluster& cluster)
31c813d5	332	{
	333	Float_t pos[3];
	334	cluster.GetPosition(pos);
	335	TVector3 cp(pos);
43dd5a38	336	float myphi = cp.Phi();
	337	myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
	338	return TMath::Abs(cp.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
	339	&& myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
	340	&& myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
31c813d5	341	}
32503dac	342	UInt_t AliAnalysisEtSelectorPhos::GetLabel(const AliESDCaloCluster cluster, AliStack stack){
	343	//Finds primary and estimates if it unique one?
	344	Int_t n=cluster->GetNLabels() ;
32503dac	345	Int_t iMax=0;
beb92504	346
	347	if (n > 0) {
	348	Double_t* Ekin= new Double_t[n] ;
	349	for(Int_t i=0; i<n; i++){
	350	TParticle* p= stack->Particle(cluster->GetLabelAt(i)) ;
	351	Ekin[i]=p->P() ; // estimate of kinetic energy
	352	if(p->GetPdgCode()==-2212 \|\| p->GetPdgCode()==-2112){
	353	Ekin[i]+=1.8 ; //due to annihilation
	354	}
	355	}
	356	Double_t eMax=0.;//eSubMax=0. ;
	357	for(Int_t i=0; i<n; i++){
	358	if(Ekin[i]>eMax){
	359	// eSubMax=eMax;
	360	eMax=Ekin[i];
	361	iMax=i;
	362	}
32503dac	363	}
32503dac	364	// if(eSubMax>0.8*eMax)//not obvious primary
	365	// sure=kFALSE;
	366	// else
	367	// sure=kTRUE;
18875e70	368	delete [] Ekin;
beb92504	369
	370	} // n>0
	371	return cluster->GetLabelAt(iMax) ; // DS: should this line be inside n>0 check, and return another value if n<=0 ?
32503dac	372
32503dac	373	}