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