PWGLF/totEt/AliAnalysisEtMonteCarloPhos.cxx

   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 //_________________________________________________________________________
   9 #include "AliAnalysisEtMonteCarloPhos.h"
  10 #include "AliAnalysisEtSelectorPhos.h"
  11 #include "AliAnalysisEtCuts.h"
  12 #include "AliESDtrack.h"
  13 #include <iostream>
  14 #include "AliPHOSGeoUtils.h"
  15 #include "TFile.h"
  16 #include "TH2I.h"
  17 #include <AliPHOSGeometry.h>
  18
  19 using namespace std;
  20
  21 ClassImp(AliAnalysisEtMonteCarloPhos);
  22
  23
  24 AliAnalysisEtMonteCarloPhos::AliAnalysisEtMonteCarloPhos():AliAnalysisEtMonteCarlo()
  25 ,fBadMapM2(0)
  26 ,fBadMapM3(0)
  27 ,fBadMapM4(0)
  28 {
  29    fHistogramNameSuffix = TString("PhosMC");
  30 }
  31
  32 AliAnalysisEtMonteCarloPhos::~AliAnalysisEtMonteCarloPhos()
  33 {
  34 }
  35
  36
  37 void AliAnalysisEtMonteCarloPhos::Init()
  38 { // Init
  39   AliAnalysisEtMonteCarlo::Init();
  40   fSelector = new AliAnalysisEtSelectorPhos(fCuts);
  41
  42     fSelector->Init(137366);
  43    fClusterType = fCuts->GetPhosClusterType();
  44   fDetectorRadius = fCuts->GetGeometryPhosDetectorRadius();
  45   fEtaCutAcc = fCuts->GetGeometryPhosEtaAccCut();
  46   fPhiCutAccMax = fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180.;
  47   fPhiCutAccMin = fCuts->GetGeometryPhosPhiAccMinCut() * TMath::Pi()/180.;
  48   fClusterEnergyCut = fCuts->GetReconstructedPhosClusterEnergyCut();
  49   fSingleCellEnergyCut = fCuts->GetReconstructedPhosSingleCellEnergyCut();
  50
  51   fTrackDistanceCut = fCuts->GetPhosTrackDistanceCut();
  52   fTrackDxCut = fCuts->GetPhosTrackDxCut();
  53   fTrackDzCut = fCuts->GetPhosTrackDzCut();
  54
  55   fDetector = fCuts->GetDetectorPhos();
  56
  57  // ifstream f("badchannels.txt", ios::in);
  58   TFile *f = TFile::Open("badchannels.root", "READ");
  59
  60   fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
  61    fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
  62    fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
  63 //
  64 //  fGeoUtils = new AliPHOSGeoUtils("PHOS", "noCPV");
  65    fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
  66
  67 }
  68
  69
  70 Bool_t AliAnalysisEtMonteCarloPhos::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
  71 {
  72
  73     Float_t gPos[3];
  74     cluster.GetPosition(gPos);
  75     Int_t relId[4];
  76     TVector3 glVec(gPos);
  77     fGeoUtils->GlobalPos2RelId(glVec, relId);
  78
  79     TVector3 locVec;
  80     fGeoUtils->Global2Local(locVec, glVec, relId[0]);
  81 //    std::cout << fGeoUtils << std::endl;
  82     //std::cout << relId[0] << " " << cluster.IsPHOS() <<  std::endl;
  83     //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
  84     for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
  85     {
  86         for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
  87         {
  88             if (relId[0] == 3)
  89             {
  90                 if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
  91                 {
  92                     Int_t tmpRel[4];
  93                     tmpRel[0] = 3;
  94                     tmpRel[1] = 0;
  95                     tmpRel[2] = x+1;
  96                     tmpRel[3] = z+1;
  97
  98                     Float_t tmpX;
  99                     Float_t tmpZ;
 100                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 101
 102                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 103                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 104                     if (distance < fCuts->GetPhosBadDistanceCut())
 105                     {
 106 //                    std::cout << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
 107
 108                         return kTRUE;
 109                     }
 110                 }
 111             }
 112             if (relId[0] == 2)
 113             {
 114                 if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
 115                 {
 116                     Int_t tmpRel[4];
 117                     tmpRel[0] = 2;
 118                     tmpRel[1] = 0;
 119                     tmpRel[2] = x+1;
 120                     tmpRel[3] = z+1;
 121
 122                     Float_t tmpX;
 123                     Float_t tmpZ;
 124                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 125
 126                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 127
 128 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 129                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 130                     if (distance < fCuts->GetPhosBadDistanceCut())
 131                     {
 132 //                      std::cout << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
 133
 134                         return kTRUE;
 135                     }
 136                 }
 137             }
 138             if (relId[0] == 1)
 139             {
 140                 if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
 141                 {
 142                     Int_t tmpRel[4];
 143                     tmpRel[0] = 1;
 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
 154 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 155                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 156                     if (distance < fCuts->GetPhosBadDistanceCut())
 157                     {
 158 //                      std::cout << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
 159
 160                         return kTRUE;
 161                     }
 162                 }
 163             }
 164
 165         }
 166     }
 167
 168     return kFALSE;
 169 }
 170
 171