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 "AliESDCaloCluster.h"
  14 #include "TVector3.h"
  15 #include "AliVEvent.h"
  16 #include "AliESDEvent.h"
  17 #include "AliVParticle.h"
  18 #include <iostream>
  19 #include "TH2F.h"
  20
  21 AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
  22         AliAnalysisEt()
  23         ,fTrackDistanceCut(0)
  24         ,fClusterType(0)
  25 {
  26
  27 }
  28
  29 AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
  30 {
  31 }
  32
  33 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
  34 { // analyse ESD event
  35     ResetEventValues();
  36     AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
  37
  38     for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++)
  39     {
  40         AliVParticle *track = event->GetTrack(iTrack);
  41         if (!track)
  42         {
  43             Printf("ERROR: Could not get track %d", iTrack);
  44             continue;
  45         }
  46
  47         fMultiplicity++;
  48
  49         Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0);
  50         Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1);
  51
  52         Float_t massPart = 0;
  53
  54         const Double_t *pidWeights = track->PID();
  55         if (pidWeights)
  56         {
  57             Int_t maxpid = -1;
  58             Float_t maxpidweight = 0;
  59             for (Int_t p =0; p < AliPID::kSPECIES; p++)
  60             {
  61                 if (pidWeights[p] > maxpidweight)
  62                 {
  63                     maxpidweight = pidWeights[p];
  64                     maxpid = p;
  65                 }
  66             }
  67             if (maxpid == AliPID::kProton)
  68             {
  69                 //     massPart = -0.938*track->Charge();
  70             }
  71
  72         }
  73
  74         Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
  75         // 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
  76
  77         if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() )
  78         {
  79             fTotChargedEt +=  et;
  80             fChargedMultiplicity++;
  81
  82             if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
  83             {
  84                 fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
  85             }
  86         }
  87
  88         if (TrackHitsCalorimeter(track, event->GetMagneticField()))
  89         {
  90           Double_t phi = track->Phi();
  91           Double_t pt = track->Pt();
  92           // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
  93           if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
  94           else fHistPhivsPtNeg->Fill(phi, pt);
  95         }
  96     }
  97
  98     for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
  99     {
 100         AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
 101         if (!cluster)
 102         {
 103             Printf("ERROR: Could not get cluster %d", iCluster);
 104             continue;
 105         }
 106
 107         // printf("Rec Cluster: iCluster %03d E %4.3f type %d NCells %d\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells()); // tmp/debug printout
 108         if (cluster->GetType() != fClusterType) continue;
 109
 110         if (cluster->E() < fClusterEnergyCut) continue;
 111         Float_t pos[3];
 112         TVector3 cp(pos);
 113         cluster->GetPosition(pos);
 114         //if (pos[0] < -(32.0*2.2)) continue; //Ensure that modules 0 and 1 are not used
 115         // if(cp.Phi() < 260.*TMath::Pi()/180.) continue;
 116         fHistTMDeltaR->Fill(cluster->GetEmcCpvDistance());
 117         if (cluster->GetEmcCpvDistance() < fTrackDistanceCut)
 118         {
 119             continue;
 120             //AliVParticle *matchedTrack = event->GetTrack(cluster->GetTrackMatched());
 121 //          if(CheckGoodVertex(matchedTrack))
 122 //          {
 123 //             totChargedEnergy +=  matchedTrack->E();;
 124 //             totChargedEt += matchedTrack->E()*TMath::Sin(matchedTrack);
 125 //          }
 126         }
 127
 128         if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
 129
 130         cluster->GetPosition(pos);
 131
 132         // TODO: replace with TVector3, too lazy now...
 133
 134         float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 135
 136         float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 137         // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
 138         fTotNeutralEt += cluster->E() * TMath::Sin(theta);
 139         fNeutralMultiplicity++;
 140
 141         fMultiplicity++;
 142     }
 143
 144     fTotNeutralEtAcc = fTotNeutralEt;
 145
 146     fTotEt = fTotChargedEt + fTotNeutralEt;
 147     fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
 148
 149     std::cout << fTotChargedEtAcc << std::endl;
 150     // Fill the histograms...
 151     FillHistograms();
 152
 153     return 0;
 154 }
 155
 156 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
 157 { // check vertex
 158
 159     Float_t bxy = 999.;
 160     Float_t bz = 999.;
 161     dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz);
 162
 163     bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
 164       (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
 165       (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
 166       (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
 167       (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
 168
 169     return status;
 170 }
 171
 172 void AliAnalysisEtReconstructed::Init()
 173 { // Init
 174     AliAnalysisEt::Init();
 175 }
 176
 177 bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
 178 { // propagate track to detector radius
 179
 180    AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 181    // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 182
 183     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 184
 185     // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 186     return prop &&
 187                    TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
 188                    esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
 189                    esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
 190 }
 191