1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Class AliHFEemcalPIDqa
18 // Monitoring EMCAL PID in the HFE PID montioring framework. The following
19 // quantities are monitored:
20 // TPC Signal distribution
21 // (Always as function of momentum, particle species and centrality
22 // before and after cut)
23 // More information about the PID monitoring framework can be found in
24 // AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
28 // Shingo Sakai <ssakai@lbl.gov>
31 #include <THnSparse.h>
35 #include "AliESDInputHandler.h"
36 //#include "AliVCluster.h"
37 //#include "AliVCaloCells.h"
38 //#include "AliVEvent.h"
40 #include "AliPIDResponse.h"
44 #include "AliVParticle.h"
45 //#include "AliVTrack.h"
46 //#include "AliESDtrack.h"
47 #include "AliHFEcollection.h"
48 #include "AliHFEpid.h"
49 #include "AliHFEpidBase.h"
50 #include "AliHFEpidQAmanager.h"
51 #include "AliHFEpidTPC.h"
52 #include "AliHFEpidEMCAL.h"
53 #include "AliHFEtools.h"
54 #include "AliHFEemcalPIDqa.h"
55 #include "AliTrackerBase.h"
57 ClassImp(AliHFEemcalPIDqa)
59 //_________________________________________________________
60 AliHFEemcalPIDqa::AliHFEemcalPIDqa():
69 //_________________________________________________________
70 AliHFEemcalPIDqa::AliHFEemcalPIDqa(const char* name):
71 AliHFEdetPIDqa(name, "QA for EMCAL")
75 // Default constructor
79 //_________________________________________________________
80 AliHFEemcalPIDqa::AliHFEemcalPIDqa(const AliHFEemcalPIDqa &o):
90 //_________________________________________________________
91 AliHFEemcalPIDqa &AliHFEemcalPIDqa::operator=(const AliHFEemcalPIDqa &o){
95 AliHFEdetPIDqa::operator=(o);
96 if(&o != this) o.Copy(*this);
100 //_________________________________________________________
101 AliHFEemcalPIDqa::~AliHFEemcalPIDqa(){
105 if(fHistos) delete fHistos;
108 //_________________________________________________________
109 void AliHFEemcalPIDqa::Copy(TObject &o) const {
113 AliHFEemcalPIDqa &target = dynamic_cast<AliHFEemcalPIDqa &>(o);
115 delete target.fHistos;
116 target.fHistos = NULL;
118 if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
121 //_________________________________________________________
122 Long64_t AliHFEemcalPIDqa::Merge(TCollection *coll){
124 // Merge with other objects
127 if(coll->IsEmpty()) return 1;
130 AliHFEemcalPIDqa *refQA = NULL;
135 refQA = dynamic_cast<AliHFEemcalPIDqa *>(o);
138 listHistos.Add(refQA->fHistos);
141 fHistos->Merge(&listHistos);
145 //_________________________________________________________
146 void AliHFEemcalPIDqa::Initialize(){
151 fHistos = new AliHFEcollection("emcalqahistos", "Collection of EMCAL QA histograms");
153 // Make common binning
154 const Int_t kCentralityBins = 11;
155 const Double_t kMinP = 0.;
156 const Double_t kMaxP = 50.;
157 const Double_t kTPCSigMim = 40.;
158 const Double_t kTPCSigMax = 140.;
160 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
161 Int_t nBins[16] = {AliPID::kSPECIES + 1, 500, 500, 400, 300, 100, 600, 300, 100, 40, 300, 300, 300, kCentralityBins, 6, 2};
162 Double_t min[16] = {-1, kMinP, kMinP, kTPCSigMim, 0., -1.0, -8.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -3.0, 0.};
163 Double_t max[16] = {AliPID::kSPECIES, kMaxP, kMaxP, kTPCSigMax, 6.0, 1.0, 4.0, 3.0, 0.1, 40, 3.0, 3.0, 3.0, 11., 3.0, 2.};
164 fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; phi ; eta ; nSig ; E/p ; Rmatch ; Ncell ; M02 ; M20 ; Disp ; Centrality;charge ;PID Step; ", 16, nBins, min, max);
166 //2nd histogram: EMCAL signal - E/p
167 //Int_t nBins2[7] = {AliPID::kSPECIES + 1, 500, 500, 500, 100, 600, 2};
168 //Double_t min2[7] = {-1, kMinP, kMinP, 0, 0, -8., 0.};
169 //Double_t max2[7] = {AliPID::kSPECIES, kMaxP, kMaxP, 5, 0.1, 4., 2.};
170 //fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; TPCnsigma; PID Step", 7, nBins2, min2, max2);
177 //_________________________________________________________
178 void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t step){
180 // Fill TPC histograms
182 //AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
183 Float_t centrality = track->GetCentrality();
185 //const AliVTrack *vtrack = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
186 //const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(vtrack);
187 const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track->GetRecTrack());
189 Int_t species = track->GetAbInitioPID();
190 if(species >= AliPID::kSPECIES) species = -1;
192 // Get TPC nSigma (based on AliHFEtpcPIDqa)
193 AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
194 const AliPIDResponse *pidResponse = tpcpid ? tpcpid->GetPIDResponse() : NULL;
195 Float_t nSigmatpc = pidResponse ? pidResponse->NumberOfSigmasTPC(static_cast<const AliVTrack *>(track->GetRecTrack()), AliPID::kElectron) : 0.;
196 //printf("nSigmatpc = %f\n",nSigmatpc);
197 AliDebug(2, Form("nSigmatpc = %f\n",nSigmatpc));
200 double charge = esdtrack->Charge();
201 //printf("charge %f\n",charge);
204 //printf("phi %f, eta %f\n;",phi,eta);
207 Double_t showerEMC[4];
208 TVector3 emcsignal = MomentumEnergyMatchV2(esdtrack,showerEMC);
209 AliDebug(2, Form("shower shape ; %f %f %f %f \n",showerEMC[0],showerEMC[1],showerEMC[2],showerEMC[3]));
212 Double_t contentSignal2[7];
213 contentSignal2[0] = species;
214 contentSignal2[1] = track->GetRecTrack()->P();
215 contentSignal2[2] = track->GetRecTrack()->Pt();
216 contentSignal2[3] = emcsignal(0);//e over p
217 contentSignal2[4] = emcsignal(1);//residual
218 contentSignal2[5] = nSigmatpc;
219 contentSignal2[6] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
222 Double_t contentSignal[16];
223 contentSignal[0] = species;
224 contentSignal[1] = track->GetRecTrack()->P();
225 contentSignal[2] = track->GetRecTrack()->Pt();
226 contentSignal[3] = esdtrack->GetTPCsignal(); //?
227 double phi = track->GetRecTrack()->Phi();
228 double eta = track->GetRecTrack()->Eta();
229 contentSignal[4] = phi;
230 contentSignal[5] = eta;
231 contentSignal[6] = nSigmatpc;
232 contentSignal[7] = emcsignal(0);
233 contentSignal[8] = emcsignal(2);
234 contentSignal[9] = showerEMC[0];
235 contentSignal[10] = showerEMC[1];
236 contentSignal[11] = showerEMC[2];
237 contentSignal[12] = showerEMC[3];
238 contentSignal[13] = centrality;
239 contentSignal[14] = charge;
240 contentSignal[15] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
243 //printf("ProcessTrack ; Print Content %g; %g; %g; %g \n",contentSignal[0],contentSignal[1],contentSignal[2],contentSignal[3]);
244 fHistos->Fill("EMCAL_TPCdedx", contentSignal);
245 //fHistos->Fill("EMCAL_Signal", contentSignal2);
248 //_________________________________________________________
249 TH1 *AliHFEemcalPIDqa::GetHistogram(const char *name){
253 if(!fHistos) return NULL;
254 TObject *histo = fHistos->Get(name);
255 if(!histo->InheritsFrom("TH1")) return NULL;
256 return dynamic_cast<TH1 *>(histo);
260 //___________________________________________________________________________
261 TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV2(const AliESDtrack *esdtrack, Double_t *shower) const
262 { // Returns e/p & Rmatch
266 Double_t fMass=0.139;
267 Double_t fStep=10; //This is taken from EMCAL tender, hardcoded!
269 TVector3 refVec(-9999,-9999,-9999);
270 Double_t matchclsE = 9999.9;
271 for(int i=0; i<4; i++)shower[i] = -9999;
273 const AliESDEvent *evt = esdtrack->GetESDEvent();
275 //Proper trackmatching, added by Tomas
276 //Get matched cluster
277 Int_t icl = esdtrack->GetEMCALcluster(); //From tender
278 AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
279 if(cluster==NULL) return refVec;
280 if(!cluster->IsEMCAL()) return refVec;
281 //printf("EMCal cluster pointer ? %p\n",(void*)cluster);
282 //printf("EMCal cluster ? %d\n",cluster->IsEMCAL());
285 double delphi = cluster->GetTrackDx();
286 double deleta = cluster->GetTrackDz();
287 double rmatch1 = sqrt(pow(delphi,2)+pow(deleta,2));
289 matchclsE = cluster->E();
291 //double feop = -9999.9;
293 double feop = matchclsE/esdtrack->P();
295 shower[0] = cluster->GetNCells(); // number of cells in cluster
296 shower[1] = cluster->GetM02(); // long axis
297 shower[2] = cluster->GetM20(); // short axis
298 shower[3] = cluster->GetDispersion(); // dispersion
300 // if(feop!=-9999.9)printf("%f\n",feop) ;
301 TVector3 emcsignal(feop,0.0,rmatch1);