407aa5a44d4b046ee46a52d4b1397b8bf4cdb540
[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         ,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