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 ,fGeoUtils(0)
  29 {
  30    fHistogramNameSuffix = TString("PhosMC");
  31 }
  32
  33 AliAnalysisEtMonteCarloPhos::~AliAnalysisEtMonteCarloPhos()
  34 { // dtor
  35   delete fBadMapM2;
  36   delete fBadMapM3;
  37   delete fBadMapM4;
  38   delete fGeoUtils;
  39 }
  40
  41
  42 void AliAnalysisEtMonteCarloPhos::Init()
  43 { // Init
  44   AliAnalysisEtMonteCarlo::Init();
  45   fSelector = new AliAnalysisEtSelectorPhos(fCuts);
  46
  47   fDetectorRadius = fCuts->GetGeometryPhosDetectorRadius();
  48   fSingleCellEnergyCut = fCuts->GetReconstructedPhosSingleCellEnergyCut();
  49
  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 //
  57    fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
  58
  59 }
  60
  61
  62 Bool_t AliAnalysisEtMonteCarloPhos::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
  63 { // too close to bad channel?
  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
 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 }