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 #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 Bool_t AliAnalysisEtSelectorPhos::PassMinEnergyCut(Double_t e) const
 110 {
 111   return e > fCuts->GetReconstructedEmcalClusterEnergyCut();
 112 }
 113
 114
 115 Bool_t AliAnalysisEtSelectorPhos::PassDistanceToBadChannelCut(const AliESDCaloCluster& cluster) const
 116 { // cut distance to bad channel
 117   if(!fMatrixInitialized)
 118   {
 119     Printf("Misalignment matrices are not initialized");
 120     return kFALSE;
 121   }
 122     Float_t gPos[3];
 123     cluster.GetPosition(gPos);
 124     Int_t relId[4];
 125     TVector3 glVec(gPos);
 126     fGeoUtils->GlobalPos2RelId(glVec, relId);
 127
 128     //std::cout << "In phos distance to bad channel cut!" << std::endl;
 129     TVector3 locVec;
 130     fGeoUtils->Global2Local(locVec, glVec, relId[0]);
 131 //    std::cout << fGeoUtils << std::endl;
 132     //std::cout << relId[0] << " " << cluster.IsPHOS() <<  std::endl;
 133     //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
 134     for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
 135     {
 136         for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
 137         {
 138             if (relId[0] == 3)
 139             {
 140                 if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
 141                 {
 142                     Int_t tmpRel[4];
 143                     tmpRel[0] = 3;
 144                     tmpRel[1] = 0;
 145                     tmpRel[2] = x+1;
 146                     tmpRel[3] = z+1;
 147
 148                     Float_t tmpX;
 149                     Float_t tmpZ;
 150                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 151
 152                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 153                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 154
 155                     if (distance < fCuts->GetPhosBadDistanceCut())
 156                     {
 157 //                    std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 158                         return kFALSE;
 159                     }
 160                 }
 161             }
 162             if (relId[0] == 2)
 163             {
 164                 if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
 165                 {
 166                     Int_t tmpRel[4];
 167                     tmpRel[0] = 2;
 168                     tmpRel[1] = 0;
 169                     tmpRel[2] = x+1;
 170                     tmpRel[3] = z+1;
 171
 172                     Float_t tmpX;
 173                     Float_t tmpZ;
 174                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 175
 176                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 177
 178 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 179                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 180                     if (distance < fCuts->GetPhosBadDistanceCut())
 181                     {
 182 //                    std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 183                         return kFALSE;
 184                     }
 185                 }
 186             }
 187             if (relId[0] == 1)
 188             {
 189                 if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
 190                 {
 191                     Int_t tmpRel[4];
 192                     tmpRel[0] = 1;
 193                     tmpRel[1] = 0;
 194                     tmpRel[2] = x+1;
 195                     tmpRel[3] = z+1;
 196
 197                     Float_t tmpX;
 198                     Float_t tmpZ;
 199                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 200
 201                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 202
 203 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 204                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 205                     if (distance < fCuts->GetPhosBadDistanceCut())
 206                     {
 207 //                      std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 208                         return kFALSE;
 209                     }
 210                 }
 211             }
 212
 213         }
 214     }
 215
 216     return kTRUE;
 217
 218 }
 219
 220 Bool_t AliAnalysisEtSelectorPhos::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
 221 { // cut track matching
 222
 223   if(!fMatrixInitialized)
 224   {
 225     Printf("Misalignment matrices are not initialized");
 226     return kFALSE;
 227   }
 228
 229   // cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
 230
 231   Int_t nTracksMatched = cluster.GetNTracksMatched();
 232   if(nTracksMatched == 0) return kTRUE;
 233
 234   Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
 235   if(trackMatchedIndex < 0) return kTRUE;
 236
 237   AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
 238   if(track->Pt()<0.5) return kTRUE;//Track matches below about 500 MeV are mostly random.  It takes ~460 MeV to reach the EMCal
 239   Double_t dx = cluster.GetTrackDx();
 240   Double_t dz = cluster.GetTrackDz();
 241   Double_t pt = track->Pt();
 242   Int_t charge = track->Charge();
 243
 244   Double_t meanX=0;
 245   Double_t meanZ=0.;
 246   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
 247               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);
 248   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) ;
 249
 250   Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --
 251
 252   if(mf<0.){ //field --
 253     meanZ = -0.468318 ;
 254     if(charge>0)
 255       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)) ;
 256     else
 257       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)) ;
 258   }
 259   else{ //Field ++
 260     meanZ= -0.468318;
 261     if(charge>0)
 262       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)) ;
 263     else
 264       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)) ;
 265   }
 266
 267   Double_t rz=(dz-meanZ)/sz ;
 268   Double_t rx=(dx-meanX)/sx ;
 269   Double_t r = TMath::Sqrt(rx*rx+rz*rz);
 270   if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
 271
 272   return kTRUE;
 273
 274 }
 275
 276 int AliAnalysisEtSelectorPhos::LoadGeometry()
 277 { // load geometry
 278
 279   fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
 280     // ifstream f("badchannels.txt", ios::in);
 281   return 0;
 282 }
 283
 284 int AliAnalysisEtSelectorPhos::LoadBadMaps()
 285 { // load bad maps
 286 TFile *f = TFile::Open("badchannels.root", "READ");
 287
 288     if(!f)
 289     {
 290       std::cout << "Could not open badchannels.root" << std::endl;
 291       return -1;
 292     }
 293
 294     fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
 295     if(!fBadMapM2)
 296     {
 297       std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
 298     }
 299     fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
 300     if(!fBadMapM2)
 301     {
 302       std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
 303     }
 304
 305     fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
 306     if(!fBadMapM4)
 307     {
 308       std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
 309     }
 310
 311
 312     return 0;
 313
 314 }
 315
 316
 317 Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const TParticle& part)
 318 {
 319   float myphi = part.Phi();
 320   myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
 321   return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
 322     && myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
 323     && myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
 324 }
 325
 326 Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track)
 327 {
 328   float myphi = track.Phi();
 329   myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
 330   return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
 331     && myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
 332     && myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
 333 }
 334
 335
 336 Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliESDCaloCluster& cluster)
 337 {
 338   Float_t pos[3];
 339   cluster.GetPosition(pos);
 340   TVector3 cp(pos);
 341   float myphi = cp.Phi();
 342   myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
 343   return TMath::Abs(cp.Eta()) < fCuts->GetGeometryPhosEtaAccCut()
 344     && myphi < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
 345     && myphi > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
 346 }
 347 UInt_t AliAnalysisEtSelectorPhos::GetLabel(const AliESDCaloCluster *cluster, AliStack *stack){
 348   //Finds primary and estimates if it unique one?
 349   Int_t n=cluster->GetNLabels() ;
 350   Int_t iMax=0;
 351
 352   if (n > 0) {
 353     Double_t*  Ekin=  new  Double_t[n] ;
 354     for(Int_t i=0;  i<n;  i++){
 355       TParticle*  p=  stack->Particle(cluster->GetLabelAt(i)) ;
 356       Ekin[i]=p->P() ;  // estimate of kinetic energy
 357       if(p->GetPdgCode()==-2212  ||  p->GetPdgCode()==-2112){
 358         Ekin[i]+=1.8  ;  //due to annihilation
 359       }
 360     }
 361     Double_t eMax=0.;//eSubMax=0. ;
 362     for(Int_t i=0;  i<n;  i++){
 363       if(Ekin[i]>eMax){
 364         //      eSubMax=eMax;
 365         eMax=Ekin[i];
 366         iMax=i;
 367       }
 368     }
 369 //   if(eSubMax>0.8*eMax)//not obvious primary
 370 //     sure=kFALSE;
 371 //   else
 372 //     sure=kTRUE;
 373     delete [] Ekin;
 374
 375   } // n>0
 376   UInt_t correctLabel = cluster->GetLabelAt(iMax);
 377   correctLabel = GetFirstMotherNotFromDetectorCover(correctLabel,*stack);
 378 //   //Now we want to see if this particle is really just something that converted in the cover of the detector and if so, override the label
 379 //   if( stack->IsSecondaryFromMaterial(correctLabel) && correctLabel>0){//if this is flagged as a secondary then we look to see where it really came from
 380 //     TParticle *hitParticle = stack->Particle(correctLabel);
 381 //     if(hitParticle){
 382 //       Bool_t partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) >400);
 383 //       if(partVtxSecondary){//at this point we have something which converted near the detector.  Let's find the mother particle
 384 //      UInt_t mothIdx = stack->Particle(correctLabel)->GetMother(0);
 385 //      if(mothIdx>0){
 386 //        TParticle *mother = stack->Particle(mothIdx);
 387 //        if(mother){
 388 //          partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) >400);
 389 //          if(!partVtxSecondary) return mothIdx;
 390 //          else{
 391 //          }
 392             //      if(AliAnalysisEtSelector::CutGeometricalAcceptance(*mother)){//and the mother is in the acceptance
 393             //        if( !(mother->GetPdgCode()==fgPi0Code)){//some of these are decays that just happen this far out
 394             //          //cout<<"I am declaring that "<<hitParticle->GetName()<<" with a vertex of "<< TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) <<" is actually "<<mother->GetName()<<endl;
 395             //          //cout<<"ID check "<<mothIdx<<" vs "<<mother->GetUniqueID()<<endl;
 396             //          //so now we know that the particle originated near the cover and within the acceptance of the detector
 397             //          return mothIdx;
 398             //        }
 399             //      }
 400 //        }
 401 //      }
 402 //       }
 403
 404 //     }
 405 //   }
 406   return  correctLabel ; // DS: should this line be inside n>0 check, and return another value if n<=0 ?
 407
 408 }