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