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