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