]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PHOS/AliPHOSAnalyze.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Timing added to reconstruction
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.cxx
CommitLineData
6ad0bfa0 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
b2a60966 16/* $Id$ */
17
6ad0bfa0 18//_________________________________________________________________________
efad3788 19// Algorythm class to analyze PHOS events. In this class we demostrate,
20// how to handle reconstructed objects with AliPHSOIndexToObject.
21// As an example we propose sulotions for four most frequently used tasks:
22// DrawRecon(...) - to draw reconstructed objects in the PHOS plane,
23// very usefull in the debuging
24// InvarianMass(...) - to calculate "REAL" and "MIXED" photon pairs
25// invariant mass distributions
26// EnergyResoluition(...) -\ Energy and position resolutions of the
27// PositionResolution(...)-/ reconstructed photons
28// Contamination(...) - calculates contamination of the photon spectrum and
29// pobability of reconstruction of several primaries as
30// kGAMMA,kELECTRON etc.
31//
32// User Case:
e4761a49 33// root [0] AliPHOSAnalyze * a = new AliPHOSAnalyze("galice.root")
efad3788 34// // set the file you want to analyse
35// root [1] a->DrawRecon(1,3)
36// // plot RecObjects, made in event 1, PHOS module 3
37// root [2] a->DrawRecon(1,3,"PHOSRP","another PID")
38// // plot RecObjets made in the event 1, PHOS module 3,
39// // produced in the another reconstruction pass,
40// // which produced PHOS RecParticles ("PHOSRP") with
41// // title "another PID".
42// root [3] a->InvariantMass()
43// // Calculates "REAL" and "MIXED" invariant mass
44// // distributions of kGAMMA and (kGAMMA+kNEUTRALEM)
45// // and APPENDS this to the file "invmass.root"
46// root [4] a->PositionResolution()
47// // calculates two dimentional histos: energy of the primary
48// // photon vs distance betwin incedence point and reconstructed
49// // poisition. One can analyse the produced file position.root
50// // with macro PhotonPosition.C
51// root [5] a->EnergyResolution()
52// // calculates two dimentional histos: energy of the primary
53// // photon vs energy of the reconstructed particle. One can
54// // analyse the produced file energy.root
55// // with macro PhotonEnergy.C
56// root [6] a->Contamination()
57// // fills spectra of primary photons and several kinds of
58// // reconstructed particles, so that analyzing them one can
59// // estimate conatmination, efficiency of registration etc.
a3dfe79c 60//*--
efad3788 61//*-- Author: Dmitri Peressounko (SUBATECH & RRC Kurchatov Institute)
6ad0bfa0 62//////////////////////////////////////////////////////////////////////////////
63
baef0810 64
6ad0bfa0 65// --- ROOT system ---
66
92862013 67#include "TFile.h"
68#include "TH1.h"
6ad0bfa0 69#include "TH2.h"
83448140 70#include "TH2.h"
6ad0bfa0 71#include "TParticle.h"
72#include "TClonesArray.h"
73#include "TTree.h"
74#include "TMath.h"
a88eea6f 75#include "TROOT.h"
76#include "TFolder.h"
6ad0bfa0 77
78// --- Standard library ---
79
efad3788 80#include <iomanip.h>
92862013 81
6ad0bfa0 82// --- AliRoot header files ---
83
84#include "AliRun.h"
83448140 85#include "AliPHOSv1.h"
6ad0bfa0 86#include "AliPHOSAnalyze.h"
6ad0bfa0 87#include "AliPHOSDigit.h"
efad3788 88#include "AliPHOSSDigitizer.h"
6ad0bfa0 89#include "AliPHOSTrackSegment.h"
90#include "AliPHOSRecParticle.h"
867ede0e 91#include "AliPHOSCpvRecPoint.h"
7b7c1533 92#include "AliPHOSGetter.h"
6ad0bfa0 93
baef0810 94
6ad0bfa0 95ClassImp(AliPHOSAnalyze)
96
6ad0bfa0 97//____________________________________________________________________________
98 AliPHOSAnalyze::AliPHOSAnalyze()
99{
b2a60966 100 // default ctor (useless)
efad3788 101 fCorrection = 1.2 ; //Value calculated for default parameters of reconstruction
6ad0bfa0 102}
103
104//____________________________________________________________________________
efad3788 105AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName)
6ad0bfa0 106{
83974468 107 // ctor: analyze events from root file "name"
efad3788 108 ffileName = fileName ;
7b7c1533 109 fCorrection = 1.05 ; //Value calculated for default parameters of reconstruction
6ad0bfa0 110}
111
88714635 112//____________________________________________________________________________
113AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
114{
115 // copy ctor
116 ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
117}
118
6ad0bfa0 119//____________________________________________________________________________
120AliPHOSAnalyze::~AliPHOSAnalyze()
121{
122 // dtor
123
6ad0bfa0 124}
6ad0bfa0 125//____________________________________________________________________________
efad3788 126void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,const char* branchTitle){
127 //Draws pimary particles and reconstructed
128 //digits, RecPoints, RecPartices etc
129 //for event Nevent in the module Nmod.
9688c1dd 130
7b7c1533 131 //========== Create ObjectGetter
132 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
133 const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
134 gime->Event(Nevent);
9688c1dd 135
136 Int_t nx = phosgeom->GetNPhi() ;
137 Int_t nz = phosgeom->GetNZ() ;
138 Float_t * cri= phosgeom->GetEMCAGeometry()->GetCrystalHalfSize() ;
139 Float_t x = nx*cri[0] ;
140 Float_t z = nz*cri[2] ;
141 Int_t nxCPV = (Int_t) (nx*phosgeom->GetPadSizePhi()/(2.*cri[0])) ;
142 Int_t nzCPV = (Int_t) (nz*phosgeom->GetPadSizeZ()/(2.*cri[2])) ;
143
144 TH2F * emcDigits = new TH2F("emcDigits","EMC digits", nx,-x,x,nz,-z,z);
145 TH2F * emcSdigits = new TH2F("emcSdigits","EMC sdigits", nx,-x,x,nz,-z,z);
146 TH2F * emcRecPoints = new TH2F("emcRecPoints","EMC RecPoints",nx,-x,x,nz,-z,z);
147 TH2F * cpvSdigits = new TH2F("cpvSdigits","CPV sdigits", nx,-x,x,nz,-z,z);
148 TH2F * cpvDigits = new TH2F("cpvDigits","CPV digits", nxCPV,-x,x,nzCPV,-z,z) ;
149 TH2F * cpvRecPoints = new TH2F("cpvCl","CPV RecPoints", nxCPV,-x,x,nzCPV,-z,z) ;
150 TH2F * nbar = new TH2F("nbar","Primary nbar", nx,-x,x,nz,-z,z);
151 TH2F * phot = new TH2F("phot","Primary Photon", nx,-x,x,nz,-z,z);
152 TH2F * charg = new TH2F("charg","Primary charged",nx,-x,x,nz,-z,z);
153 TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",nx,-x,x,nz,-z,z);
154 TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", nx,-x,x,nz,-z,z);
155
ed4205d8 156
efad3788 157 //Plot Primary Particles
7b7c1533 158 const TParticle * primary ;
efad3788 159 Int_t iPrimary ;
7b7c1533 160 for ( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++)
efad3788 161 {
7b7c1533 162 primary = gime->Primary(iPrimary) ;
efad3788 163 Int_t primaryType = primary->GetPdgCode() ;
9688c1dd 164 if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212)
165 ||(primaryType == 11)||(primaryType == -11) ) {
efad3788 166 Int_t moduleNumber ;
167 Double_t primX, primZ ;
7b7c1533 168 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
efad3788 169 if(moduleNumber==Nmod)
170 charg->Fill(primZ,primX,primary->Energy()) ;
171 }
172 if( primaryType == 22 ) {
173 Int_t moduleNumber ;
174 Double_t primX, primZ ;
7b7c1533 175 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
176 if(moduleNumber==Nmod)
efad3788 177 phot->Fill(primZ,primX,primary->Energy()) ;
178 }
179 else{
180 if( primaryType == -2112 ) {
181 Int_t moduleNumber ;
182 Double_t primX, primZ ;
7b7c1533 183 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
efad3788 184 if(moduleNumber==Nmod)
185 nbar->Fill(primZ,primX,primary->Energy()) ;
186 }
187 }
188 }
9688c1dd 189
ed4205d8 190
efad3788 191 Int_t iSDigit ;
cce01dd5 192 AliPHOSDigit * sdigit ;
9688c1dd 193 TClonesArray * sdigits = gime->SDigits() ;
194 if(sdigits){
195 for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++)
efad3788 196 {
9688c1dd 197 sdigit = (AliPHOSDigit *) sdigits->At(iSDigit) ;
efad3788 198 Int_t relid[4];
7b7c1533 199 phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
efad3788 200 Float_t x,z ;
7b7c1533 201 phosgeom->RelPosInModule(relid,x,z) ;
9688c1dd 202 Float_t e = gime->SDigitizer()->Calibrate(sdigit->GetAmp()) ;
efad3788 203 if(relid[0]==Nmod){
204 if(relid[1]==0) //EMC
9688c1dd 205 emcSdigits->Fill(x,z,e) ;
206 if( relid[1]!=0 )
207 cpvSdigits->Fill(x,z,e) ;
efad3788 208 }
209 }
9688c1dd 210 }
211
efad3788 212 //Plot digits
213 Int_t iDigit ;
cce01dd5 214 AliPHOSDigit * digit ;
9688c1dd 215 TClonesArray * digits = gime->Digits();
216 if(digits) {
cce01dd5 217 for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++)
218 {
9688c1dd 219 digit = (AliPHOSDigit *) digits->At(iDigit) ;
cce01dd5 220 Int_t relid[4];
221 phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
222 Float_t x,z ;
223 phosgeom->RelPosInModule(relid,x,z) ;
9688c1dd 224 Float_t e = gime->SDigitizer()->Calibrate(digit->GetAmp()) ;
cce01dd5 225 if(relid[0]==Nmod){
226 if(relid[1]==0) //EMC
9688c1dd 227 emcDigits->Fill(x,z,e) ;
228 if( relid[1]!=0 )
229 cpvDigits->Fill(x,z,e) ;
cce01dd5 230 }
efad3788 231 }
9688c1dd 232 }
cce01dd5 233
234
efad3788 235 //Plot RecPoints
236 Int_t irecp ;
237 TVector3 pos ;
9688c1dd 238 TObjArray * emcrp = gime->EmcRecPoints() ;
239 if(emcrp) {
cce01dd5 240 for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){
9688c1dd 241 AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
cce01dd5 242 if(emc->GetPHOSMod()==Nmod){
243 emc->GetLocalPosition(pos) ;
9688c1dd 244 emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy());
cce01dd5 245 }
efad3788 246 }
9688c1dd 247 }
cce01dd5 248
9688c1dd 249 TObjArray * cpvrp = gime->CpvRecPoints() ;
250 if(cpvrp) {
cce01dd5 251 for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){
9688c1dd 252 AliPHOSRecPoint * cpv = (AliPHOSCpvRecPoint *) cpvrp->At(irecp) ;
253 if(cpv->GetPHOSMod()==Nmod){
254 cpv->GetLocalPosition(pos) ;
255 cpvRecPoints->Fill(pos.X(),pos.Z(),cpv->GetEnergy());
efad3788 256 }
257 }
9688c1dd 258 }
259
efad3788 260 //Plot RecParticles
9688c1dd 261 AliPHOSRecParticle * recParticle ;
efad3788 262 Int_t iRecParticle ;
9688c1dd 263 TClonesArray * rp = gime->RecParticles() ;
264 TClonesArray * ts = gime->TrackSegments() ;
265 if(rp && ts && emcrp) {
266 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ )
cce01dd5 267 {
9688c1dd 268 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
cce01dd5 269 Int_t moduleNumberRec ;
270 Double_t recX, recZ ;
271 phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
272 if(moduleNumberRec == Nmod){
ed4205d8 273
cce01dd5 274 Double_t minDistance = 5. ;
9688c1dd 275 Int_t closestPrimary = -1 ;
276
cce01dd5 277 //extract list of primaries: it is stored at EMC RecPoints
9688c1dd 278 Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
cce01dd5 279 Int_t numberofprimaries ;
9688c1dd 280 Int_t * listofprimaries = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
cce01dd5 281 Int_t index ;
282 const TParticle * primary ;
283 Double_t distance = minDistance ;
284
285 for ( index = 0 ; index < numberofprimaries ; index++){
286 primary = gime->Primary(listofprimaries[index]) ;
287 Int_t moduleNumber ;
288 Double_t primX, primZ ;
289 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
290 if(moduleNumberRec == moduleNumber)
291 distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
292 if(minDistance > distance)
293 {
294 minDistance = distance ;
295 closestPrimary = listofprimaries[index] ;
296 }
297 }
298
299 if(closestPrimary >=0 ){
300
301 Int_t primaryType = gime->Primary(closestPrimary)->GetPdgCode() ;
302
303 if(primaryType==22)
304 recPhot->Fill(recZ,recX,recParticle->Energy()) ;
305 else
306 if(primaryType==-2112)
307 recNbar->Fill(recZ,recX,recParticle->Energy()) ;
308 }
efad3788 309 }
310 }
9688c1dd 311
312 }
ed4205d8 313
efad3788 314 //Plot made histograms
9688c1dd 315 emcSdigits->Draw("box") ;
316 emcDigits->SetLineColor(5) ;
317 emcDigits->Draw("box") ;
318 emcRecPoints->SetLineColor(2) ;
319 emcRecPoints->Draw("boxsame") ;
320 cpvSdigits->SetLineColor(1) ;
321 cpvSdigits->Draw("box") ;
322 charg->SetLineColor(2) ;
323 charg->Draw("boxsame") ;
efad3788 324 nbar->SetLineColor(6) ;
325 nbar->Draw("boxsame") ;
2bed9e3e 326
2bed9e3e 327}
e4761a49 328//____________________________________________________________________________
329void AliPHOSAnalyze::Ls(){
330 //lists branches and titles of PHOS-related branches of TreeR, TreeD, TreeS
331
7b7c1533 332 AliPHOSGetter::GetInstance(ffileName.Data()) ;
134ce69a 333
e4761a49 334 Int_t ibranch;
335 TObjArray * branches;
336
337 branches = gAlice->TreeS()->GetListOfBranches() ;
338
339 cout << "TreeS: " << endl ;
340 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
341 TBranch * branch=(TBranch *) branches->At(ibranch) ;
342 if(strstr(branch->GetName(),"PHOS") )
343 cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ;
344 }
345 cout << endl ;
346
347 branches = gAlice->TreeD()->GetListOfBranches() ;
348
349 cout << "TreeD: " << endl ;
350 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
351 TBranch * branch=(TBranch *) branches->At(ibranch) ;
352 if(strstr(branch->GetName(),"PHOS") )
353 cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ;
354 }
355 cout << endl ;
356
357
358 branches = gAlice->TreeR()->GetListOfBranches() ;
359
360 cout << "TreeR: " << endl ;
361 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
362 TBranch * branch=(TBranch *) branches->At(ibranch) ;
363 if(strstr(branch->GetName(),"PHOS") )
364 cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ;
365 }
366 cout << endl ;
367
368
369}
2bed9e3e 370//____________________________________________________________________________
efad3788 371 void AliPHOSAnalyze::InvariantMass(const char* branchTitle)
2bed9e3e 372{
efad3788 373 // Calculates Real and Mixed invariant mass distributions
7b7c1533 374 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
8c0cd6e9 375
efad3788 376 Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
2bed9e3e 377
efad3788 378 //opening file
379 TFile * mfile = new TFile("invmass.root","update");
69183710 380
efad3788 381 //========== Reading /Booking Histograms
382 TH2D * hRealEM = 0 ;
383 hRealEM = (TH2D*) mfile->Get("hRealEM") ;
384 if(hRealEM == 0)
385 hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ;
386 TH2D * hRealPhot = 0 ;
387
388 hRealPhot = (TH2D*)mfile->Get("hRealPhot");
389 if(hRealPhot == 0)
390 hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
391
392 TH2D * hMixedEM = 0 ;
393 hMixedEM = (TH2D*) mfile->Get("hMixedEM") ;
394 if(hMixedEM == 0)
395 hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ;
396
397 TH2D * hMixedPhot = 0 ;
398 hMixedPhot = (TH2D*) mfile->Get("hMixedPhot") ;
399 if(hMixedPhot == 0)
400 hMixedPhot = new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
69183710 401
efad3788 402
403 //reading event and copyng it to TConesArray of all photons
404
405 TClonesArray * allRecParticleList = new TClonesArray("AliPHOSRecParticle", 1000) ;
fd0eb6da 406 Int_t * nRecParticles = new Int_t[nMixedEvents] ; // to mark boundaries of each event in the total list
efad3788 407 for(Int_t index = 0; index < nMixedEvents; index ++)
408 nRecParticles[index] = 0 ;
409 Int_t iRecPhot = 0 ; // number of EM particles in total list
69183710 410
efad3788 411 //scan over all events
412 Int_t event ;
154fc86d 413 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
7b7c1533 414 // for(event = 0; event < gime->MaxEvent(); event++ ){
415 for(event = 0; event < maxevent; event++ ){
9688c1dd 416 gime->Event(event,"R"); //will read only TreeR
efad3788 417
418 //copy EM RecParticles to the "total" list
7b7c1533 419 const AliPHOSRecParticle * recParticle ;
efad3788 420 Int_t iRecParticle ;
9688c1dd 421 TClonesArray * rp = gime->RecParticles() ;
422 if(!rp){
423 cout << "AliPHOSAnalyze::InvariantMass --> Can't find RecParticles " << endl ;
424 return ;
425 }
426
427 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ )
428 {
429 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
430 if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
431 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW))
432 new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
433 }
efad3788 434
435 Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle
436 nRecParticles[mevent] = iRecPhot-1 ;
9688c1dd 437
efad3788 438 //check, if it is time to calculate invariant mass?
154fc86d 439 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
7b7c1533 440 if((mevent == 0) && (event +1 == maxevent)){
9688c1dd 441
442 // if((mevent == 0) && (event +1 == gime->MaxEvent())){
69183710 443
efad3788 444 //calculate invariant mass:
445 Int_t irp1,irp2 ;
446 Int_t nCurEvent = 0 ;
69183710 447
efad3788 448 for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
449 AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
450
451 for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
452 AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
453
454 Double_t invMass ;
455 invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
456 (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
457 (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
458 (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
459
460 if(invMass> 0)
461 invMass = TMath::Sqrt(invMass);
462
463 Double_t pt ;
464 pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +
465 (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) );
466
467 if(irp1 > nRecParticles[nCurEvent])
468 nCurEvent++;
469
470 if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
471 hRealEM->Fill(invMass,pt);
9688c1dd 472 if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
473 (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
efad3788 474 hRealPhot->Fill(invMass,pt);
475 }
476 else{
477 hMixedEM->Fill(invMass,pt);
9688c1dd 478 if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
479 (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
efad3788 480 hMixedPhot->Fill(invMass,pt);
481 } //real-mixed
482
483 } //loop over second rp
484 }//loop over first rp
ed4205d8 485
efad3788 486 //Make some cleanings
487 for(Int_t index = 0; index < nMixedEvents; index ++)
488 nRecParticles[index] = 0 ;
489 iRecPhot = 0 ;
490 allRecParticleList->Clear() ;
69183710 491
efad3788 492 }
493 }
494 delete allRecParticleList ;
69183710 495
efad3788 496 //writing output
497 mfile->cd();
69183710 498
efad3788 499 hRealEM->Write(0,kOverwrite) ;
500 hRealPhot->Write(0,kOverwrite) ;
501 hMixedEM->Write(0,kOverwrite) ;
502 hMixedPhot->Write(0,kOverwrite) ;
69183710 503
efad3788 504 mfile->Write();
505 mfile->Close();
506 delete mfile ;
fd0eb6da 507 delete nRecParticles;
c3b9b3f9 508
509}
510
ed4205d8 511//____________________________________________________________________________
efad3788 512 void AliPHOSAnalyze::EnergyResolution(const char * branchTitle)
ed4205d8 513{
efad3788 514 //fills two dimentional histo: energy of primary vs. energy of reconstructed
ed4205d8 515
efad3788 516 TH2F * hAllEnergy = 0 ; //all reconstructed with primary photon
517 TH2F * hPhotEnergy= 0 ; //kGamma with primary photon
518 TH2F * hEMEnergy = 0 ; //electromagnetic with primary photon
ed4205d8 519
efad3788 520 //opening file and reading histograms if any
521 TFile * efile = new TFile("energy.root","update");
ed4205d8 522
efad3788 523 hAllEnergy = (TH2F*)efile->Get("hAllEnergy") ;
524 if(hAllEnergy == 0)
525 hAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
526
527 hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ;
528 if(hPhotEnergy == 0)
529 hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
530
531 hEMEnergy =(TH2F*) efile->Get("hEMEnergy");
532 if(hEMEnergy == 0)
533 hEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.);
ed4205d8 534
ed4205d8 535
7b7c1533 536 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
537 const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
efad3788 538
539 Int_t ievent;
154fc86d 540 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
7b7c1533 541 for ( ievent=0; ievent < maxevent ; ievent++){
efad3788 542
543 //read the current event
7b7c1533 544 gime->Event(ievent) ;
efad3788 545
7b7c1533 546 const AliPHOSRecParticle * recParticle ;
efad3788 547 Int_t iRecParticle ;
9688c1dd 548 TClonesArray * rp = gime->RecParticles() ;
549 if(!rp) {
550 cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
551 << ", Can't find RecParticles " << endl ;
552 return ;
553 }
554 TClonesArray * ts = gime->TrackSegments() ;
555 if(!ts) {
556 cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
557 << ", Can't find TrackSegments " << endl ;
558 return ;
559 }
560 TObjArray * emcrp = gime->EmcRecPoints() ;
561 if(!emcrp){
562 cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
563 << ", Can't find EmcRecPoints " << endl ;
564 return ;
565 }
566
567 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){
568 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
569
570 //find the closest primary
571 Int_t moduleNumberRec ;
572 Double_t recX, recZ ;
573 phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
574
575 Double_t minDistance = 100. ;
576 Int_t closestPrimary = -1 ;
577
578 //extract list of primaries: it is stored at EMC RecPoints
579 Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
580 Int_t numberofprimaries ;
581 Int_t * listofprimaries = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
582
583 Int_t index ;
584 const TParticle * primary ;
585 Double_t distance = minDistance ;
586 Double_t dX, dZ;
587 Double_t dXmin = 0.;
588 Double_t dZmin = 0. ;
589 for ( index = 0 ; index < numberofprimaries ; index++){
590 primary = gime->Primary(listofprimaries[index]) ;
591 Int_t moduleNumber ;
592 Double_t primX, primZ ;
593 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
594 if(moduleNumberRec == moduleNumber) {
595 dX = recX - primX;
596 dZ = recZ - primZ;
597 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
598 if(minDistance > distance) {
599 minDistance = distance ;
600 dXmin = dX;
601 dZmin = dZ;
602 closestPrimary = listofprimaries[index] ;
efad3788 603 }
604 }
9688c1dd 605 }
606
607 //if found primary, fill histograms
608 if(closestPrimary >=0 ){
609 const TParticle * primary = gime->Primary(closestPrimary) ;
610 if(primary->GetPdgCode() == 22){
611 hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
612 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
613 hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
614 hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
cce01dd5 615 }
9688c1dd 616 else
617 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
618 hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
efad3788 619 }
620 }
9688c1dd 621 }
efad3788 622 }
9688c1dd 623
efad3788 624 //write filled histograms
625 efile->cd() ;
626 hAllEnergy->Write(0,kOverwrite) ;
627 hPhotEnergy->Write(0,kOverwrite) ;
628 hEMEnergy->Write(0,kOverwrite) ;
629 // efile->Write() ;
630 efile->Close() ;
631 delete efile ;
9688c1dd 632
efad3788 633}
fc879520 634//____________________________________________________________________________
efad3788 635void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
fc879520 636{
efad3788 637 //fills two dimentional histo: energy vs. primary - reconstructed distance
638
639
640
641 TH2F * hAllPosition = 0; // Position of any RP with primary photon
642 TH2F * hPhotPosition = 0; // Position of kGAMMA with primary photon
643 TH2F * hEMPosition = 0; // Position of EM with primary photon
644
645 TH1F * hAllPositionX = 0; // X-Position Resolution of photons with photon primary
646 TH1F * hAllPositionZ = 0; // Z-Position Resolution of photons with photon primary
647
648
649 //opening file and reading histograms if any
650 TFile * pfile = new TFile("position.root","update");
651
652 hAllPosition = (TH2F*)pfile->Get("hAllPosition");
653 if(hAllPosition == 0)
654 hAllPosition = new TH2F("hAllPosition",
655 "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
656 hPhotPosition= (TH2F*)pfile->Get("hPhotPosition");
657 if(hPhotPosition == 0)
658 hPhotPosition = new TH2F("hPhotPosition",
659 "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
660 hEMPosition= (TH2F*)pfile->Get("hEMPosition") ;
661 if(hEMPosition == 0)
662 hEMPosition = new TH2F("hEMPosition",
663 "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.);
664 hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ;
665 if(hAllPositionX == 0)
666 hAllPositionX = new TH1F("hAllPositionX",
667 "Delta X of any RP with primary photon",100, -2., 2.);
668 hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ;
669 if(hAllPositionZ == 0)
670 hAllPositionZ = new TH1F("hAllPositionZ",
671 "Delta X of any RP with primary photon",100, -2., 2.);
672
673
7b7c1533 674 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
675 const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
a88eea6f 676
efad3788 677 Int_t ievent;
9688c1dd 678 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
7b7c1533 679 for ( ievent=0; ievent < maxevent ; ievent++){
9688c1dd 680
efad3788 681 //read the current event
7b7c1533 682 gime->Event(ievent) ;
9688c1dd 683 TClonesArray * rp = gime->RecParticles() ;
684 if(!rp) {
685 cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
686 << ", Can't find RecParticles " << endl ;
687 return ;
688 }
689 TClonesArray * ts = gime->TrackSegments() ;
690 if(!ts) {
691 cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
692 << ", Can't find TrackSegments " << endl ;
693 return ;
694 }
695 TObjArray * emcrp = gime->EmcRecPoints() ;
696 if(!emcrp){
697 cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
698 << ", Can't find EmcRecPoints " << endl ;
699 return ;
700 }
efad3788 701
9688c1dd 702
7b7c1533 703 const AliPHOSRecParticle * recParticle ;
efad3788 704 Int_t iRecParticle ;
9688c1dd 705 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
706 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
707
708 //find the closest primary
709 Int_t moduleNumberRec ;
710 Double_t recX, recZ ;
711 phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
712
713 Double_t minDistance = 100. ;
714 Int_t closestPrimary = -1 ;
715
716 //extract list of primaries: it is stored at EMC RecPoints
717 Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
718 Int_t numberofprimaries ;
719 Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
720
721 Int_t index ;
722 const TParticle * primary ;
723 Double_t distance = minDistance ;
724 Double_t dX = 1000; // incredible number
725 Double_t dZ = 1000; // for the case if no primary will be found
726 Double_t dXmin = 0.;
727 Double_t dZmin = 0. ;
728 for ( index = 0 ; index < numberofprimaries ; index++){
729 primary = gime->Primary(listofprimaries[index]) ;
730 Int_t moduleNumber ;
731 Double_t primX, primZ ;
732 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
733 if(moduleNumberRec == moduleNumber) {
734 dX = recX - primX;
735 dZ = recZ - primZ;
736 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
737 if(minDistance > distance) {
738 minDistance = distance ;
739 dXmin = dX;
740 dZmin = dZ;
741 closestPrimary = listofprimaries[index] ;
efad3788 742 }
743 }
9688c1dd 744 }
745
746 //if found primary, fill histograms
747 if(closestPrimary >=0 ){
748 const TParticle * primary = gime->Primary(closestPrimary) ;
749 if(primary->GetPdgCode() == 22){
750 hAllPosition->Fill(primary->Energy(), minDistance) ;
751 hAllPositionX->Fill(primary->Energy(), dX) ;
752 hAllPositionZ->Fill(primary->Energy(), dZ) ;
753 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
754 hPhotPosition->Fill(primary->Energy(), minDistance ) ;
755 hEMPosition->Fill(primary->Energy(), minDistance ) ;
cce01dd5 756 }
9688c1dd 757 else
758 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
759 hEMPosition->Fill(primary->Energy(), minDistance ) ;
efad3788 760 }
761 }
9688c1dd 762 }
efad3788 763 }
764
765 //Write output histgrams
766 pfile->cd() ;
767 hAllPosition->Write(0,kOverwrite) ;
768 hAllPositionX->Write(0,kOverwrite) ;
769 hAllPositionZ->Write(0,kOverwrite) ;
770 hPhotPosition->Write(0,kOverwrite) ;
771 hEMPosition->Write(0,kOverwrite) ;
772 pfile->Write() ;
773 pfile->Close() ;
774 delete pfile ;
9688c1dd 775
776
efad3788 777}
92862013 778//____________________________________________________________________________
efad3788 779void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
780// fills spectra of primary photons and several kinds of
781// reconstructed particles, so that analyzing them one can
782// estimate conatmination, efficiency of registration etc.
783
784 //define several general histograms
785 TH1F * hPrimary = 0; //spectrum (P_t distribution) of primary photons
786 TH1F * hAllRP = 0; //spectrum of all RecParticles in PHOS
787 TH1F * hPhot = 0; //spectrum of kGAMMA RecParticles
efad3788 788 TH1F * hShape = 0; //spectrum of all EM RecParticles
789 TH1F * hVeto = 0; //spectrum of all neutral RecParticles
790
791 //Now separate histograms in accoradance with primary
792 //primary - photon
793 TH1F * hPhotReg = 0; //Registeres as photon
794 TH1F * hPhotEM = 0; //Registered as EM
efad3788 795
796 //primary - n
797 TH1F * hNReg = 0; //Registeres as photon
798 TH1F * hNEM = 0; //Registered as EM
efad3788 799
800 //primary - nBar
801 TH1F * hNBarReg = 0; //Registeres as photon
802 TH1F * hNBarEM = 0; //Registered as EM
efad3788 803
804 //primary - charged hadron (pBar excluded)
805 TH1F * hChargedReg = 0; //Registeres as photon
806 TH1F * hChargedEM = 0; //Registered as EM
efad3788 807
808 //primary - pBar
809 TH1F * hPbarReg = 0; //Registeres as photon
810 TH1F * hPbarEM = 0; //Registered as EM
efad3788 811
812
813 //Reading histograms from the file
814 TFile * cfile = new TFile("contamination.root","update") ;
815
816 //read general histograms
817 hPrimary = (TH1F*) cfile->Get("hPrimary") ;
818 if(hPrimary == 0)
819 hPrimary= new TH1F("hPrimary", "Primary photon spectrum", 100, 0., 5.);
820 hAllRP = (TH1F*)cfile->Get("hAllRP") ;
821 if(hAllRP == 0)
822 hAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
823 hPhot = (TH1F*)cfile->Get("hPhot") ;
824 if(hPhot == 0)
825 hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.);
efad3788 826 hShape = (TH1F*) cfile->Get("hShape") ;
827 if(hShape == 0)
828 hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.);
829 hVeto= (TH1F*)cfile->Get("hVeto") ;
830 if(hVeto == 0)
831 hVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.);
832
833
834 //primary - photon
835 hPhotReg = (TH1F*)cfile->Get("hPhotReg");
836 if(hPhotReg == 0)
837 hPhotReg = new TH1F("hPhotReg","Photon registered as photon",100, 0., 5.);
838 hPhotEM =(TH1F*)cfile->Get("hPhotEM");
839 if(hPhotEM== 0)
840 hPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
efad3788 841
842 //primary - n
843 hNReg = (TH1F*)cfile->Get("hNReg");
844 if(hNReg== 0)
845 hNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.);
846 hNEM = (TH1F*)cfile->Get("hNEM");
847 if(hNEM== 0)
848 hNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
efad3788 849
850 //primary - nBar
851 hNBarReg =(TH1F*)cfile->Get("hNBarReg");
852 if(hNBarReg== 0)
853 hNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.);
854 hNBarEM =(TH1F*)cfile->Get("hNBarEM");
855 if(hNBarEM== 0)
856 hNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
efad3788 857
858 //primary - charged hadron (pBar excluded)
859 hChargedReg = (TH1F*)cfile->Get("hChargedReg");
860 if(hChargedReg== 0)
861 hChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
862 hChargedEM = (TH1F*)cfile->Get("hChargedEM");
863 if(hChargedEM== 0)
864 hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
9688c1dd 865
efad3788 866 //primary - pBar
867 hPbarReg = (TH1F*)cfile->Get("hPbarReg");
868 if(hPbarReg== 0)
869 hPbarReg= new TH1F("hPbarReg", "pBar registered as photon",100, 0., 5.);
870 hPbarEM = (TH1F*)cfile->Get("hPbarEM");
871 if(hPbarEM== 0)
872 hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.);
eecb6765 873
92862013 874
efad3788 875 //Now make some initializations
876
877 Int_t counter[8][5] ; //# of registered particles
878 Int_t i1,i2 ;
879 for(i1 = 0; i1<8; i1++)
880 for(i2 = 0; i2<5; i2++)
881 counter[i1][i2] = 0 ;
882
883
884
7b7c1533 885 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),RecPointsTitle) ;
886 const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
154fc86d 887
efad3788 888 Int_t ievent;
154fc86d 889 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
7b7c1533 890 for ( ievent=0; ievent < maxevent ; ievent++){
efad3788 891
7b7c1533 892 gime->Event(ievent) ;
efad3788 893
9688c1dd 894 TClonesArray * rp = gime->RecParticles() ;
895 if(!rp) {
896 cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
897 << ", Can't find RecParticles " << endl ;
898 return ;
899 }
900 TClonesArray * ts = gime->TrackSegments() ;
901 if(!ts) {
902 cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
903 << ", Can't find TrackSegments " << endl ;
904 return ;
905 }
906 TObjArray * emcrp = gime->EmcRecPoints() ;
907 if(!emcrp){
908 cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
909 << ", Can't find EmcRecPoints " << endl ;
910 return ;
911 }
912
913
efad3788 914 //=========== Make spectrum of the primary photons
7b7c1533 915 const TParticle * primary ;
efad3788 916 Int_t iPrimary ;
7b7c1533 917 for( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++){
918 primary = gime->Primary(iPrimary) ;
efad3788 919 Int_t primaryType = primary->GetPdgCode() ;
920 if( primaryType == 22 ) {
921 //check, if photons folls onto PHOS
922 Int_t moduleNumber ;
923 Double_t primX, primZ ;
7b7c1533 924 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
efad3788 925 if(moduleNumber)
926 hPrimary->Fill(primary->Energy()) ;
927
928 }
929
930 }
9688c1dd 931
efad3788 932 //========== Now scan over RecParticles
7b7c1533 933 const AliPHOSRecParticle * recParticle ;
efad3788 934 Int_t iRecParticle ;
9688c1dd 935 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
936 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
937 //fill histo spectrum of all RecParticles
938 hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
939
940 //==========find the closest primary
941 Int_t moduleNumberRec ;
942 Double_t recX, recZ ;
943 phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
944
945 Double_t minDistance = 100. ;
946 Int_t closestPrimary = -1 ;
947
948 //extract list of primaries: it is stored at EMC RecPoints
949 Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
950 Int_t numberofprimaries ;
951 Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
952 Int_t index ;
953 const TParticle * primary ;
954 Double_t distance = minDistance ;
955 Double_t dX, dZ;
956 Double_t dXmin = 0.;
957 Double_t dZmin = 0. ;
958 for ( index = 0 ; index < numberofprimaries ; index++){
959 primary = gime->Primary(listofprimaries[index]) ;
960 Int_t moduleNumber ;
961 Double_t primX, primZ ;
962 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
963 if(moduleNumberRec == moduleNumber) {
964 dX = recX - primX;
965 dZ = recZ - primZ;
966 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
967 if(minDistance > distance) {
968 minDistance = distance ;
969 dXmin = dX;
970 dZmin = dZ;
971 closestPrimary = listofprimaries[index] ;
efad3788 972 }
973 }
9688c1dd 974 }
975
976 //===========define the "type" of closest primary
977 if(closestPrimary >=0 ){
978 Int_t primaryCode = -1;
979 const TParticle * primary = gime->Primary(closestPrimary) ;
980 Int_t primaryType = primary->GetPdgCode() ;
981 if(primaryType == 22) // photon ?
982 primaryCode = 0 ;
983 else
984 if(primaryType == 2112) // neutron
985 primaryCode = 1 ;
efad3788 986 else
9688c1dd 987 if(primaryType == -2112) // Anti neutron
988 primaryCode = 2 ;
efad3788 989 else
9688c1dd 990 if(primaryType == -2122) //Anti proton
991 primaryCode = 4 ;
992 else {
993 TParticle tempo(*primary) ;
994 if(tempo.GetPDG()->Charge())
995 primaryCode = 3 ;
996 }
997
998 //==========Now look at the type of RecParticle
999 Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
1000 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
1001 hPhot->Fill(energy ) ;
1002 switch(primaryCode){
1003 case 0:
1004 hPhotReg->Fill(energy ) ;
1005 break ;
1006 case 1:
1007 hNReg->Fill(energy ) ;
1008 break ;
1009 case 2:
1010 hNBarReg->Fill(energy ) ;
1011 break ;
1012 case 3:
1013 hChargedReg->Fill(energy ) ;
1014 break ;
1015 case 4:
1016 hPbarReg->Fill(energy ) ;
1017 break ;
1018 default:
1019 break ;
efad3788 1020 }
9688c1dd 1021 }
1022 if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
1023 (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)||
1024 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)||
1025 (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower
1026 hShape->Fill(energy ) ;
1027 switch(primaryCode){
1028 case 0:
1029 hPhotEM->Fill(energy ) ;
1030 break ;
1031 case 1:
1032 hNEM->Fill(energy ) ;
1033 break ;
1034 case 2:
1035 hNBarEM->Fill(energy ) ;
1036 break ;
1037 case 3:
1038 hChargedEM->Fill(energy ) ;
1039 break ;
1040 case 4:
1041 hPbarEM->Fill(energy ) ;
1042 break ;
1043 default:
1044 break ;
efad3788 1045 }
1046 }
1047
9688c1dd 1048 if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
1049 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) ||
1050 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) ||
1051 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral
1052 hVeto->Fill(energy ) ;
1053
1054 //fill number of primaries identified as ...
1055 if(primaryCode >= 0) // Primary code defined
1056 counter[recParticle->GetType()][primaryCode]++ ;
1057
1058 }
1059
1060 } // no closest primary found
efad3788 1061 }
46b146ca 1062
46b146ca 1063
efad3788 1064 //=================== SaveHistograms
1065 cfile->cd() ;
1066 hPrimary->Write(0,kOverwrite);
1067 hAllRP->Write(0,kOverwrite);
1068 hPhot->Write(0,kOverwrite);
efad3788 1069 hShape->Write(0,kOverwrite);
1070 hVeto->Write(0,kOverwrite);
1071 hPhotReg->Write(0,kOverwrite);
1072 hPhotEM->Write(0,kOverwrite);
efad3788 1073 hNReg ->Write(0,kOverwrite);
1074 hNEM ->Write(0,kOverwrite);
efad3788 1075 hNBarReg ->Write(0,kOverwrite);
1076 hNBarEM ->Write(0,kOverwrite);
efad3788 1077 hChargedReg ->Write(0,kOverwrite);
1078 hChargedEM ->Write(0,kOverwrite);
efad3788 1079 hPbarReg ->Write(0,kOverwrite);
1080 hPbarEM ->Write(0,kOverwrite);
46b146ca 1081
efad3788 1082 cfile->Write(0,kOverwrite);
1083 cfile->Close();
1084 delete cfile ;
1085
69183710 1086
efad3788 1087 //print Final Table
154fc86d 1088 maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
7b7c1533 1089 // cout << "Resolutions: Analyzed " << gime->MaxEvent() << " event(s)" << endl ;
2aad621e 1090
7b7c1533 1091 cout << "Resolutions: Analyzed " << maxevent << " event(s)" << endl ;
efad3788 1092 cout << endl ;
b2a60966 1093
efad3788 1094 cout << " Primary: Photon Neutron Antineutron Charged hadron AntiProton" << endl ;
1095 cout << "--------------------------------------------------------------------------------" << endl ;
1096 cout << " kGAMMA: "
1097 << setw(8) << counter[2][0] << setw(9) << counter[2][1] << setw(13) << counter[2][2]
1098 << setw(15)<< counter[2][3] << setw(13) << counter[2][4] << endl ;
1099 cout << " kGAMMAHA: "
1100 << setw(8) << counter[3][0] << setw(9) << counter[3][1] << setw(13) << counter[3][2]
1101 << setw(15)<< counter[3][3] << setw(13) << counter[3][4] << endl ;
1102 cout << " kNEUTRALEM: "
1103 << setw(8) << counter[0][0] << setw(9) << counter[0][1] << setw(13) << counter[0][2]
1104 << setw(15)<< counter[0][3] << setw(13) << counter[0][4] << endl ;
1105 cout << " kNEUTRALHA: "
1106 << setw(8) << counter[1][0] << setw(9) << counter[1][1] << setw(13) << counter[1][2]
1107 << setw(15)<< counter[1][3] << setw(13) << counter[1][4] << endl ;
1108 cout << " kABSURDEM: "
1109 << setw(8) << counter[4][0] << setw(9) << counter[4][1] << setw(13) << counter[4][2]
1110 << setw(15)<< counter[4][3] << setw(13) << counter[4][4] << endl ;
1111 cout << " kABSURDHA: "
1112 << setw(8) << counter[5][0] << setw(9) << counter[5][1] << setw(13) << counter[5][2]
1113 << setw(15)<< counter[5][3] << setw(13) << counter[5][4] << endl ;
1114 cout << " kELECTRON: "
1115 << setw(8) << counter[6][0] << setw(9) << counter[6][1] << setw(13) << counter[6][2]
1116 << setw(15)<< counter[6][3] << setw(13) << counter[6][4] << endl ;
1117 cout << " kCHARGEDHA: "
1118 << setw(8) << counter[7][0] << setw(9) << counter[7][1] << setw(13) << counter[7][2]
1119 << setw(15)<< counter[7][3] << setw(13) << counter[7][4] << endl ;
1120 cout << "--------------------------------------------------------------------------------" << endl ;
1121
1122 Int_t totalInd = 0 ;
1123 for(i1 = 0; i1<8; i1++)
1124 for(i2 = 0; i2<5; i2++)
1125 totalInd+=counter[i1][i2] ;
1126 cout << "Indentified particles: " << totalInd << endl ;
46b146ca 1127
69183710 1128}
1129
46b146ca 1130
efad3788 1131