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.687)
  55         ,fRecEffCorrection(1.0)
  56         ,fClusterPosition(0)
  57         ,fClusterEnergy(0)
  58         ,fClusterEt(0)
  59         ,fHistChargedEnergyRemoved(0)
  60         ,fHistNeutralEnergyRemoved(0)
  61         ,fHistGammaEnergyAdded(0)
  62 {
  63
  64 }
  65
  66 AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
  67 {//destructor
  68     delete fCorrections;
  69     delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
  70     delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
  71     delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
  72     delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
  73     delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
  74
  75     delete fHistRemovedEnergy; // removed energy
  76     delete fClusterPosition;
  77     delete fClusterEnergy;
  78     delete fClusterEt;
  79     delete fHistChargedEnergyRemoved;
  80     delete fHistNeutralEnergyRemoved;
  81     delete fHistGammaEnergyAdded;
  82
  83 }
  84
  85 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
  86 {
  87
  88     //AliAnalysisEt::AnalyseEvent(ev);
  89     // analyse ESD event
  90     ResetEventValues();
  91     if (!ev) {
  92         AliFatal("ERROR: Event does not exist");
  93         return 0;
  94     }
  95
  96     AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
  97     if (!event) {
  98         AliFatal("ERROR: ESD Event does not exist");
  99         return 0;
 100     }
 101     if(!fSelector){
 102         AliFatal("ERROR: fSelector does not exist");
 103         return 0;
 104     }
 105     fSelector->SetEvent(event);
 106
 107     Int_t cent = -1;
 108     if (fCentrality && cent)
 109     {
 110         cent = fCentrality->GetCentralityClass10("V0M");
 111         fCentClass = fCentrality->GetCentralityClass10("V0M");
 112     }
 113
 114     for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
 115     {
 116         AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
 117         if (!cluster)
 118         {
 119             AliError(Form("ERROR: Could not get cluster %d", iCluster));
 120             continue;
 121         }
 122         int x = 0;
 123         fCutFlow->Fill(x++);
 124         if(!fSelector->IsDetectorCluster(*cluster)) continue;
 125         fCutFlow->Fill(x++);
 126         if(!fSelector->PassMinEnergyCut(*cluster)) continue;
 127         fCutFlow->Fill(x++);
 128         if (!fSelector->PassDistanceToBadChannelCut(*cluster)) continue;
 129         fCutFlow->Fill(x++);
 130
 131         Float_t pos[3];
 132
 133         cluster->GetPosition(pos);
 134         TVector3 cp(pos);
 135
 136         Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();
 137
 138         Bool_t matched = false;
 139
 140
 141         matched = !fSelector->PassTrackMatchingCut(*cluster);
 142
 143         if (matched)
 144         {
 145
 146             if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
 147             {
 148                 AliVTrack *track = event->GetTrack(trackMatchedIndex);
 149                 if (!track) {
 150                     AliError("Error: track does not exist");
 151                 }
 152                 else {
 153                     const Double_t *pidWeights = track->PID();
 154
 155                     Double_t maxpidweight = 0;
 156                     Int_t maxpid = 0;
 157
 158                     if (pidWeights)
 159                     {
 160                         for (Int_t p =0; p < AliPID::kSPECIES; p++)
 161                         {
 162                             if (pidWeights[p] > maxpidweight)
 163                             {
 164                                 maxpidweight = pidWeights[p];
 165                                 maxpid = p;
 166                             }
 167                         }
 168                         if (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
 169                         {
 170                             fEnergyDeposited = cluster->E();
 171                             fMomentumTPC = track->P();
 172                             fCharge = track->Charge();
 173                             fParticlePid = maxpid;
 174                             fPidProb = maxpidweight;
 175                             AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 176                             if (!esdTrack) {
 177                                 AliError("Error: track does not exist");
 178                             }
 179                             else {
 180                                 if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
 181                                 fDepositTree->Fill();
 182                             }
 183                         }
 184
 185                         if (maxpidweight > fPidCut)
 186                         {
 187                             //Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 188
 189                             //Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 190
 191                             //Float_t et = cluster->E() * TMath::Sin(theta);
 192                             if (maxpid == AliPID::kProton)
 193                             {
 194
 195                                 if (track->Charge() == 1)
 196                                 {
 197                                     fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
 198                                 }
 199                                 else if (track->Charge() == -1)
 200                                 {
 201                                     fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
 202                                 }
 203                             }
 204                             else if (maxpid == AliPID::kPion)
 205                             {
 206                                 fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
 207                             }
 208                             else if (maxpid == AliPID::kKaon)
 209                             {
 210                                 fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
 211                             }
 212                             else if (maxpid == AliPID::kMuon)
 213                             {
 214                                 fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
 215                             }
 216                         }
 217                     }
 218                 }
 219             }
 220             //continue;
 221         } // distance
 222         else
 223         {
 224           fCutFlow->Fill(x++);
 225           //std::cout << x++ << std::endl;
 226
 227             //if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
 228             //if (cluster->E() < fClusterEnergyCut) continue;
 229             cluster->GetPosition(pos);
 230
 231             TVector3 p2(pos);
 232
 233             fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());
 234             fClusterEnergy->Fill(cluster->E());
 235             fClusterEt->Fill(TMath::Sin(p2.Theta())*cluster->E());
 236
 237             fTotNeutralEt += CorrectForReconstructionEfficiency(*cluster);
 238             fNeutralMultiplicity++;
 239         }
 240         fMultiplicity++;
 241     }
 242
 243     fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
 244     fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
 245     fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
 246     fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
 247
 248     fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
 249     fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
 250
 251     Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) + GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity);
 252     fHistRemovedEnergy->Fill(removedEnergy);
 253
 254     fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
 255     fTotEt = fTotChargedEt + fTotNeutralEt;
 256 // Fill the histograms...0
 257     FillHistograms();
 258     //std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
 259     return 0;
 260 }
 261
 262 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
 263 { // check vertex
 264
 265     Float_t bxy = 999.;
 266     Float_t bz = 999.;
 267     if (!track) {
 268         AliError("ERROR: no track");
 269         return kFALSE;
 270     }
 271     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 272     if (!esdTrack) {
 273         AliError("ERROR: no track");
 274         return kFALSE;
 275     }
 276     esdTrack->GetImpactParametersTPC(bxy,bz);
 277
 278
 279     bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
 280                   (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
 281                   (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
 282                   (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
 283                   (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
 284
 285     return status;
 286 }
 287
 288 void AliAnalysisEtReconstructed::Init()
 289 { // Init
 290     AliAnalysisEt::Init();
 291     fPidCut = fCuts->GetReconstructedPidCut();
 292     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
 293     if (!fCorrections) {
 294         cout<<"Warning!  You have not set corrections.  Your code will crash.  You have to set the corrections."<<endl;
 295     }
 296 }
 297
 298 bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
 299 { // propagate track to detector radius
 300
 301     if (!track) {
 302         cout<<"Warning: track empty"<<endl;
 303         return kFALSE;
 304     }
 305     AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track);
 306     if (!esdTrack) {
 307         AliError("ERROR: no ESD track");
 308         return kFALSE;
 309     }
 310     // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 311
 312     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 313
 314     // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 315     return prop && fSelector->CutGeometricalAcceptance(*esdTrack);
 316 }
 317
 318 void AliAnalysisEtReconstructed::FillOutputList(TList* list)
 319 { // add some extra histograms to the ones from base class
 320     AliAnalysisEt::FillOutputList(list);
 321
 322     list->Add(fHistChargedPionEnergyDeposit);
 323     list->Add(fHistProtonEnergyDeposit);
 324     list->Add(fHistAntiProtonEnergyDeposit);
 325     list->Add(fHistChargedKaonEnergyDeposit);
 326     list->Add(fHistMuonEnergyDeposit);
 327
 328     list->Add(fHistRemovedEnergy);
 329     list->Add(fClusterPosition);
 330     list->Add(fClusterEnergy);
 331     list->Add(fClusterEt);
 332
 333     list->Add(fHistChargedEnergyRemoved);
 334     list->Add(fHistNeutralEnergyRemoved);
 335     list->Add(fHistGammaEnergyAdded);
 336 }
 337
 338 void AliAnalysisEtReconstructed::CreateHistograms()
 339 { // add some extra histograms to the ones from base class
 340     AliAnalysisEt::CreateHistograms();
 341
 342     Int_t nbinsEt = 1000;
 343     Double_t minEt = 0;
 344     Double_t maxEt = 10;
 345
 346     // possibly change histogram limits
 347 //     if (fCuts) {
 348 //         nbinsEt = fCuts->GetHistNbinsParticleEt();
 349 //         minEt = fCuts->GetHistMinParticleEt();
 350 //         maxEt = fCuts->GetHistMaxParticleEt();
 351 //     }
 352
 353     TString histname;
 354     histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
 355     fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 356     fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 357     fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
 358
 359     histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
 360     fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 361     fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 362     fHistProtonEnergyDeposit->SetYTitle("Energy of track");
 363
 364     histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
 365     fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 366     fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 367     fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
 368
 369     histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
 370     fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 371     fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 372     fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
 373
 374     histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
 375     fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
 376     fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 377     fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 378
 379     histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
 380     fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
 381     //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
 382     //fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 383
 384     histname = "fClusterPosition" + fHistogramNameSuffix;
 385     fClusterPosition = new TH2D(histname.Data(), "Position of accepted neutral clusters",1000, -2.0, -.5, 1000, -.13 , 0.13);
 386     fClusterPosition->SetXTitle("Energy deposited in calorimeter");
 387     fClusterPosition->SetYTitle("Energy of track");
 388
 389     histname = "fClusterEnergy" + fHistogramNameSuffix;
 390     fClusterEnergy = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
 391     fClusterEnergy->SetXTitle("Number of clusters");
 392     fClusterEnergy->SetYTitle("Energy of cluster");
 393
 394     histname = "fClusterEt" + fHistogramNameSuffix;
 395     fClusterEt = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
 396     fClusterEt->SetXTitle("Number of clusters");
 397     fClusterEt->SetYTitle("E_{T} of cluster");
 398
 399     histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
 400     fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 401
 402     histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
 403     fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 404
 405     histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
 406     fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 407
 408
 409 }