[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::PassMinEnergyCut(Double_t e) const
{
  return e > fCuts->GetReconstructedEmcalClusterEnergyCut();
}


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
  UInt_t correctLabel = cluster->GetLabelAt(iMax);
  correctLabel = GetFirstMotherNotFromDetectorCover(correctLabel,*stack);
//   //Now we want to see if this particle is really just something that converted in the cover of the detector and if so, override the label
//   if( stack->IsSecondaryFromMaterial(correctLabel) && correctLabel>0){//if this is flagged as a secondary then we look to see where it really came from
//     TParticle *hitParticle = stack->Particle(correctLabel);
//     if(hitParticle){
//       Bool_t partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) >400);
//       if(partVtxSecondary){//at this point we have something which converted near the detector.  Let's find the mother particle
//  	UInt_t mothIdx = stack->Particle(correctLabel)->GetMother(0);
//  	if(mothIdx>0){
//  	  TParticle *mother = stack->Particle(mothIdx);
//  	  if(mother){
// 	    partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) >400);
// 	    if(!partVtxSecondary) return mothIdx;
// 	    else{
// 	    }
	    // 	    if(AliAnalysisEtSelector::CutGeometricalAcceptance(*mother)){//and the mother is in the acceptance
	    // 	      if( !(mother->GetPdgCode()==fgPi0Code)){//some of these are decays that just happen this far out
	    // 		//cout<<"I am declaring that "<<hitParticle->GetName()<<" with a vertex of "<< TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) <<" is actually "<<mother->GetName()<<endl;
	    // 		//cout<<"ID check "<<mothIdx<<" vs "<<mother->GetUniqueID()<<endl;
	    // 		//so now we know that the particle originated near the cover and within the acceptance of the detector
	    // 		return mothIdx;
	    // 	      }
	    // 	    }
//  	  }
//  	}
//       }
      
//     }
//   }
  return  correctLabel ; // 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	}
f61cec2f	108
02d47689	109	Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(Double_t e) const
	110	{
	111	return e > fCuts->GetReconstructedEmcalClusterEnergyCut();
	112	}
	113
f61cec2f	114
86e7d5db	115	Bool_t AliAnalysisEtSelectorPhos::PassDistanceToBadChannelCut(const AliESDCaloCluster& cluster) const
4d376d01	116	{ // cut distance to bad channel
f61cec2f	117	if(!fMatrixInitialized)
	118	{
	119	Printf("Misalignment matrices are not initialized");
	120	return kFALSE;
	121	}
ef647350	122	Float_t gPos[3];
	123	cluster.GetPosition(gPos);
	124	Int_t relId[4];
	125	TVector3 glVec(gPos);
	126	fGeoUtils->GlobalPos2RelId(glVec, relId);
	127
b2c10007	128	//std::cout << "In phos distance to bad channel cut!" << std::endl;
ef647350	129	TVector3 locVec;
	130	fGeoUtils->Global2Local(locVec, glVec, relId[0]);
	131	// std::cout << fGeoUtils << std::endl;
	132	//std::cout << relId[0] << " " << cluster.IsPHOS() << std::endl;
	133	//std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
	134	for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
	135	{
	136	for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
	137	{
	138	if (relId[0] == 3)
	139	{
	140	if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
	141	{
	142	Int_t tmpRel[4];
	143	tmpRel[0] = 3;
	144	tmpRel[1] = 0;
	145	tmpRel[2] = x+1;
	146	tmpRel[3] = z+1;
	147
	148	Float_t tmpX;
	149	Float_t tmpZ;
	150	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
	151
	152	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
	153	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
	154
	155	if (distance < fCuts->GetPhosBadDistanceCut())
	156	{
	157	// std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
	158	return kFALSE;
	159	}
	160	}
	161	}
	162	if (relId[0] == 2)
	163	{
	164	if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
	165	{
	166	Int_t tmpRel[4];
	167	tmpRel[0] = 2;
	168	tmpRel[1] = 0;
	169	tmpRel[2] = x+1;
	170	tmpRel[3] = z+1;
	171
	172	Float_t tmpX;
	173	Float_t tmpZ;
	174	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
	175
	176	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
	177
	178	// Float_t distance = TMath::Sqrt((x-locVec[0])(x-locVec[0]) + (z - locVec[2])(z-locVec[2]));
	179	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
	180	if (distance < fCuts->GetPhosBadDistanceCut())
	181	{
	182	// std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
	183	return kFALSE;
	184	}
	185	}
	186	}
	187	if (relId[0] == 1)
	188	{
	189	if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
	190	{
	191	Int_t tmpRel[4];
	192	tmpRel[0] = 1;
193	tmpRel[1] = 0;
194	tmpRel[2] = x+1;
195	tmpRel[3] = z+1;
196
197	Float_t tmpX;
198	Float_t tmpZ;
199	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
200
201	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
202
203	// Float_t distance = TMath::Sqrt((x-locVec[0])(x-locVec[0]) + (z - locVec[2])(z-locVec[2]));
204	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
205	if (distance < fCuts->GetPhosBadDistanceCut())
206	{
207	// std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
208	return kFALSE;
209	}
210	}
211	}
212
213	}
214	}
215
216	return kTRUE;
217
218	}
219
86e7d5db	220	Bool_t AliAnalysisEtSelectorPhos::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
4d376d01	221	{ // cut track matching
ef647350	222
f61cec2f	223	if(!fMatrixInitialized)
	224	{
	225	Printf("Misalignment matrices are not initialized");
	226	return kFALSE;
	227	}
ef647350	228
ef647350	229	// cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
6a152780	230
ef647350	231	Int_t nTracksMatched = cluster.GetNTracksMatched();
b2c10007	232	if(nTracksMatched == 0) return kTRUE;
ef647350	233
	234	Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
	235	if(trackMatchedIndex < 0) return kTRUE;
	236
	237	AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
ea3ce170	238	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	239	Double_t dx = cluster.GetTrackDx();
	240	Double_t dz = cluster.GetTrackDz();
	241	Double_t pt = track->Pt();
	242	Int_t charge = track->Charge();
	243
	244	Double_t meanX=0;
	245	Double_t meanZ=0.;
	246	Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
	247	6.58365e-015.91917e-015.91917e-01/((pt-9.61306e-01)(pt-9.61306e-01)+5.91917e-015.91917e-01)+1.59219);
	248	Double_t sz=TMath::Min(2.75,4.90341e+021.91456e-021.91456e-02/(ptpt+1.91456e-021.91456e-02)+1.60) ;
	249
	250	Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --
	251
	252	if(mf<0.){ //field --
	253	meanZ = -0.468318 ;
	254	if(charge>0)
	255	meanX=TMath::Min(7.3, 3.899941.206791.20679/(ptpt+1.206791.20679)+0.249029+2.49088e+07TMath::Exp(-pt3.33650e+01)) ;
	256	else
	257	meanX=-TMath::Min(7.7,3.860400.9124990.912499/(ptpt+0.9124990.912499)+1.23114+4.48277e+05TMath::Exp(-pt2.57070e+01)) ;
	258	}
	259	else{ //Field ++
	260	meanZ= -0.468318;
	261	if(charge>0)
	262	meanX=-TMath::Min(8.0,3.860401.313571.31357/(ptpt+1.313571.31357)+0.880579+7.56199e+06TMath::Exp(-pt3.08451e+01)) ;
	263	else
	264	meanX= TMath::Min(6.85, 3.899941.162401.16240/(ptpt+1.162401.16240)-0.120787+2.20275e+05TMath::Exp(-pt2.40913e+01)) ;
	265	}
	266
	267	Double_t rz=(dz-meanZ)/sz ;
	268	Double_t rx=(dx-meanX)/sx ;
f61cec2f	269	Double_t r = TMath::Sqrt(rxrx+rzrz);
ef647350	270	if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
	271
	272	return kTRUE;
	273
	274	}
	275
	276	int AliAnalysisEtSelectorPhos::LoadGeometry()
4d376d01	277	{ // load geometry
ef647350	278
	279	fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
	280	// ifstream f("badchannels.txt", ios::in);
	281	return 0;
	282	}
	283
	284	int AliAnalysisEtSelectorPhos::LoadBadMaps()
4d376d01	285	{ // load bad maps
ef647350	286	TFile *f = TFile::Open("badchannels.root", "READ");
	287
	288	if(!f)
	289	{
	290	std::cout << "Could not open badchannels.root" << std::endl;
	291	return -1;
	292	}
f61cec2f	293
ef647350	294	fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
f61cec2f	295	if(!fBadMapM2)
ef647350	296	{
	297	std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
	298	}
	299	fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
f61cec2f	300	if(!fBadMapM2)
ef647350	301	{
	302	std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
	303	}
	304
	305	fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
f61cec2f	306	if(!fBadMapM4)
ef647350	307	{
	308	std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
	309	}
ef647350	310
ef647350	311
f61cec2f	312	return 0;
ef647350	313
ef647350	314	}
f61cec2f	315
f61cec2f	316
43dd5a38	317	Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const TParticle& part)
f61cec2f	318	{
43dd5a38	319	float myphi = part.Phi();
43dd5a38	320	myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
f61cec2f	321	return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
43dd5a38	322	&& myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
43dd5a38	323	&& myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
f61cec2f	324	}
f61cec2f	325
43dd5a38	326	Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track)
f61cec2f	327	{
43dd5a38	328	float myphi = track.Phi();
	329	myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
	330	return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
	331	&& myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
	332	&& myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
f61cec2f	333	}
b2c10007	334
31c813d5	335
43dd5a38	336	Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliESDCaloCluster& cluster)
31c813d5	337	{
	338	Float_t pos[3];
	339	cluster.GetPosition(pos);
	340	TVector3 cp(pos);
43dd5a38	341	float myphi = cp.Phi();
	342	myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
	343	return TMath::Abs(cp.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
	344	&& myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
	345	&& myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
31c813d5	346	}
32503dac	347	UInt_t AliAnalysisEtSelectorPhos::GetLabel(const AliESDCaloCluster cluster, AliStack stack){
	348	//Finds primary and estimates if it unique one?
	349	Int_t n=cluster->GetNLabels() ;
32503dac	350	Int_t iMax=0;
beb92504	351
	352	if (n > 0) {
	353	Double_t* Ekin= new Double_t[n] ;
	354	for(Int_t i=0; i<n; i++){
	355	TParticle* p= stack->Particle(cluster->GetLabelAt(i)) ;
	356	Ekin[i]=p->P() ; // estimate of kinetic energy
	357	if(p->GetPdgCode()==-2212 \|\| p->GetPdgCode()==-2112){
	358	Ekin[i]+=1.8 ; //due to annihilation
	359	}
	360	}
	361	Double_t eMax=0.;//eSubMax=0. ;
	362	for(Int_t i=0; i<n; i++){
	363	if(Ekin[i]>eMax){
	364	// eSubMax=eMax;
	365	eMax=Ekin[i];
	366	iMax=i;
	367	}
32503dac	368	}
32503dac	369	// if(eSubMax>0.8*eMax)//not obvious primary
	370	// sure=kFALSE;
	371	// else
	372	// sure=kTRUE;
18875e70	373	delete [] Ekin;
beb92504	374
beb92504	375	} // n>0
6844c491	376	UInt_t correctLabel = cluster->GetLabelAt(iMax);
	377	correctLabel = GetFirstMotherNotFromDetectorCover(correctLabel,*stack);
	378	// //Now we want to see if this particle is really just something that converted in the cover of the detector and if so, override the label
	379	// if( stack->IsSecondaryFromMaterial(correctLabel) && correctLabel>0){//if this is flagged as a secondary then we look to see where it really came from
	380	// TParticle *hitParticle = stack->Particle(correctLabel);
	381	// if(hitParticle){
	382	// Bool_t partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()hitParticle->Vx() + hitParticle->Vy()hitParticle->Vy()) >400);
	383	// if(partVtxSecondary){//at this point we have something which converted near the detector. Let's find the mother particle
	384	// UInt_t mothIdx = stack->Particle(correctLabel)->GetMother(0);
	385	// if(mothIdx>0){
	386	// TParticle *mother = stack->Particle(mothIdx);
	387	// if(mother){
	388	// partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()hitParticle->Vx() + hitParticle->Vy()hitParticle->Vy()) >400);
	389	// if(!partVtxSecondary) return mothIdx;
	390	// else{
	391	// }
	392	// if(AliAnalysisEtSelector::CutGeometricalAcceptance(*mother)){//and the mother is in the acceptance
	393	// if( !(mother->GetPdgCode()==fgPi0Code)){//some of these are decays that just happen this far out
	394	// //cout<<"I am declaring that "<<hitParticle->GetName()<<" with a vertex of "<< TMath::Sqrt(hitParticle->Vx()hitParticle->Vx() + hitParticle->Vy()hitParticle->Vy()) <<" is actually "<<mother->GetName()<<endl;
	395	// //cout<<"ID check "<<mothIdx<<" vs "<<mother->GetUniqueID()<<endl;
	396	// //so now we know that the particle originated near the cover and within the acceptance of the detector
	397	// return mothIdx;
	398	// }
	399	// }
	400	// }
	401	// }
	402	// }
	403
	404	// }
	405	// }
	406	return correctLabel ; // DS: should this line be inside n>0 check, and return another value if n<=0 ?
32503dac	407
32503dac	408	}