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