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