318bf5a7cbc059e8293f009b5f0dbd78f6695646
[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 "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 }