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