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

//_________________________________________________________________________
//  Utility Class for transverse energy studies
//  Base class for MC analysis, for PHOS
//  - MC output
//  implementation file 
//
//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
//_________________________________________________________________________
#include "AliAnalysisEtMonteCarloPhos.h"
#include "AliAnalysisEtSelectorPhos.h"
#include "AliAnalysisEtCuts.h"
#include "AliESDtrack.h"
#include <iostream>
#include "AliPHOSGeoUtils.h"
#include "TFile.h"
#include "TH2I.h"
#include <AliPHOSGeometry.h>

using namespace std;

ClassImp(AliAnalysisEtMonteCarloPhos);


AliAnalysisEtMonteCarloPhos::AliAnalysisEtMonteCarloPhos():AliAnalysisEtMonteCarlo()
,fBadMapM2(0)
,fBadMapM3(0)
,fBadMapM4(0)
,fGeoUtils(0)
{
   fHistogramNameSuffix = TString("PhosMC");
}

AliAnalysisEtMonteCarloPhos::~AliAnalysisEtMonteCarloPhos()
{ // dtor
  delete fBadMapM2;
  delete fBadMapM3;
  delete fBadMapM4;
  delete fGeoUtils;
}


void AliAnalysisEtMonteCarloPhos::Init()
{ // Init
  AliAnalysisEtMonteCarlo::Init();
  fSelector = new AliAnalysisEtSelectorPhos(fCuts);
    
  fDetectorRadius = fCuts->GetGeometryPhosDetectorRadius();
  fSingleCellEnergyCut = fCuts->GetReconstructedPhosSingleCellEnergyCut();

 // ifstream f("badchannels.txt", ios::in);
  TFile *f = TFile::Open("badchannels.root", "READ");
  
  fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
   fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
   fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
// 
   fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
  
}


Bool_t AliAnalysisEtMonteCarloPhos::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
{ // too close to bad channel?

    Float_t gPos[3];
    cluster.GetPosition(gPos);
    Int_t relId[4];
    TVector3 glVec(gPos);
    fGeoUtils->GlobalPos2RelId(glVec, relId);

    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 << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
                    
                        return kTRUE;
                    }
                }
            }
            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 << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
                    
                        return kTRUE;
                    }
                }
            }
            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 << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
                    
                        return kTRUE;
                    }
                }
            }

        }
    }

    return kFALSE;
}


void AliAnalysisEtMonteCarloPhos::CreateHistograms()
{ // add some extra histograms & objects to the ones from base class
  if(!fSelector){
    cout<<__FILE__<<" "<<"Creating new fSelector"<<endl;
    fSelector = new AliAnalysisEtSelectorPhos(fCuts);
  }
  AliAnalysisEtMonteCarlo::CreateHistograms();
}
Commit	Line	Data
cf6522d1	1	//_________________________________________________________________________
	2	// Utility Class for transverse energy studies
	3	// Base class for MC analysis, for PHOS
	4	// - MC output
	5	// implementation file
	6	//
	7	//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
	8	//_________________________________________________________________________
2fbf38ac	9	#include "AliAnalysisEtMonteCarloPhos.h"
ef647350	10	#include "AliAnalysisEtSelectorPhos.h"
2fbf38ac	11	#include "AliAnalysisEtCuts.h"
2fbf38ac	12	#include "AliESDtrack.h"
ef647350	13	#include <iostream>
	14	#include "AliPHOSGeoUtils.h"
	15	#include "TFile.h"
	16	#include "TH2I.h"
	17	#include <AliPHOSGeometry.h>
2fbf38ac	18
16abb579	19	using namespace std;
	20
	21	ClassImp(AliAnalysisEtMonteCarloPhos);
	22
	23
ef647350	24	AliAnalysisEtMonteCarloPhos::AliAnalysisEtMonteCarloPhos():AliAnalysisEtMonteCarlo()
	25	,fBadMapM2(0)
	26	,fBadMapM3(0)
	27	,fBadMapM4(0)
f61cec2f	28	,fGeoUtils(0)
2fbf38ac	29	{
	30	fHistogramNameSuffix = TString("PhosMC");
	31	}
	32
cf6522d1	33	AliAnalysisEtMonteCarloPhos::~AliAnalysisEtMonteCarloPhos()
4d376d01	34	{ // dtor
311c6540	35	delete fBadMapM2;
	36	delete fBadMapM3;
	37	delete fBadMapM4;
	38	delete fGeoUtils;
cf6522d1	39	}
cf6522d1	40
2fbf38ac	41
2fbf38ac	42	void AliAnalysisEtMonteCarloPhos::Init()
cf6522d1	43	{ // Init
4998becf	44	AliAnalysisEtMonteCarlo::Init();
ef647350	45	fSelector = new AliAnalysisEtSelectorPhos(fCuts);
2fbf38ac	46
4998becf	47	fDetectorRadius = fCuts->GetGeometryPhosDetectorRadius();
4998becf	48	fSingleCellEnergyCut = fCuts->GetReconstructedPhosSingleCellEnergyCut();
f61cec2f	49
ef647350	50	// ifstream f("badchannels.txt", ios::in);
	51	TFile *f = TFile::Open("badchannels.root", "READ");
	52
	53	fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
	54	fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
	55	fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
	56	//
ef647350	57	fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
ef647350	58
2fbf38ac	59	}
ef647350	60
	61
	62	Bool_t AliAnalysisEtMonteCarloPhos::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
4d376d01	63	{ // too close to bad channel?
ef647350	64
	65	Float_t gPos[3];
	66	cluster.GetPosition(gPos);
	67	Int_t relId[4];
	68	TVector3 glVec(gPos);
	69	fGeoUtils->GlobalPos2RelId(glVec, relId);
	70
	71	TVector3 locVec;
	72	fGeoUtils->Global2Local(locVec, glVec, relId[0]);
	73	// std::cout << fGeoUtils << std::endl;
	74	//std::cout << relId[0] << " " << cluster.IsPHOS() << std::endl;
	75	//std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
	76	for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
	77	{
	78	for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
	79	{
	80	if (relId[0] == 3)
	81	{
	82	if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
	83	{
	84	Int_t tmpRel[4];
	85	tmpRel[0] = 3;
	86	tmpRel[1] = 0;
	87	tmpRel[2] = x+1;
	88	tmpRel[3] = z+1;
	89
	90	Float_t tmpX;
	91	Float_t tmpZ;
	92	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
	93
	94	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
	95	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
	96	if (distance < fCuts->GetPhosBadDistanceCut())
	97	{
	98	// std::cout << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
	99
	100	return kTRUE;
	101	}
	102	}
	103	}
	104	if (relId[0] == 2)
	105	{
	106	if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
	107	{
	108	Int_t tmpRel[4];
	109	tmpRel[0] = 2;
	110	tmpRel[1] = 0;
	111	tmpRel[2] = x+1;
	112	tmpRel[3] = z+1;
	113
	114	Float_t tmpX;
	115	Float_t tmpZ;
	116	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
	117
	118	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
	119
	120	// Float_t distance = TMath::Sqrt((x-locVec[0])(x-locVec[0]) + (z - locVec[2])(z-locVec[2]));
	121	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
	122	if (distance < fCuts->GetPhosBadDistanceCut())
	123	{
	124	// std::cout << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
	125
	126	return kTRUE;
	127	}
128	}
129	}
130	if (relId[0] == 1)
131	{
132	if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
133	{
134	Int_t tmpRel[4];
135	tmpRel[0] = 1;
136	tmpRel[1] = 0;
137	tmpRel[2] = x+1;
138	tmpRel[3] = z+1;
139
140	Float_t tmpX;
141	Float_t tmpZ;
142	fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
143
144	Float_t distance = TMath::Sqrt((tmpX-locVec[0])(tmpX-locVec[0]) + (tmpZ - locVec[2])(tmpZ-locVec[2]));
145
146	// Float_t distance = TMath::Sqrt((x-locVec[0])(x-locVec[0]) + (z - locVec[2])(z-locVec[2]));
147	//Float_t distance = TMath::Sqrt((x-relId[3])(x-relId[3]) + (z - relId[2])(z-relId[2]));
148	if (distance < fCuts->GetPhosBadDistanceCut())
149	{
150	// std::cout << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
151
152	return kTRUE;
153	}
154	}
155	}
156
157	}
158	}
159
160	return kFALSE;
161	}
162
163
2aab9269	164
	165	void AliAnalysisEtMonteCarloPhos::CreateHistograms()
	166	{ // add some extra histograms & objects to the ones from base class
	167	if(!fSelector){
	168	cout<<__FILE__<<" "<<"Creating new fSelector"<<endl;
	169	fSelector = new AliAnalysisEtSelectorPhos(fCuts);
	170	}
	171	AliAnalysisEtMonteCarlo::CreateHistograms();
	172	}