PWGLF/totEt/AliAnalysisEtSelectorPhos.cxx

   1 #include "AliAnalysisEtSelectorPhos.h"
   2 #include "AliAnalysisEtCuts.h"
   3 #include "AliESDCaloCluster.h"
   4 #include <AliVEvent.h>
   5 #include "TRefArray.h"
   6 #include "AliPHOSGeometry.h"
   7 #include "TH2I.h"
   8 #include "TFile.h"
   9 #include "TMath.h"
  10 #include "AliLog.h"
  11 #include <iostream>
  12 ClassImp(AliAnalysisEtSelectorPhos)
  13 AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
  14 ,fGeoUtils(0)
  15 ,fBadMapM2(0)
  16 ,fBadMapM3(0)
  17 ,fBadMapM4(0)
  18 {
  19
  20 }
  21
  22 AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
  23 {
  24
  25 }
  26
  27 TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
  28 {
  29   return 0;
  30 }
  31
  32 Int_t AliAnalysisEtSelectorPhos::Init(Int_t runNumber)
  33 {
  34   AliAnalysisEtSelector::Init(runNumber);
  35
  36   int res = LoadGeometry();
  37   if(res) return -1;
  38   return LoadBadMaps();
  39 }
  40
  41 Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const AliESDCaloCluster& cluster) const
  42 {
  43   return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
  44 }
  45
  46 Bool_t AliAnalysisEtSelectorPhos::CutDistanceToBadChannel(const AliESDCaloCluster& cluster) const
  47 {
  48     Float_t gPos[3];
  49     cluster.GetPosition(gPos);
  50     Int_t relId[4];
  51     TVector3 glVec(gPos);
  52     fGeoUtils->GlobalPos2RelId(glVec, relId);
  53
  54     TVector3 locVec;
  55     fGeoUtils->Global2Local(locVec, glVec, relId[0]);
  56 //    std::cout << fGeoUtils << std::endl;
  57     //std::cout << relId[0] << " " << cluster.IsPHOS() <<  std::endl;
  58     //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
  59     for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
  60     {
  61         for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
  62         {
  63             if (relId[0] == 3)
  64             {
  65                 if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
  66                 {
  67                     Int_t tmpRel[4];
  68                     tmpRel[0] = 3;
  69                     tmpRel[1] = 0;
  70                     tmpRel[2] = x+1;
  71                     tmpRel[3] = z+1;
  72
  73                     Float_t tmpX;
  74                     Float_t tmpZ;
  75                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
  76
  77                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
  78                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
  79
  80                     if (distance < fCuts->GetPhosBadDistanceCut())
  81                     {
  82 //                    std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
  83                         return kFALSE;
  84                     }
  85                 }
  86             }
  87             if (relId[0] == 2)
  88             {
  89                 if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
  90                 {
  91                     Int_t tmpRel[4];
  92                     tmpRel[0] = 2;
  93                     tmpRel[1] = 0;
  94                     tmpRel[2] = x+1;
  95                     tmpRel[3] = z+1;
  96
  97                     Float_t tmpX;
  98                     Float_t tmpZ;
  99                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 100
 101                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 102
 103 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 104                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 105                     if (distance < fCuts->GetPhosBadDistanceCut())
 106                     {
 107 //                    std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 108                         return kFALSE;
 109                     }
 110                 }
 111             }
 112             if (relId[0] == 1)
 113             {
 114                 if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
 115                 {
 116                     Int_t tmpRel[4];
 117                     tmpRel[0] = 1;
 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 << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 133                         return kFALSE;
 134                     }
 135                 }
 136             }
 137
 138         }
 139     }
 140
 141     return kTRUE;
 142
 143 }
 144
 145 Bool_t AliAnalysisEtSelectorPhos::CutTrackMatching(const AliESDCaloCluster& cluster, Double_t &r) const
 146 {
 147
 148
 149   // cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
 150
 151   Int_t nTracksMatched = cluster.GetNTracksMatched();
 152   if(nTracksMatched == 0) return kFALSE;
 153
 154   Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
 155   if(trackMatchedIndex < 0) return kTRUE;
 156
 157   AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
 158   Double_t dx = cluster.GetTrackDx();
 159   Double_t dz = cluster.GetTrackDz();
 160   Double_t pt = track->Pt();
 161   Int_t charge = track->Charge();
 162
 163   Double_t meanX=0;
 164   Double_t meanZ=0.;
 165   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
 166               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);
 167   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) ;
 168
 169   Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --
 170
 171   if(mf<0.){ //field --
 172     meanZ = -0.468318 ;
 173     if(charge>0)
 174       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)) ;
 175     else
 176       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)) ;
 177   }
 178   else{ //Field ++
 179     meanZ= -0.468318;
 180     if(charge>0)
 181       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)) ;
 182     else
 183       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)) ;
 184   }
 185
 186   Double_t rz=(dz-meanZ)/sz ;
 187   Double_t rx=(dx-meanX)/sx ;
 188   r = TMath::Sqrt(rx*rx+rz*rz);
 189   if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
 190
 191   return kTRUE;
 192
 193 }
 194
 195 int AliAnalysisEtSelectorPhos::LoadGeometry()
 196 {
 197
 198   fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
 199     // ifstream f("badchannels.txt", ios::in);
 200   return 0;
 201 }
 202
 203 int AliAnalysisEtSelectorPhos::LoadBadMaps()
 204 {
 205 TFile *f = TFile::Open("badchannels.root", "READ");
 206
 207     if(!f)
 208     {
 209       std::cout << "Could not open badchannels.root" << std::endl;
 210       return -1;
 211     }
 212
 213     fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
 214     if(fBadMapM2)
 215     {
 216       std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
 217     }
 218     fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
 219     if(fBadMapM2)
 220     {
 221       std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
 222     }
 223
 224     fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
 225     if(fBadMapM4)
 226     {
 227       std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
 228     }
 229
 230     return 0;
 231
 232
 233
 234 }