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