PWGLF/totEt/AliAnalysisEtSelectorPhos.cxx

   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 //_________________________________________________________________________
   9 #include "AliAnalysisEtSelectorPhos.h"
  10 #include "AliAnalysisEtCuts.h"
  11 #include "AliESDCaloCluster.h"
  12 #include "AliESDEvent.h"
  13 #include "TRefArray.h"
  14 #include "AliPHOSGeometry.h"
  15 #include "TH2I.h"
  16 #include "TFile.h"
  17 #include "TMath.h"
  18 #include "TParticle.h"
  19 #include "AliLog.h"
  20 #include <iostream>
  21
  22 ClassImp(AliAnalysisEtSelectorPhos)
  23
  24 AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
  25 ,fGeoUtils(0)
  26 ,fBadMapM2(0)
  27 ,fBadMapM3(0)
  28 ,fBadMapM4(0)
  29 ,fMatrixInitialized(kFALSE)
  30 {
  31
  32 }
  33
  34 AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
  35 {
  36
  37 }
  38
  39 TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
  40 { // Get clusters
  41   if(!fClusterArray) fClusterArray = new TRefArray;
  42
  43   if(fClusterArray)
  44   {
  45     fEvent->GetPHOSClusters(fClusterArray);
  46   }
  47   else
  48   {
  49     Printf("Could not initialize cluster array");
  50   }
  51
  52   return fClusterArray;
  53 }
  54
  55 Int_t AliAnalysisEtSelectorPhos::Init(const AliESDEvent* event)
  56 { // Init
  57
  58   AliAnalysisEtSelector::Init(event);
  59   Printf("Initializing selector for run: %d", event->GetRunNumber());
  60   int res = LoadGeometry();
  61   if(res) return -1;
  62   if(LoadBadMaps()) return -1;
  63   fInitialized = kTRUE;
  64   if (!fMatrixInitialized)
  65     {
  66         Printf("INITIALIZING MISALIGNMENT MATRICES");
  67         for (Int_t mod=0; mod<5; mod++) {
  68
  69             if (!event->GetPHOSMatrix(mod))
  70             {
  71               Printf("Could not find geo matrix for module %d", mod);
  72               continue;
  73             }
  74             fMatrixInitialized = kTRUE;
  75             fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
  76             Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
  77         }
  78     }
  79   return 0;
  80 }
  81
  82 Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(const AliESDCaloCluster& cluster) const
  83 {
  84
  85 //    std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
  86   return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
  87 }
  88
  89 Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(const TParticle& part) const
  90 {
  91 //    std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
  92     return part.Energy() > fCuts->GetReconstructedPhosClusterEnergyCut();
  93 }
  94
  95
  96 Bool_t AliAnalysisEtSelectorPhos::PassDistanceToBadChannelCut(const AliESDCaloCluster& cluster) const
  97 { // cut distance to bad channel
  98   if(!fMatrixInitialized)
  99   {
 100     Printf("Misalignment matrices are not initialized");
 101     return kFALSE;
 102   }
 103     Float_t gPos[3];
 104     cluster.GetPosition(gPos);
 105     Int_t relId[4];
 106     TVector3 glVec(gPos);
 107     fGeoUtils->GlobalPos2RelId(glVec, relId);
 108
 109     //std::cout << "In phos distance to bad channel cut!" << std::endl;
 110     TVector3 locVec;
 111     fGeoUtils->Global2Local(locVec, glVec, relId[0]);
 112 //    std::cout << fGeoUtils << std::endl;
 113     //std::cout << relId[0] << " " << cluster.IsPHOS() <<  std::endl;
 114     //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
 115     for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
 116     {
 117         for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
 118         {
 119             if (relId[0] == 3)
 120             {
 121                 if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
 122                 {
 123                     Int_t tmpRel[4];
 124                     tmpRel[0] = 3;
 125                     tmpRel[1] = 0;
 126                     tmpRel[2] = x+1;
 127                     tmpRel[3] = z+1;
 128
 129                     Float_t tmpX;
 130                     Float_t tmpZ;
 131                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 132
 133                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 134                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 135
 136                     if (distance < fCuts->GetPhosBadDistanceCut())
 137                     {
 138 //                    std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 139                         return kFALSE;
 140                     }
 141                 }
 142             }
 143             if (relId[0] == 2)
 144             {
 145                 if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
 146                 {
 147                     Int_t tmpRel[4];
 148                     tmpRel[0] = 2;
 149                     tmpRel[1] = 0;
 150                     tmpRel[2] = x+1;
 151                     tmpRel[3] = z+1;
 152
 153                     Float_t tmpX;
 154                     Float_t tmpZ;
 155                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 156
 157                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 158
 159 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 160                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 161                     if (distance < fCuts->GetPhosBadDistanceCut())
 162                     {
 163 //                    std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 164                         return kFALSE;
 165                     }
 166                 }
 167             }
 168             if (relId[0] == 1)
 169             {
 170                 if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
 171                 {
 172                     Int_t tmpRel[4];
 173                     tmpRel[0] = 1;
 174                     tmpRel[1] = 0;
 175                     tmpRel[2] = x+1;
 176                     tmpRel[3] = z+1;
 177
 178                     Float_t tmpX;
 179                     Float_t tmpZ;
 180                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 181
 182                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 183
 184 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 185                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 186                     if (distance < fCuts->GetPhosBadDistanceCut())
 187                     {
 188 //                      std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 189                         return kFALSE;
 190                     }
 191                 }
 192             }
 193
 194         }
 195     }
 196
 197     return kTRUE;
 198
 199 }
 200
 201 Bool_t AliAnalysisEtSelectorPhos::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
 202 { // cut track matching
 203
 204   if(!fMatrixInitialized)
 205   {
 206     Printf("Misalignment matrices are not initialized");
 207     return kFALSE;
 208   }
 209
 210   // cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
 211
 212   Int_t nTracksMatched = cluster.GetNTracksMatched();
 213   if(nTracksMatched == 0) return kTRUE;
 214
 215   Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
 216   if(trackMatchedIndex < 0) return kTRUE;
 217
 218   AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
 219   Double_t dx = cluster.GetTrackDx();
 220   Double_t dz = cluster.GetTrackDz();
 221   Double_t pt = track->Pt();
 222   Int_t charge = track->Charge();
 223
 224   Double_t meanX=0;
 225   Double_t meanZ=0.;
 226   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
 227               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);
 228   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) ;
 229
 230   Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --
 231
 232   if(mf<0.){ //field --
 233     meanZ = -0.468318 ;
 234     if(charge>0)
 235       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)) ;
 236     else
 237       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)) ;
 238   }
 239   else{ //Field ++
 240     meanZ= -0.468318;
 241     if(charge>0)
 242       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)) ;
 243     else
 244       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)) ;
 245   }
 246
 247   Double_t rz=(dz-meanZ)/sz ;
 248   Double_t rx=(dx-meanX)/sx ;
 249   Double_t r = TMath::Sqrt(rx*rx+rz*rz);
 250   if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
 251
 252   return kTRUE;
 253
 254 }
 255
 256 int AliAnalysisEtSelectorPhos::LoadGeometry()
 257 { // load geometry
 258
 259   fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
 260     // ifstream f("badchannels.txt", ios::in);
 261   return 0;
 262 }
 263
 264 int AliAnalysisEtSelectorPhos::LoadBadMaps()
 265 { // load bad maps
 266 TFile *f = TFile::Open("badchannels.root", "READ");
 267
 268     if(!f)
 269     {
 270       std::cout << "Could not open badchannels.root" << std::endl;
 271       return -1;
 272     }
 273
 274     fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
 275     if(!fBadMapM2)
 276     {
 277       std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
 278     }
 279     fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
 280     if(!fBadMapM2)
 281     {
 282       std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
 283     }
 284
 285     fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
 286     if(!fBadMapM4)
 287     {
 288       std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
 289     }
 290
 291
 292     return 0;
 293
 294 }
 295
 296
 297 Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const TParticle& part) const
 298 {
 299   return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
 300           && part.Phi() < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
 301           && part.Phi() > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
 302 }
 303
 304 Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track) const
 305 {
 306   return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut() &&
 307            track.Phi() > fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180. &&
 308            track.Phi() < fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
 309 }
 310