PWG4/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 "AliAnalysisHadEtCorrections.h"
  28 #include "AliLog.h"
  29 #include "AliCentrality.h"
  30
  31
  32
  33
  34 using namespace std;
  35
  36 ClassImp(AliAnalysisEtReconstructed);
  37
  38
  39 AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
  40         AliAnalysisEt()
  41         ,fCorrections(0)
  42         ,fPidCut(0)
  43         ,fClusterType(0)
  44         ,fHistChargedPionEnergyDeposit(0)
  45         ,fHistProtonEnergyDeposit(0)
  46         ,fHistAntiProtonEnergyDeposit(0)
  47         ,fHistChargedKaonEnergyDeposit(0)
  48         ,fHistMuonEnergyDeposit(0)
  49         ,fHistRemovedEnergy(0)
  50         ,fGeomCorrection(1.0)
  51         ,fEMinCorrection(1.0)
  52 {
  53
  54 }
  55
  56 AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
  57 {
  58     delete fCorrections;
  59 }
  60
  61 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
  62 {
  63     // analyse ESD event
  64     ResetEventValues();
  65
  66     if (!ev) {
  67         AliFatal("ERROR: Event does not exist");
  68         return 0;
  69     }
  70
  71     AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
  72     if (!event) {
  73         AliFatal("ERROR: ESD Event does not exist");
  74         return 0;
  75     }
  76
  77     Int_t cent = -1;
  78     if (fCentrality)
  79     {
  80         cent = fCentrality->GetCentralityClass10("V0M");
  81         fCentClass = fCentrality->GetCentralityClass10("V0M");
  82     }
  83
  84     Double_t protonMass = fgProtonMass;
  85
  86     //for PID
  87     AliESDpid pID;
  88     pID.MakePID(event);
  89     TObjArray* list = fEsdtrackCutsTPC->GetAcceptedTracks(event);
  90
  91     Int_t nGoodTracks = list->GetEntries();
  92     // printf("nGoodTracks %d nCaloClusters %d\n", nGoodTracks, event->GetNumberOfCaloClusters());
  93
  94     for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
  95     {
  96         AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
  97         if (!track)
  98         {
  99             AliError(Form("ERROR: Could not get track %d", iTrack));
 100             continue;
 101         }
 102
 103         fMultiplicity++;
 104
 105
 106         Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
 107         nSigmaPion = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kPion));
 108         nSigmaProton = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kProton));
 109         nSigmaKaon = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kKaon));
 110         nSigmaElectron = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kElectron));
 111         /*
 112         bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
 113         bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
 114         bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0);
 115         bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0);
 116         */
 117
 118         Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0);
 119         Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1);
 120
 121         Float_t massPart = 0;
 122
 123         const Double_t *pidWeights = track->PID();
 124         Int_t maxpid = -1;
 125         Double_t maxpidweight = 0;
 126
 127         if (pidWeights)
 128         {
 129             for (Int_t p =0; p < AliPID::kSPECIES; p++)
 130             {
 131                 if (pidWeights[p] > maxpidweight)
 132                 {
 133                     maxpidweight = pidWeights[p];
 134                     maxpid = p;
 135                 }
 136             }
 137             if (maxpid == AliPID::kProton)
 138             {
 139                 //by definition of ET
 140                 massPart = -protonMass*track->Charge();
 141             }
 142
 143         }
 144
 145         Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
 146
 147         fSparseTracks[0] = maxpid;
 148         fSparseTracks[1] = track->Charge();
 149         fSparseTracks[2] = track->M();
 150         fSparseTracks[3] = et;
 151         fSparseTracks[4] = track->Pt();
 152         fSparseTracks[5] = track->Eta();
 153         fSparseTracks[6] = cent;
 154         fSparseHistTracks->Fill(fSparseTracks);
 155         //printf("Rec track: iTrack %03d eta %4.3f phi %4.3f nITSCl %d nTPCCl %d\n", iTrack, track->Eta(), track->Phi(), nItsClusters, nTPCClusters); // tmp/debug printout
 156
 157         if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() )
 158         {
 159             fTotChargedEt +=  et;
 160             fChargedMultiplicity++;
 161             if (maxpid != -1)
 162             {
 163                 if (maxpid == AliPID::kProton)
 164                 {
 165                     fProtonEt += et;
 166                 }
 167                 if (maxpid == AliPID::kPion)
 168                 {
 169                     fPionEt += et;
 170                 }
 171                 if (maxpid == AliPID::kKaon)
 172                 {
 173                     fChargedKaonEt += et;
 174                 }
 175                 if (maxpid == AliPID::kMuon)
 176                 {
 177                     fMuonEt += et;
 178                 }
 179                 if (maxpid == AliPID::kElectron)
 180                 {
 181                     fElectronEt += et;
 182                 }
 183             }
 184
 185             if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
 186             {
 187                 fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
 188                 if (maxpid != -1)
 189                 {
 190                     if (maxpid == AliPID::kProton)
 191                     {
 192                         fProtonEtAcc += et;
 193                     }
 194                     if (maxpid == AliPID::kPion)
 195                     {
 196                         fPionEtAcc += et;
 197                     }
 198                     if (maxpid == AliPID::kKaon)
 199                     {
 200                         fChargedKaonEtAcc += et;
 201                     }
 202                     if (maxpid == AliPID::kMuon)
 203                     {
 204                         fMuonEtAcc += et;
 205                     }
 206                     if (maxpid == AliPID::kElectron)
 207                     {
 208                         fElectronEtAcc += et;
 209                     }
 210                 }
 211
 212             }
 213         }
 214
 215         if (TrackHitsCalorimeter(track, event->GetMagneticField()))
 216         {
 217             Double_t phi = track->Phi();
 218             Double_t pt = track->Pt();
 219             // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
 220             if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
 221             else fHistPhivsPtNeg->Fill(phi, pt);
 222         }
 223     }
 224
 225     for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
 226     {
 227         AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
 228         if (!cluster)
 229         {
 230             AliError(Form("ERROR: Could not get cluster %d", iCluster));
 231             continue;
 232         }
 233         if (cluster->GetType() != fClusterType) continue;
 234
 235         //if(cluster->GetTracksMatched() > 0)
 236 //      printf("Rec Cluster: iCluster %03d E %4.3f type %.qd NCells %d, nmatched: %d, distance to closest: %f\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells(), cluster->GetNTracksMatched(), cluster->GetEmcCpvDistance()); // tmp/debug printout
 237
 238
 239         if (cluster->E() < fClusterEnergyCut) continue;
 240
 241         Float_t pos[3];
 242
 243         cluster->GetPosition(pos);
 244         TVector3 cp(pos);
 245
 246         Double_t distance = cluster->GetEmcCpvDistance();
 247         Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();
 248         if ( cluster->IsEMCAL() ) {
 249             distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event);
 250         }
 251
 252         fSparseClusters[0] = 0;
 253         fSparseClusters[1] = 0;
 254         fSparseClusters[2] = 0;
 255         fSparseClusters[6] = 0;
 256         fSparseClusters[7] = 0;
 257         fSparseClusters[8] = 0;
 258         fSparseClusters[9] = cent;
 259         fSparseClusters[10] = 0;
 260
 261
 262         if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex > -1)
 263         {
 264             AliVTrack *tmptrack = event->GetTrack(trackMatchedIndex);
 265             if (!tmptrack)
 266             {
 267                 AliError("Error: track does not exist");
 268                 return -1;
 269             }
 270             const Double_t *pidWeights = tmptrack->PID();
 271
 272             Double_t maxpidweight = 0;
 273             Int_t maxpid = 0;
 274             Double_t massPart = 0;
 275             if (pidWeights)
 276             {
 277                 for (Int_t p =0; p < AliPID::kSPECIES; p++)
 278                 {
 279                     if (pidWeights[p] > maxpidweight)
 280                     {
 281                         maxpidweight = pidWeights[p];
 282                         maxpid = p;
 283                     }
 284                 }
 285                 if (maxpid == AliPID::kProton)
 286                 {
 287                     //by definition of ET
 288                     massPart = -protonMass*tmptrack->Charge();
 289                 }
 290             }
 291             fSparseClusters[0] = maxpid;
 292             fSparseClusters[1] = tmptrack->Charge();
 293             fSparseClusters[2] = tmptrack->M();
 294             fSparseClusters[6] = tmptrack->E() * TMath::Sin(tmptrack->Theta()) + massPart;;
 295             fSparseClusters[7] = tmptrack->Pt();
 296             fSparseClusters[8] = tmptrack->Eta();
 297         }
 298
 299         fSparseClusters[10] = distance;
 300
 301         fHistTMDeltaR->Fill(distance);
 302         fHistTMDxDz->Fill(cluster->GetTrackDx(), cluster->GetTrackDz());
 303
 304 //        Float_t clusteret = cluster->E() * TMath::Sin(cp.Theta());
 305
 306         Bool_t matched = false;
 307
 308         if (cluster->IsEMCAL()) matched = distance < fTrackDistanceCut;
 309         else matched = (TMath::Abs(cluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(cluster->GetTrackDz()) < fTrackDzCut);
 310
 311         if (matched)
 312         {
 313             if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
 314             {
 315                 AliVTrack *track = event->GetTrack(trackMatchedIndex);
 316                 if (!track) {
 317                     AliError("Error: track does not exist");
 318                 }
 319                 else {
 320                     const Double_t *pidWeights = track->PID();
 321
 322                     Double_t maxpidweight = 0;
 323                     Int_t maxpid = 0;
 324
 325                     if (pidWeights)
 326                     {
 327                         for (Int_t p =0; p < AliPID::kSPECIES; p++)
 328                         {
 329                             if (pidWeights[p] > maxpidweight)
 330                             {
 331                                 maxpidweight = pidWeights[p];
 332                                 maxpid = p;
 333                             }
 334                         }
 335                         if (fCuts->GetHistMakeTreeDeposit() && fTreeDeposit)
 336                         {
 337                             fEnergyDeposited = cluster->E();
 338                             fEnergyTPC = track->E();
 339                             fCharge = track->Charge();
 340                             fParticlePid = maxpid;
 341                             fPidProb = maxpidweight;
 342                             AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 343                             if (!esdTrack) {
 344                                 AliError("Error: track does not exist");
 345                             }
 346                             else {
 347                                 if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
 348                                 fTreeDeposit->Fill();
 349                             }
 350                         }
 351
 352                         if (maxpidweight > fPidCut)
 353                         {
 354                             Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 355
 356                             Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 357
 358                             Float_t et = cluster->E() * TMath::Sin(theta);
 359                             if (maxpid == AliPID::kProton)
 360                             {
 361
 362                                 if (track->Charge() == 1)
 363                                 {
 364                                     fBaryonEt += et;
 365                                     fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
 366                                 }
 367                                 else if (track->Charge() == -1)
 368                                 {
 369                                     fAntiBaryonEt += et;
 370                                     fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
 371                                 }
 372                             }
 373                             else if (maxpid == AliPID::kPion)
 374                             {
 375                                 fMesonEt += et;
 376                                 fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
 377                             }
 378                             else if (maxpid == AliPID::kKaon)
 379                             {
 380                                 fMesonEt += et;
 381                                 fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
 382                             }
 383                             else if (maxpid == AliPID::kMuon)
 384                             {
 385                                 fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
 386                             }
 387                         }
 388                     }
 389                 }
 390             }
 391             //continue;
 392         } // distance
 393         else
 394         {
 395             fSparseClusters[0] = AliPID::kPhoton;
 396             fSparseClusters[1] = 0;
 397
 398             if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
 399             if (cluster->E() < fClusterEnergyCut) continue;
 400             cluster->GetPosition(pos);
 401
 402             // TODO: replace with TVector3, too lazy now...
 403
 404             float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 405
 406             float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 407             // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
 408             fTotNeutralEt += cluster->E() * TMath::Sin(theta);
 409             fNeutralMultiplicity++;
 410         }
 411
 412         cluster->GetPosition(pos);
 413
 414         // TODO: replace with TVector3
 415
 416         float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 417         float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 418         float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
 419         fSparseClusters[3] = cluster->E() * TMath::Sin(theta);
 420         fSparseClusters[4] = cluster->E();
 421         fSparseClusters[5] = eta;
 422
 423         fSparseHistClusters->Fill(fSparseClusters);
 424
 425         fMultiplicity++;
 426     }
 427     Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity);
 428     fHistRemovedEnergy->Fill(removedEnergy);
 429 //  std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
 430     fTotNeutralEt = fGeomCorrection * fEMinCorrection * fTotNeutralEt - removedEnergy;
 431     fTotNeutralEtAcc = fTotNeutralEt/fGeomCorrection;
 432     fTotEt = fTotChargedEt + fTotNeutralEt;
 433     fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
 434     fSparseEt[0] = fTotEt;
 435     fSparseEt[1] = fTotNeutralEt;
 436     fSparseEt[2] = fTotChargedEtAcc;
 437     fSparseEt[3] = fMultiplicity;
 438     fSparseEt[4] = fNeutralMultiplicity;
 439     fSparseEt[5] = fChargedMultiplicity;
 440     fSparseEt[6] = cent;
 441     // Fill the histograms...
 442     FillHistograms();
 443
 444     return 0;
 445 }
 446
 447 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
 448 { // check vertex
 449
 450     Float_t bxy = 999.;
 451     Float_t bz = 999.;
 452     if (!track) {
 453         AliError("ERROR: no track");
 454         return kFALSE;
 455     }
 456     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 457     if (!esdTrack) {
 458         AliError("ERROR: no track");
 459         return kFALSE;
 460     }
 461     esdTrack->GetImpactParametersTPC(bxy,bz);
 462
 463
 464     bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
 465                   (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
 466                   (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
 467                   (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
 468                   (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
 469
 470     return status;
 471 }
 472
 473 void AliAnalysisEtReconstructed::Init()
 474 { // Init
 475     AliAnalysisEt::Init();
 476     fPidCut = fCuts->GetReconstructedPidCut();
 477     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
 478     if (!fCorrections) {
 479         cout<<"Warning!  You have not set corrections.  Your code will crash.  You have to set the corrections."<<endl;
 480     }
 481 }
 482
 483 bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
 484 { // propagate track to detector radius
 485
 486     if (!track) {
 487         cout<<"Warning: track empty"<<endl;
 488         return kFALSE;
 489     }
 490     AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track);
 491     if (!esdTrack) {
 492         AliError("ERROR: no ESD track");
 493         return kFALSE;
 494     }
 495     // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 496
 497     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 498
 499     // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 500     return prop &&
 501            TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
 502            esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
 503            esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
 504 }
 505
 506 void AliAnalysisEtReconstructed::FillOutputList(TList* list)
 507 { // add some extra histograms to the ones from base class
 508     AliAnalysisEt::FillOutputList(list);
 509
 510     list->Add(fHistChargedPionEnergyDeposit);
 511     list->Add(fHistProtonEnergyDeposit);
 512     list->Add(fHistAntiProtonEnergyDeposit);
 513     list->Add(fHistChargedKaonEnergyDeposit);
 514     list->Add(fHistMuonEnergyDeposit);
 515
 516     list->Add(fHistRemovedEnergy);
 517 }
 518
 519 void AliAnalysisEtReconstructed::CreateHistograms()
 520 { // add some extra histograms to the ones from base class
 521     AliAnalysisEt::CreateHistograms();
 522
 523     Int_t nbinsEt = 1000;
 524     Double_t minEt = 0;
 525     Double_t maxEt = 10;
 526
 527     // possibly change histogram limits
 528     if (fCuts) {
 529         nbinsEt = fCuts->GetHistNbinsParticleEt();
 530         minEt = fCuts->GetHistMinParticleEt();
 531         maxEt = fCuts->GetHistMaxParticleEt();
 532     }
 533
 534     TString histname;
 535     histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
 536     fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 537     fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 538     fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
 539
 540     histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
 541     fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 542     fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 543     fHistProtonEnergyDeposit->SetYTitle("Energy of track");
 544
 545     histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
 546     fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 547     fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 548     fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
 549
 550     histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
 551     fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 552     fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 553     fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
 554
 555     histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
 556     fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 557     fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 558     fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 559
 560     histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
 561     fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
 562     //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 563     //fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 564
 565
 566
 567 }
 568
 569 Double_t
 570 AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3],
 571         Int_t *trkMatchId,
 572         const AliESDEvent *event)
 573 { // calculate distance between cluster and closest track
 574
 575     Double_t trkPos[3] = {0,0,0};
 576
 577     Int_t bestTrkMatchId = -1;
 578     Double_t distance = 9999; // init to a big number
 579
 580     Double_t dist = 0;
 581     Double_t distX = 0, distY = 0, distZ = 0;
 582
 583     for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
 584         AliESDtrack *track = event->GetTrack(iTrack);
 585         if (!track) {
 586             AliError(Form("ERROR: Could not get track %d", iTrack));
 587             continue;
 588         }
 589
 590         // check for approx. eta and phi range before we propagate..
 591         // TBD
 592
 593         AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
 594         if (!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
 595             continue;
 596         }
 597         emctrack->GetXYZ(trkPos);
 598         delete emctrack;
 599
 600         distX = clsPos[0]-trkPos[0];
 601         distY = clsPos[1]-trkPos[1];
 602         distZ = clsPos[2]-trkPos[2];
 603         dist = TMath::Sqrt(distX*distX + distY*distY + distZ*distZ);
 604
 605         if (dist < distance) {
 606             distance = dist;
 607             bestTrkMatchId = iTrack;
 608         }
 609     } // iTrack
 610
 611     // printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
 612     *trkMatchId = bestTrkMatchId;
 613     return distance;
 614 }
 615