48bbfc8319f4356969f400ced46bd0f2536d8241
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidEMCAL.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, 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 // Class for EMCAL PID
17 // Implements the abstract base class AliHFEpidBase
18 // IsInitialized() does the PID decision with energy and 
19 // momentum matching (e/p) 
20 // 
21 // Authors:
22 //  Shingo Sakai 
23 //   
24 //
25 //
26 #include <TMath.h>
27 #include "AliESDInputHandler.h"
28 #include "AliAODpidUtil.h"
29 #include "AliAODTrack.h"
30 #include "AliESDpid.h"
31
32 #include "AliHFEdetPIDqa.h"
33 //#include "AliHFEpidEMCAL.h"
34 #include "AliHFEpidQAmanager.h"
35
36 //#include "AliEMCALRecoUtils.h"
37 //#include "AliEMCALGeometry.h"
38 #include "AliVCluster.h"
39 #include "AliVCaloCells.h"
40 #include "AliVEvent.h"
41 #include "AliLog.h"
42 #include "AliEMCALPIDUtils.h"
43 #include "AliPID.h"
44 #include "AliESDEvent.h"
45 #include "AliESDtrack.h"
46 //#include "AliEMCALTrack.h"
47 //#include "AliEMCALTracker.h"
48 #include "AliHFEpidEMCAL.h"
49
50 ClassImp(AliHFEpidEMCAL)
51
52 //___________________________________________________________________
53 AliHFEpidEMCAL::AliHFEpidEMCAL():
54   AliHFEpidBase()
55   , fPID(NULL)
56   , feopMim(0.8)
57   , feopMax(1.4)
58 {
59   //
60   // Constructor
61   //
62
63
64 //___________________________________________________________________
65 AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name):
66   AliHFEpidBase(name)
67   , fPID(NULL)
68   , feopMim(0.8)
69   , feopMax(1.4)
70 {
71   //
72   // Constructor
73   //
74 }
75
76 //___________________________________________________________________
77 AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
78   AliHFEpidBase("")
79   , fPID(NULL)
80   , feopMim(0.8)
81   , feopMax(1.4)
82 {  
83   // 
84   // Copy operator
85   //
86
87   c.Copy(*this);
88 }
89 //___________________________________________________________________
90 AliHFEpidEMCAL &AliHFEpidEMCAL::operator=(const AliHFEpidEMCAL &ref){
91   //
92   // Assignment operator
93   //
94
95   if(this != &ref){
96     ref.Copy(*this);
97   }
98
99   return *this;
100 }
101 //___________________________________________________________________
102 AliHFEpidEMCAL::~AliHFEpidEMCAL(){
103   //
104   // Destructor
105   //
106   if(fPID) delete fPID;
107 }
108 //___________________________________________________________________
109 void AliHFEpidEMCAL::Copy(TObject &ref) const {
110   //
111   // Performs the copying of the object
112   //
113   AliHFEpidEMCAL &target = dynamic_cast<AliHFEpidEMCAL &>(ref);
114
115   target.fPID = fPID;          
116   target.feopMax = feopMax;
117   target.feopMim = feopMim;
118
119   AliHFEpidBase::Copy(ref);
120 }
121 //___________________________________________________________________
122 Bool_t AliHFEpidEMCAL::InitializePID(){
123   //
124   // InitializePID: EMCAL experts have to implement code here
125   //
126   return kTRUE;
127 }
128
129 //___________________________________________________________________
130 Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager */*pidqa*/) const
131 { //Function to a return a code indicating whether or not an electron candidate is selected
132   //
133   //
134   if((!fESDpid && track->IsESDanalysis()) || (!fAODpid && track->IsAODanalysis())) return 0;
135   AliDebug(2, "PID object available");
136   // EMCal not fESDpid  (s.s Feb. 11)
137
138   /*
139   AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
140   if(!((dynamic_cast<const AliVTrack *>(track->GetRecTrack()))->GetStatus() & AliESDtrack::kEMCALpid)) return 0;
141   AliDebug(2, "Track Has EMCAL PID");
142   */
143   //if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
144   // not QA for EMCal, will be added (s.s Feb. 11)
145
146   Double_t eop = 0.;//MomentumEnergyMatch(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
147   AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
148   Int_t pdg = 0;
149   if(eop>feopMim && eop<feopMax){
150     pdg = 11;
151     //if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
152   }
153   else
154   {
155    pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
156   } 
157
158     printf("eID %g ; %d \n",eop,pdg);  
159
160   return pdg;
161 }
162 /*
163 Double_t AliHFEpidEMCAL::MomentumEnergyMatch(const AliVParticle *track) const
164 { // Returns e/p if an electron is matched
165
166   Float_t  clsPos[3];
167   Double_t trkPos[3];
168   Double_t matchclsE = 9999.9;
169
170   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
171   AliESDEvent *evt = esdtrack->GetESDEvent();
172   Double_t  magF = evt->GetMagneticField();
173   Double_t magSign = 1.0;
174   if(magF<0)magSign = -1.0;
175   //printf("magF ; %g ; %g \n", magF,magSign);
176  
177   if (!TGeoGlobalMagField::Instance()->GetField()) {
178           printf("Loading field map...\n");
179           //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
180           AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
181           TGeoGlobalMagField::Instance()->SetField(field);
182                     }
183
184   AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
185   Double_t fieldB[3]; 
186   emctrack->GetBxByBz(fieldB);
187   //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
188
189   for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
190    {
191
192    AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
193    if(!cluster->IsEMCAL()) continue;
194
195    cluster->GetPosition(clsPos);
196    if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
197    emctrack->GetXYZ(trkPos);
198
199    TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
200    TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
201
202    Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi();  // track cluster matching
203    Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
204
205    double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
206
207    if(rmatch<0.02)
208       {
209        matchclsE = cluster->E();
210       }
211   }
212    delete emctrack;
213
214    double feop = -9999.9;
215    if(matchclsE<9999) feop = matchclsE/esdtrack->P();
216
217    //   if(feop!=-9999.9)printf("%f\n",feop) ; 
218    return feop;
219
220 }
221 */
222