Transition to NewIO
[u/mrichter/AliRoot.git] / PHOS / AliPHOSFastGlobalReconstruction.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 // --- AliRoot header files ---
19
20 // Fast global reconstruction class.
21 // It performes fast reconstruction for charged particles only,
22 // assuming that they were detected by all central ALICE detectors but PHOS.
23 // This class acts as a filter for primary particles, selects them
24 // and deteriorates their 4-momenta.
25 // The filter can recognize the primary particle list from different
26 // event generators: Pythia, Hijing, HIJINGPara.
27 //!
28 // Usage:
29 // rec = new AliPHOSGlobalReconstruction("galice.root");
30 // rec->FastReconstruction(ievent);
31 // TClonesArray *recList = rec->GetRecParticles();
32 //!
33 // See also $ALICE_ROOT/macros/TestGlobalReconstruction.C
34 //!
35 // Author: Yuri Kharlov. 17 April 2003
36
37 #include "AliRun.h"
38 #include "AliGenerator.h"
39 #include "AliPHOSGetter.h"
40 #include "AliPHOSFastGlobalReconstruction.h"
41
42 ClassImp(AliPHOSFastGlobalReconstruction)
43
44 //____________________________________________________________________________
45 AliPHOSFastGlobalReconstruction::AliPHOSFastGlobalReconstruction(const char* headerFile )
46 {
47   // Constructor of fast global reconstruction:
48   // create an instance of the PHOS getter,
49   // create an array or reconstructed particles.
50
51   gime = AliPHOSGetter::Instance(headerFile);
52   fGenerator = gAlice->Generator();
53   fParticles = new TClonesArray("TParticle",100);
54   fNParticles = 0;
55 }
56
57 //____________________________________________________________________________
58 AliPHOSFastGlobalReconstruction::~AliPHOSFastGlobalReconstruction()
59 {
60   // Destructor of fast global reconstruction:
61   // delete the array of reconstructed particles
62   if (fParticles != 0) {
63     delete fParticles;
64     fParticles  = 0;
65     fNParticles = 0;
66   }
67 }
68
69 //____________________________________________________________________________
70 void AliPHOSFastGlobalReconstruction::FastReconstruction(Int_t event)
71 {
72   // Perform a fast global reconstruction of event numbered "event".
73   // Reconstructed particles will be stored into array recParticles
74
75   TParticle *primary;
76   TLorentzVector p,v;
77   Int_t kf,ks,imom1,imom2,idaug1,idaug2;
78
79   gime->Event(event,"X") ;
80   fParticles  ->Clear();
81   fNParticles = 0;
82   Int_t        nPrimaries = gime->NPrimaries();
83   TClonesArray *primaries = gime->Primaries();
84
85   for (Int_t iprim=0; iprim<nPrimaries; iprim++) {
86     primary = (TParticle*)primaries->At(iprim);
87     if ((strcmp(fGenerator->GetName(),"Pythia")    ==0 && primary->GetStatusCode() == 1) ||
88         (strcmp(fGenerator->GetName(),"HIJINGpara")==0 && primary->GetFirstMother()==-1) ||
89         (strcmp(fGenerator->GetName(),"Hijing")    ==0 && primary->GetStatusCode() == 3)) {
90       if (Detected(primary)) {
91         primary->Momentum(p);
92         primary->ProductionVertex(v);
93         kf     = primary->GetPdgCode();
94         ks     = primary->GetStatusCode();
95         imom1  = primary->GetFirstMother();
96         imom2  = primary->GetSecondMother();
97         idaug1 = primary->GetFirstDaughter();
98         idaug2 = primary->GetLastDaughter();
99         SmearMomentum(p);
100         new((*fParticles)[fNParticles]) TParticle(kf,ks,imom1,imom2,idaug1,idaug2,p,v);
101         fNParticles++;
102       }
103     }
104   }
105 }
106
107 //____________________________________________________________________________
108 Bool_t AliPHOSFastGlobalReconstruction::Detected(TParticle *particle)
109 {
110   // Returns kTRUE is a particle is reconstructed, kFALSE otherwise.
111   // A particle is reconstructed if it is charged and accepted with the
112   // probability Efficiency(pt,eta) depending on pt and eta.
113
114   Bool_t detected = kFALSE;
115   if (particle->GetPDG()->Charge() != 0) {
116     Float_t pt  = particle->Pt();
117     Float_t eta = particle->Eta();
118     if (gRandom->Rndm() < Efficiency(pt,eta)) detected = kTRUE;
119   }
120   return detected;
121 }
122
123 //____________________________________________________________________________
124 Float_t AliPHOSFastGlobalReconstruction::Efficiency(Float_t pt, Float_t eta)
125 {
126   // Detection probability vs. pt and eta, i.e. a probability to detect
127   // a particle with transverse momentum pt and speudorapidity eta.
128   // For the moment assume that charged particles are detected with
129   // 80% efficiency within |eta|<0.9 and pt>0.15 GeV, and with 0% efficiency
130   // beyond that acceptance
131
132   const Float_t kEtaLimit = 0.9;
133   const Float_t kPtLimit  = 0.15;
134   Float_t efficiency = 0.0;
135   if (TMath::Abs(eta) < kEtaLimit && pt > kPtLimit) efficiency = 0.80;
136   return efficiency;
137 }
138
139 //____________________________________________________________________________
140 void AliPHOSFastGlobalReconstruction::SmearMomentum(TLorentzVector &p)
141 {
142   // Smear 4-momentum according to known resolution (2% for the moment)
143
144   const Float_t kAngleResolution    = 0.02;
145   const Float_t kMomentumResolution = 0.02;
146   Double_t mass = p.M();
147   for (Int_t i=0; i<3; i++) {
148     p[i] *= gRandom->Gaus(1.,kAngleResolution   );
149     p[i] *= gRandom->Gaus(1.,kMomentumResolution);
150   }
151   p[3] = TMath::Sqrt(p.P()*p.P() + mass*mass);
152 }