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