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