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