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