PWGLF/totEt/AliAnalysisEtReconstructed.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies
   3 //  Base class for ESD analysis
   4 //  - reconstruction output
   5 //  implementation file
   6 //
   7 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
   8 //_________________________________________________________________________
   9
  10 #include "AliAnalysisEtReconstructed.h"
  11 #include "AliAnalysisEtCuts.h"
  12 #include "AliESDtrack.h"
  13 #include "AliEMCALTrack.h"
  14 #include "AliESDCaloCluster.h"
  15 #include "TVector3.h"
  16 #include "TGeoGlobalMagField.h"
  17 #include "AliMagF.h"
  18 #include "AliVEvent.h"
  19 #include "AliESDEvent.h"
  20 #include "AliESDtrackCuts.h"
  21 #include "AliVParticle.h"
  22 #include "TDatabasePDG.h"
  23 #include "TList.h"
  24 #include "AliESDpid.h"
  25 #include <iostream>
  26 #include "TH2F.h"
  27 #include "TH2I.h"
  28 #include "TH1I.h"
  29 #include "TFile.h"
  30 #include "AliAnalysisHadEtCorrections.h"
  31 #include "AliAnalysisEtSelector.h"
  32 #include "AliLog.h"
  33 #include "AliCentrality.h"
  34 #include "AliPHOSGeoUtils.h"
  35 #include "AliPHOSGeometry.h"
  36
  37
  38 using namespace std;
  39
  40 ClassImp(AliAnalysisEtReconstructed);
  41
  42
  43 AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
  44     AliAnalysisEt()
  45     ,fCorrections(0)
  46     ,fPidCut(0)
  47     ,fHistChargedPionEnergyDeposit(0)
  48     ,fHistProtonEnergyDeposit(0)
  49     ,fHistAntiProtonEnergyDeposit(0)
  50     ,fHistChargedKaonEnergyDeposit(0)
  51     ,fHistMuonEnergyDeposit(0)
  52     ,fHistRemovedEnergy(0)
  53     ,fGeomCorrection(1.0)
  54     ,fEMinCorrection(1.0/0.89)
  55     ,fRecEffCorrection(1.0)
  56     ,fGeoUtils(0)
  57     ,fBadMapM2(0)
  58     ,fBadMapM3(0)
  59     ,fBadMapM4(0)
  60     ,fClusterPosition(0)
  61     ,fHistChargedEnergyRemoved(0)
  62     ,fHistNeutralEnergyRemoved(0)
  63     ,fHistGammaEnergyAdded(0)
  64 {
  65
  66 }
  67
  68 AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
  69 {   //destructor
  70     delete fCorrections;
  71     delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
  72     delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
  73     delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
  74     delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
  75     delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
  76
  77     delete fHistRemovedEnergy; // removed energy
  78
  79 }
  80
  81 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
  82 {
  83     // analyse ESD event
  84     ResetEventValues();
  85
  86     if (!ev) {
  87         AliFatal("ERROR: Event does not exist");
  88         return 0;
  89     }
  90
  91     AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
  92     if (!event) {
  93         AliFatal("ERROR: ESD Event does not exist");
  94         return 0;
  95     }
  96     fSelector->SetEvent(event);
  97     if (!fMatrixInitialized)
  98     {
  99         for (Int_t mod=0; mod<5; mod++) {
 100             if (!event->GetPHOSMatrix(mod)) continue;
 101             fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
 102 //          std::cout << event->GetPHOSMatrix(mod) << std::endl;
 103             Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
 104         }
 105         fMatrixInitialized = kTRUE;
 106     }
 107
 108     Int_t cent = -1;
 109     if (fCentrality)
 110     {
 111         cent = fCentrality->GetCentralityClass10("V0M");
 112         fCentClass = fCentrality->GetCentralityClass10("V0M");
 113     }
 114
 115     //Double_t protonMass = fgProtonMass;
 116
 117     for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
 118     {
 119         AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
 120         if (!cluster)
 121         {
 122             AliError(Form("ERROR: Could not get cluster %d", iCluster));
 123             continue;
 124         }
 125         int x = 0;
 126         fCutFlow->Fill(x++);
 127         if(cluster->IsEMCAL()) continue;
 128         fCutFlow->Fill(x++);
 129         if(!fSelector->CutMinEnergy(*cluster)) continue;
 130         fCutFlow->Fill(x++);
 131         if (!fSelector->CutDistanceToBadChannel(*cluster)) continue;
 132         fCutFlow->Fill(x++);
 133
 134         Float_t pos[3];
 135
 136         cluster->GetPosition(pos);
 137         TVector3 cp(pos);
 138
 139         Double_t distance = cluster->GetEmcCpvDistance();
 140         Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();
 141         if ( cluster->IsEMCAL() ) {
 142             distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event);
 143         }
 144
 145         Bool_t matched = false;
 146
 147         Double_t r = 0;
 148
 149         matched = !fSelector->CutTrackMatching(*cluster, r);
 150
 151         if (matched)
 152         {
 153
 154             if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
 155             {
 156                 AliVTrack *track = event->GetTrack(trackMatchedIndex);
 157                 if (!track) {
 158                     AliError("Error: track does not exist");
 159                 }
 160                 else {
 161                     const Double_t *pidWeights = track->PID();
 162
 163                     Double_t maxpidweight = 0;
 164                     Int_t maxpid = 0;
 165
 166                     if (pidWeights)
 167                     {
 168                         for (Int_t p =0; p < AliPID::kSPECIES; p++)
 169                         {
 170                             if (pidWeights[p] > maxpidweight)
 171                             {
 172                                 maxpidweight = pidWeights[p];
 173                                 maxpid = p;
 174                             }
 175                         }
 176                         if (fCuts->GetHistMakeTreeDeposit() && fTreeDeposit)
 177                         {
 178                             fEnergyDeposited = cluster->E();
 179                             fEnergyTPC = track->E();
 180                             fCharge = track->Charge();
 181                             fParticlePid = maxpid;
 182                             fPidProb = maxpidweight;
 183                             AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 184                             if (!esdTrack) {
 185                                 AliError("Error: track does not exist");
 186                             }
 187                             else {
 188                                 if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
 189                                 fTreeDeposit->Fill();
 190                             }
 191                         }
 192
 193                         if (maxpidweight > fPidCut)
 194                         {
 195                             Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 196
 197                             Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 198
 199                             Float_t et = cluster->E() * TMath::Sin(theta);
 200                             if (maxpid == AliPID::kProton)
 201                             {
 202
 203                                 if (track->Charge() == 1)
 204                                 {
 205                                     fBaryonEt += et;
 206                                     fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
 207                                 }
 208                                 else if (track->Charge() == -1)
 209                                 {
 210                                     fAntiBaryonEt += et;
 211                                     fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
 212                                 }
 213                             }
 214                             else if (maxpid == AliPID::kPion)
 215                             {
 216                                 fMesonEt += et;
 217                                 fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
 218                             }
 219                             else if (maxpid == AliPID::kKaon)
 220                             {
 221                                 fMesonEt += et;
 222                                 fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
 223                             }
 224                             else if (maxpid == AliPID::kMuon)
 225                             {
 226                                 fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
 227                             }
 228                         }
 229                     }
 230                 }
 231             }
 232             //continue;
 233         } // distance
 234         else
 235         {
 236             fCutFlow->Fill(x++);
 237             //std::cout << x++ << std::endl;
 238             fSparseClusters[0] = AliPID::kPhoton;
 239             fSparseClusters[1] = 0;
 240
 241             //if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
 242             //if (cluster->E() < fClusterEnergyCut) continue;
 243             cluster->GetPosition(pos);
 244
 245             TVector3 p2(pos);
 246
 247             fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());
 248
 249             fTotNeutralEt += CalculateTransverseEnergy(cluster);
 250             fNeutralMultiplicity++;
 251         }
 252         fMultiplicity++;
 253     }
 254
 255     fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
 256     fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
 257     fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
 258     fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
 259
 260     fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
 261     fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
 262
 263     Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity);
 264     fHistRemovedEnergy->Fill(removedEnergy);
 265
 266     fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
 267     fTotNeutralEtAcc = fTotNeutralEt;
 268     fTotEt = fTotChargedEt + fTotNeutralEt;
 269     fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
 270     if(fMakeSparse) {
 271         fSparseEt[0] = fTotEt;
 272         fSparseEt[1] = fTotNeutralEt;
 273         fSparseEt[2] = fTotChargedEtAcc;
 274         fSparseEt[3] = fMultiplicity;
 275         fSparseEt[4] = fNeutralMultiplicity;
 276         fSparseEt[5] = fChargedMultiplicity;
 277         fSparseEt[6] = cent;
 278     }
 279     // Fill the histograms...
 280     FillHistograms();
 281
 282     return 0;
 283 }
 284
 285 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
 286 {   // check vertex
 287
 288     Float_t bxy = 999.;
 289     Float_t bz = 999.;
 290     if (!track) {
 291         AliError("ERROR: no track");
 292         return kFALSE;
 293     }
 294     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 295     if (!esdTrack) {
 296         AliError("ERROR: no track");
 297         return kFALSE;
 298     }
 299     esdTrack->GetImpactParametersTPC(bxy,bz);
 300
 301
 302     bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
 303                   (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
 304                   (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
 305                   (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
 306                   (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
 307
 308     return status;
 309 }
 310
 311 void AliAnalysisEtReconstructed::Init()
 312 {   // Init
 313     AliAnalysisEt::Init();
 314     fPidCut = fCuts->GetReconstructedPidCut();
 315     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
 316     if (!fCorrections) {
 317         cout<<"Warning!  You have not set corrections.  Your code will crash.  You have to set the corrections."<<endl;
 318     }
 319     //fGeoUtils = new AliPHOSGeoUtils("PHOS", "noCPV");
 320     fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
 321     // ifstream f("badchannels.txt", ios::in);
 322     TFile *f = TFile::Open("badchannels.root", "READ");
 323
 324     fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
 325     fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
 326     fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
 327 }
 328
 329 bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
 330 {   // propagate track to detector radius
 331
 332     if (!track) {
 333         cout<<"Warning: track empty"<<endl;
 334         return kFALSE;
 335     }
 336     AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track);
 337     if (!esdTrack) {
 338         AliError("ERROR: no ESD track");
 339         return kFALSE;
 340     }
 341     // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 342
 343     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 344
 345     // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 346     return prop &&
 347            TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
 348            esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
 349            esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
 350 }
 351
 352 void AliAnalysisEtReconstructed::FillOutputList(TList* list)
 353 {   // add some extra histograms to the ones from base class
 354     AliAnalysisEt::FillOutputList(list);
 355
 356     list->Add(fHistChargedPionEnergyDeposit);
 357     list->Add(fHistProtonEnergyDeposit);
 358     list->Add(fHistAntiProtonEnergyDeposit);
 359     list->Add(fHistChargedKaonEnergyDeposit);
 360     list->Add(fHistMuonEnergyDeposit);
 361
 362     list->Add(fHistRemovedEnergy);
 363     list->Add(fClusterPosition);
 364
 365     list->Add(fHistChargedEnergyRemoved);
 366     list->Add(fHistNeutralEnergyRemoved);
 367     list->Add(fHistGammaEnergyAdded);
 368 }
 369
 370 void AliAnalysisEtReconstructed::CreateHistograms()
 371 {   // add some extra histograms to the ones from base class
 372     AliAnalysisEt::CreateHistograms();
 373
 374     Int_t nbinsEt = 1000;
 375     Double_t minEt = 0;
 376     Double_t maxEt = 10;
 377
 378     // possibly change histogram limits
 379     if (fCuts) {
 380         nbinsEt = fCuts->GetHistNbinsParticleEt();
 381         minEt = fCuts->GetHistMinParticleEt();
 382         maxEt = fCuts->GetHistMaxParticleEt();
 383     }
 384
 385     TString histname;
 386     histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
 387     fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 388     fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 389     fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
 390
 391     histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
 392     fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 393     fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 394     fHistProtonEnergyDeposit->SetYTitle("Energy of track");
 395
 396     histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
 397     fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 398     fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 399     fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
 400
 401     histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
 402     fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 403     fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 404     fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
 405
 406     histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
 407     fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 408     fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 409     fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 410
 411     histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
 412     fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
 413     //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 414     //fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 415
 416     histname = "fClusterPosition" + fHistogramNameSuffix;
 417     fClusterPosition = new TH2D(histname.Data(), "Position of accepted neutral clusters",1000, -2.0, -.5, 1000, -.13 , 0.13);
 418     fClusterPosition->SetXTitle("Energy deposited in calorimeter");
 419     fClusterPosition->SetYTitle("Energy of track");
 420
 421
 422     histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
 423     fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 424
 425     histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
 426     fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 427
 428     histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
 429     fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 430
 431
 432 }
 433
 434 Double_t
 435 AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3],
 436         Int_t *trkMatchId,
 437         const AliESDEvent *event)
 438 {   // calculate distance between cluster and closest track
 439
 440     Double_t trkPos[3] = {0,0,0};
 441
 442     Int_t bestTrkMatchId = -1;
 443     Double_t distance = 9999; // init to a big number
 444
 445     Double_t dist = 0;
 446     Double_t distX = 0, distY = 0, distZ = 0;
 447
 448     for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
 449         AliESDtrack *track = event->GetTrack(iTrack);
 450         if (!track) {
 451             AliError(Form("ERROR: Could not get track %d", iTrack));
 452             continue;
 453         }
 454
 455         // check for approx. eta and phi range before we propagate..
 456         // TBD
 457
 458         AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
 459         if (!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
 460             continue;
 461         }
 462         emctrack->GetXYZ(trkPos);
 463         delete emctrack;
 464
 465         distX = clsPos[0]-trkPos[0];
 466         distY = clsPos[1]-trkPos[1];
 467         distZ = clsPos[2]-trkPos[2];
 468         dist = TMath::Sqrt(distX*distX + distY*distY + distZ*distZ);
 469
 470         if (dist < distance) {
 471             distance = dist;
 472             bestTrkMatchId = iTrack;
 473         }
 474     } // iTrack
 475
 476     // printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
 477     *trkMatchId = bestTrkMatchId;
 478     return distance;
 479 }
 480
 481
 482 Bool_t AliAnalysisEtReconstructed::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
 483 {
 484
 485     Float_t gPos[3];
 486     cluster.GetPosition(gPos);
 487     Int_t relId[4];
 488     TVector3 glVec(gPos);
 489     fGeoUtils->GlobalPos2RelId(glVec, relId);
 490
 491     TVector3 locVec;
 492     fGeoUtils->Global2Local(locVec, glVec, relId[0]);
 493 //    std::cout << fGeoUtils << std::endl;
 494     //std::cout << relId[0] << " " << cluster.IsPHOS() <<  std::endl;
 495     //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
 496     for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
 497     {
 498         for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
 499         {
 500             if (relId[0] == 3)
 501             {
 502                 if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
 503                 {
 504                     Int_t tmpRel[4];
 505                     tmpRel[0] = 3;
 506                     tmpRel[1] = 0;
 507                     tmpRel[2] = x+1;
 508                     tmpRel[3] = z+1;
 509
 510                     Float_t tmpX;
 511                     Float_t tmpZ;
 512                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 513
 514                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 515                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 516
 517                     if (distance < fCuts->GetPhosBadDistanceCut())
 518                     {
 519 //                    std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 520                         return kTRUE;
 521                     }
 522                 }
 523             }
 524             if (relId[0] == 2)
 525             {
 526                 if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
 527                 {
 528                     Int_t tmpRel[4];
 529                     tmpRel[0] = 2;
 530                     tmpRel[1] = 0;
 531                     tmpRel[2] = x+1;
 532                     tmpRel[3] = z+1;
 533
 534                     Float_t tmpX;
 535                     Float_t tmpZ;
 536                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 537
 538                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 539
 540 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 541                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 542                     if (distance < fCuts->GetPhosBadDistanceCut())
 543                     {
 544 //                    std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 545                         return kTRUE;
 546                     }
 547                 }
 548             }
 549             if (relId[0] == 1)
 550             {
 551                 if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
 552                 {
 553                     Int_t tmpRel[4];
 554                     tmpRel[0] = 1;
 555                     tmpRel[1] = 0;
 556                     tmpRel[2] = x+1;
 557                     tmpRel[3] = z+1;
 558
 559                     Float_t tmpX;
 560                     Float_t tmpZ;
 561                     fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
 562
 563                     Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
 564
 565 //                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 566                     //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
 567                     if (distance < fCuts->GetPhosBadDistanceCut())
 568                     {
 569 //                      std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
 570                         return kTRUE;
 571                     }
 572                 }
 573             }
 574
 575         }
 576     }
 577
 578     return kFALSE;
 579
 580
 581 }
 582
 583
 584
 585
 586 /*
 587 Bool_t AliAnalysisEtReconstructed::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
 588 {
 589
 590     Float_t gPos[3];
 591
 592     cluster.GetPosition(gPos);
 593     Int_t relId[4];
 594     TVector3 glVec(gPos);
 595     fGeoUtils->GlobalPos2RelId(glVec, relId);
 596     TVector3 locVec;
 597     fGeoUtils->Global2Local(locVec, glVec, relId[0]);
 598
 599     std::vector<Int_t>::const_iterator badIt;
 600
 601     for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
 602     {
 603         for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
 604         {
 605             if (relId[0] == 3)
 606             {
 607                 if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
 608                 {
 609
 610                     Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 611                     if (distance < fCuts->GetPhosBadDistanceCut())
 612                     {
 613                         return kTRUE;
 614                     }
 615                 }
 616             }
 617             if (relId[0] == 2)
 618             {
 619                 if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
 620                 {
 621
 622                     Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 623                     if (distance < fCuts->GetPhosBadDistanceCut())
 624                     {
 625                         return kTRUE;
 626                     }
 627                 }
 628             }
 629             if (relId[0] == 1)
 630             {
 631                 if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
 632                 {
 633
 634                     Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
 635                     if (distance < fCuts->GetPhosBadDistanceCut())
 636                     {
 637                         return kTRUE;
 638                     }
 639                 }
 640             }
 641         }
 642     }
 643
 644     return kFALSE;
 645 }
 646 */