]> 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
New .det file in a new format for AliITSv5
[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//_________________________________________________________________________
5f20d3fb 19// Algorythm class to analyze PHOSv1 events:
b2a60966 20// Construct histograms and displays them.
21// Use the macro EditorBar.C for best access to the functionnalities
a3dfe79c 22//*--
2aad621e 23//*-- Author: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH)
6ad0bfa0 24//////////////////////////////////////////////////////////////////////////////
25
26// --- ROOT system ---
27
92862013 28#include "TFile.h"
29#include "TH1.h"
30#include "TPad.h"
6ad0bfa0 31#include "TH2.h"
32#include "TParticle.h"
33#include "TClonesArray.h"
34#include "TTree.h"
35#include "TMath.h"
36#include "TCanvas.h"
2bed9e3e 37#include "TStyle.h"
6ad0bfa0 38
39// --- Standard library ---
40
de9ec31b 41#include <iostream.h>
42#include <stdio.h>
92862013 43
6ad0bfa0 44// --- AliRoot header files ---
45
46#include "AliRun.h"
47#include "AliPHOSAnalyze.h"
48#include "AliPHOSClusterizerv1.h"
49#include "AliPHOSTrackSegmentMakerv1.h"
26d4b141 50#include "AliPHOSPIDv1.h"
6ad0bfa0 51#include "AliPHOSReconstructioner.h"
52#include "AliPHOSDigit.h"
53#include "AliPHOSTrackSegment.h"
54#include "AliPHOSRecParticle.h"
83974468 55#include "AliPHOSIndexToObject.h"
ed4205d8 56#include "AliPHOSHit.h"
867ede0e 57#include "AliPHOSCpvRecPoint.h"
6ad0bfa0 58
59ClassImp(AliPHOSAnalyze)
60
6ad0bfa0 61//____________________________________________________________________________
62 AliPHOSAnalyze::AliPHOSAnalyze()
63{
b2a60966 64 // default ctor (useless)
6ad0bfa0 65
66 fRootFile = 0 ;
67}
68
69//____________________________________________________________________________
70AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
71{
83974468 72 // ctor: analyze events from root file "name"
6ad0bfa0 73
92862013 74 Bool_t ok = OpenRootFile(name) ;
75 if ( !ok ) {
6ad0bfa0 76 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
77 }
78 else {
eecb6765 79 //========== Get AliRun object from file
80 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
fc879520 81
eecb6765 82 //=========== Get the PHOS object and associated geometry from the file
788dcca4 83 fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
eecb6765 84 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
fc879520 85
eecb6765 86 //========== Initializes the Index to Object converter
87 fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
88 //========== Current event number
89 fEvt = -999 ;
fc879520 90
6ad0bfa0 91 }
2bed9e3e 92 fDebugLevel = 0;
69183710 93 fClu = 0 ;
94 fPID = 0 ;
95 fTrs = 0 ;
96 fRec = 0 ;
46b146ca 97 ResetHistograms() ;
6ad0bfa0 98}
99
88714635 100//____________________________________________________________________________
101AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
102{
103 // copy ctor
104 ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
105}
106
107//____________________________________________________________________________
191c1c41 108void AliPHOSAnalyze::Copy(TObject & obj)
88714635 109{
110 // copy an analysis into an other one
111 TObject::Copy(obj) ;
112 // I do nothing more because the copy is silly but the Code checkers requires one
113}
114
6ad0bfa0 115//____________________________________________________________________________
116AliPHOSAnalyze::~AliPHOSAnalyze()
117{
118 // dtor
119
a3dfe79c 120 if(fRootFile->IsOpen()) fRootFile->Close() ;
121 if(fRootFile) {delete fRootFile ; fRootFile=0 ;}
122 if(fPHOS) {delete fPHOS ; fPHOS =0 ;}
123 if(fClu) {delete fClu ; fClu =0 ;}
124 if(fPID) {delete fPID ; fPID =0 ;}
125 if(fRec) {delete fRec ; fRec =0 ;}
126 if(fTrs) {delete fTrs ; fTrs =0 ;}
6ad0bfa0 127
128}
6ad0bfa0 129//____________________________________________________________________________
ed4205d8 130void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
131 //Draws pimary particles and reconstructed
132 //digits, RecPoints, RecPartices etc
133 //for event Nevent in the module Nmod.
134
135 TH2F * digitOccupancy = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.);
ad8cfaf4 136 TH2F * sdigitOccupancy = new TH2F("sdigitOccupancy","EMC sdigits", 64,-71.,71.,64,-71.,71.);
ed4205d8 137 TH2F * emcOccupancy = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.);
138 TH2F * ppsdUp = new TH2F("ppsdUp","PPSD Up digits", 128,-71.,71.,128,-71.,71.) ;
139 TH2F * ppsdUpCl = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ;
140 TH2F * ppsdLow = new TH2F("ppsdLow","PPSD Low digits", 128,-71.,71.,128,-71.,71.) ;
141 TH2F * ppsdLowCl = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ;
142 TH2F * nbar = new TH2F("nbar","Primary nbar", 64,-71.,71.,64,-71.,71.);
143 TH2F * phot = new TH2F("phot","Primary Photon", 64,-71.,71.,64,-71.,71.);
144 TH2F * charg = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.);
145 TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.);
146 TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", 64,-71.,71.,64,-71.,71.);
147
46b146ca 148 //========== Create the Clusterizer
149 fClu = new AliPHOSClusterizerv1() ;
83974468 150
ed4205d8 151 gAlice->GetEvent(Nevent);
152
ed4205d8 153 TParticle * primary ;
154 Int_t iPrimary ;
2ab0c725 155 for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
ed4205d8 156 {
2ab0c725 157 primary = gAlice->Particle(iPrimary) ;
ed4205d8 158 Int_t primaryType = primary->GetPdgCode() ;
159 if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
160 Int_t moduleNumber ;
161 Double_t primX, primZ ;
162 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
163 if(moduleNumber==Nmod)
164 charg->Fill(primZ,primX,primary->Energy()) ;
165 }
166 if( primaryType == 22 ) {
167 Int_t moduleNumber ;
168 Double_t primX, primZ ;
169 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
170 if(moduleNumber==Nmod)
171 phot->Fill(primZ,primX,primary->Energy()) ;
172 }
173 else{
174 if( primaryType == -2112 ) {
175 Int_t moduleNumber ;
176 Double_t primX, primZ ;
177 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
178 if(moduleNumber==Nmod)
179 nbar->Fill(primZ,primX,primary->Energy()) ;
180 }
181 }
182 }
69183710 183
ad8cfaf4 184 // fPHOS->SetTreeAddress() ;
92862013 185
ad8cfaf4 186 gAlice->TreeS()->GetEvent(0) ;
187
188 Int_t iSDigit ;
189 AliPHOSDigit * sdigit ;
ed4205d8 190
ad8cfaf4 191 if(fPHOS->SDigits()){
192 for(iSDigit = 0; iSDigit < fPHOS->SDigits()->GetEntries(); iSDigit++)
193 {
194 sdigit = (AliPHOSDigit *) fPHOS->SDigits()->At(iSDigit) ;
195 Int_t relid[4];
196 fGeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
197 Float_t x,z ;
198 fGeom->RelPosInModule(relid,x,z) ;
199 Float_t e = fPHOS->Calibrate(sdigit->GetAmp()) ;
200 if(relid[0]==Nmod){
201 if(relid[1]==0) //EMC
202 sdigitOccupancy->Fill(x,z,e) ;
203 if((relid[1]>0)&&(relid[1]<17))
204 ppsdUp->Fill(x,z,e) ;
205 if(relid[1]>16)
206 ppsdLow->Fill(x,z,e) ;
207 }
208 }
209 }
210 else{
211 cout << "No SDigits read " << endl ;
212 }
ed4205d8 213
ad8cfaf4 214 gAlice->TreeD()->GetEvent(0) ;
ed4205d8 215
ad8cfaf4 216 if(fPHOS->Digits()){
217 Int_t iDigit ;
218 AliPHOSDigit * digit ;
219 for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
220 {
221 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
222 Int_t relid[4];
223 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
224 Float_t x,z ;
225 fGeom->RelPosInModule(relid,x,z) ;
226 Float_t e = fClu->Calibrate(digit->GetAmp()) ;
227 if(relid[0]==Nmod){
228 if(relid[1]==0) //EMC
229 digitOccupancy->Fill(x,z,e) ;
230 if((relid[1]>0)&&(relid[1]<17))
231 ppsdUp->Fill(x,z,e) ;
232 if(relid[1]>16)
233 ppsdLow->Fill(x,z,e) ;
234 }
ed4205d8 235 }
ad8cfaf4 236 }
237 else{
238 cout << "No Digits read " << endl ;
239 }
240
241 gAlice->TreeR()->GetEvent(0) ;
242
243 TObjArray * emcRecPoints = fPHOS->EmcRecPoints() ;
244 TObjArray * ppsdRecPoints = fPHOS->PpsdRecPoints() ;
245 TClonesArray * recParticleList = fPHOS->RecParticles() ;
246
ed4205d8 247
248 Int_t irecp ;
249 TVector3 pos ;
250
ad8cfaf4 251 if(emcRecPoints ){
252 for(irecp = 0; irecp < emcRecPoints->GetEntries() ; irecp ++){
253 AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)emcRecPoints->At(irecp) ;
254 if(emc->GetPHOSMod()==Nmod){
255 emc->GetLocalPosition(pos) ;
256 emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
257 }
ed4205d8 258 }
259 }
ad8cfaf4 260 else{
261 cout << "No EMC rec points read " << endl ;
262 }
263
264 if(ppsdRecPoints ){
265 for(irecp = 0; irecp < ppsdRecPoints->GetEntries() ; irecp ++){
266 AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)ppsdRecPoints->At(irecp) ;
267 if(ppsd->GetPHOSMod()==Nmod){
268 ppsd->GetLocalPosition(pos) ;
269 if(ppsd->GetUp())
270 ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
271 else
272 ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
273 }
ed4205d8 274 }
275 }
ad8cfaf4 276 else{
277 cout << "No PPSD/CPV rec points read " << endl ;
278 }
279
ed4205d8 280 AliPHOSRecParticle * recParticle ;
281 Int_t iRecParticle ;
ad8cfaf4 282 if(recParticleList ){
283 for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
284 {
285 recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
ed4205d8 286
ad8cfaf4 287 Int_t moduleNumberRec ;
288 Double_t recX, recZ ;
289 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
290 if(moduleNumberRec == Nmod){
11e586af 291
ad8cfaf4 292 Double_t minDistance = 5. ;
293 Int_t closestPrimary = -1 ;
ed4205d8 294
ad8cfaf4 295 Int_t numberofprimaries ;
296 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
297 Int_t index ;
298 TParticle * primary ;
299 Double_t distance = minDistance ;
ed4205d8 300
ad8cfaf4 301 for ( index = 0 ; index < numberofprimaries ; index++){
302 primary = gAlice->Particle(listofprimaries[index]) ;
303 Int_t moduleNumber ;
304 Double_t primX, primZ ;
305 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
306 if(moduleNumberRec == moduleNumber)
307 distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
308 if(minDistance > distance)
309 {
310 minDistance = distance ;
311 closestPrimary = listofprimaries[index] ;
312 }
313 }
314
315 if(closestPrimary >=0 ){
316
317 Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
318
319 if(primaryType==22)
320 recPhot->Fill(recZ,recX,recParticle->Energy()) ;
321 else
322 if(primaryType==-2112)
323 recNbar->Fill(recZ,recX,recParticle->Energy()) ;
324 }
325 }
ed4205d8 326 }
ad8cfaf4 327 }
328 else{
329 cout << "Not Rec Prticles read " << endl ;
330 }
ed4205d8 331
332 digitOccupancy->Draw("box") ;
ad8cfaf4 333 sdigitOccupancy->SetLineColor(5) ;
334 sdigitOccupancy->Draw("box") ;
ed4205d8 335 emcOccupancy->SetLineColor(2) ;
336 emcOccupancy->Draw("boxsame") ;
337 ppsdUp->SetLineColor(3) ;
338 ppsdUp->Draw("boxsame") ;
339 ppsdLow->SetLineColor(4) ;
340 ppsdLow->Draw("boxsame") ;
341 phot->SetLineColor(8) ;
342 phot->Draw("boxsame") ;
343 nbar->SetLineColor(6) ;
344 nbar->Draw("boxsame") ;
345
346}
134ce69a 347//____________________________________________________________________________
ed4205d8 348 void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )
c3b9b3f9 349{
134ce69a 350
ed4205d8 351 // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
a3dfe79c 352 // for events from FirstEvent to Nevents
2bed9e3e 353
354 Int_t ievent ;
ed4205d8 355 for ( ievent=firstEvent; ievent<nevents; ievent++) {
356 if (ievent==firstEvent) {
2bed9e3e 357 cout << "Analyze > Starting Reconstructing " << endl ;
358 //========== Create the Clusterizer
359 fClu = new AliPHOSClusterizerv1() ;
fc879520 360
2bed9e3e 361 //========== Creates the track segment maker
362 fTrs = new AliPHOSTrackSegmentMakerv1() ;
fad3e5b9 363 // fTrs->UnsetUnfoldFlag() ;
2bed9e3e 364
fad3e5b9 365 //========== Creates the particle identifier
366 fPID = new AliPHOSPIDv1() ;
367 fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
fc879520 368
2bed9e3e 369 //========== Creates the Reconstructioner
370 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
371 if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
372 }
fad3e5b9 373
2bed9e3e 374 if (fDebugLevel != 0 ||
375 (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
376 cout << "======= Analyze ======> Event " << ievent+1 << endl ;
377
037cc66d 378 fPHOS->Enable() ;
fad3e5b9 379
037cc66d 380 gAlice->Hits2Digits() ;
fad3e5b9 381
2bed9e3e 382 //=========== Do the reconstruction
ed4205d8 383 fPHOS->Reconstruction(fRec);
fad3e5b9 384
ed4205d8 385 }
386
2bed9e3e 387 if(fClu) {delete fClu ; fClu =0 ;}
388 if(fPID) {delete fPID ; fPID =0 ;}
389 if(fRec) {delete fRec ; fRec =0 ;}
390 if(fTrs) {delete fTrs ; fTrs =0 ;}
391
392}
393
394//-------------------------------------------------------------------------------------
8c0cd6e9 395void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast)
2bed9e3e 396{
fad3e5b9 397// //
398// // Read and print generated and reconstructed hits in CPV
399// // for events from EvFirst to Nevent.
400// // If only EvFirst is defined, print only this one event.
401// // Author: Yuri Kharlov
402// // 12 October 2000
403// //
404
405// if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
406// for ( Int_t ievent=EvFirst; ievent<=EvLast; ievent++) {
8c0cd6e9 407
fad3e5b9 408// //========== Event Number>
409// cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
2bed9e3e 410
fad3e5b9 411// //=========== Connects the various Tree's for evt
412// Int_t ntracks = gAlice->GetEvent(ievent);
2bed9e3e 413
fad3e5b9 414// //========== Creating branches ===================================
415// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
416// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
2bed9e3e 417
fad3e5b9 418// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
419// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
2bed9e3e 420
fad3e5b9 421// // Read and print CPV hits
a3dfe79c 422
fad3e5b9 423// AliPHOSCPVModule cpvModule;
424// TClonesArray *cpvHits;
425// Int_t nCPVhits;
426// AliPHOSCPVHit *cpvHit;
427// TLorentzVector p;
428// Float_t xgen, zgen;
429// Int_t ipart;
430// Int_t nGenHits = 0;
431// for (Int_t itrack=0; itrack<ntracks; itrack++) {
432// //=========== Get the Hits Tree for the Primary track itrack
433// gAlice->ResetHits();
434// gAlice->TreeH()->GetEvent(itrack);
435// Int_t iModule = 0 ;
436// for (iModule=0; iModule < fGeom->GetNCPVModules(); iModule++) {
437// cpvModule = fPHOS->GetCPVModule(iModule);
438// cpvHits = cpvModule.Hits();
439// nCPVhits = cpvHits->GetEntriesFast();
440// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
441// nGenHits++;
442// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
443// p = cpvHit->GetMomentum();
444// xgen = cpvHit->X();
445// zgen = cpvHit->Y();
446// ipart = cpvHit->GetIpart();
447// printf("CPV hit in module %d: ",iModule+1);
448// printf(" p = (%f, %f, %f, %f) GeV,\n",
449// p.Px(),p.Py(),p.Pz(),p.Energy());
450// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d\n",
451// xgen,zgen,ipart);
452// }
453// }
454// }
455
456// // Read and print CPV reconstructed points
457
458// //=========== Gets the Reconstruction TTree
459// gAlice->TreeR()->GetEvent(0) ;
460// printf("Recpoints: %d\n",(*fPHOS->CpvRecPoints())->GetEntries());
461// TIter nextRP(*fPHOS->CpvRecPoints() ) ;
462// AliPHOSCpvRecPoint *cpvRecPoint ;
463// Int_t nRecPoints = 0;
464// while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
465// nRecPoints++;
466// TVector3 locpos;
467// cpvRecPoint->GetLocalPosition(locpos);
468// Int_t phosModule = cpvRecPoint->GetPHOSMod();
469// printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cm\n",
470// phosModule,locpos.X(),locpos.Z());
471// }
472// printf("This event has %d generated hits and %d reconstructed points\n",
473// nGenHits,nRecPoints);
474// }
2bed9e3e 475}
134ce69a 476
2bed9e3e 477//____________________________________________________________________________
478void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
479{
037cc66d 480// //
481// // Analyzes CPV characteristics
482// // Author: Yuri Kharlov
483// // 9 October 2000
484// //
485
486// // Book histograms
487
488// TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.);
489// TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.);
490// TH1F *hDr = new TH1F("hDr" ,"CPV r-resolution@reconstruction",100, 0. , 5.);
491// TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5);
492// TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" , 21,-0.5,20.5);
493// TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" , 21,-0.5,20.5);
494
495// cout << "Start CPV Analysis"<< endl ;
496// for ( Int_t ievent=0; ievent<Nevents; ievent++) {
2bed9e3e 497
037cc66d 498// //========== Event Number>
499// // if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
500// cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
2bed9e3e 501
037cc66d 502// //=========== Connects the various Tree's for evt
503// Int_t ntracks = gAlice->GetEvent(ievent);
2bed9e3e 504
037cc66d 505// //========== Creating branches ===================================
506// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
507// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
2bed9e3e 508
037cc66d 509// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
510// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
8c0cd6e9 511
037cc66d 512// // Create and fill arrays of hits for each CPV module
8c0cd6e9 513
037cc66d 514// Int_t nOfModules = fGeom->GetNModules();
515// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
516// Int_t iModule = 0;
517// for (iModule=0; iModule < nOfModules; iModule++)
518// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
8c0cd6e9 519
037cc66d 520// AliPHOSCPVModule cpvModule;
521// TClonesArray *cpvHits;
522// Int_t nCPVhits;
523// AliPHOSCPVHit *cpvHit;
524// TLorentzVector p;
525// Float_t xzgen[2];
526// Int_t ipart;
134ce69a 527
037cc66d 528// // First go through all primary tracks and fill the arrays
529// // of hits per each CPV module
2bed9e3e 530
037cc66d 531// for (Int_t itrack=0; itrack<ntracks; itrack++) {
532// // Get the Hits Tree for the Primary track itrack
533// gAlice->ResetHits();
534// gAlice->TreeH()->GetEvent(itrack);
535// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
536// cpvModule = fPHOS->GetCPVModule(iModule);
537// cpvHits = cpvModule.Hits();
538// nCPVhits = cpvHits->GetEntriesFast();
539// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
540// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
541// p = cpvHit->GetMomentum();
542// xzgen[0] = cpvHit->X();
543// xzgen[1] = cpvHit->Y();
544// ipart = cpvHit->GetIpart();
545// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
546// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
547// }
548// cpvModule.Clear();
549// }
550// }
551// for (iModule=0; iModule < nOfModules; iModule++) {
552// Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
553// printf("Module %d has %d hits\n",iModule,nsum);
554// }
555
556// // Then go through reconstructed points and for each find
557// // the closeset hit
558// // The distance from the rec.point to the closest hit
559// // gives the coordinate resolution of the CPV
560
561// // Get the Reconstruction Tree
562// gAlice->TreeR()->GetEvent(0) ;
563// TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
564// AliPHOSCpvRecPoint *cpvRecPoint ;
565// Float_t xgen, zgen;
566// while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
567// TVector3 locpos;
568// cpvRecPoint->GetLocalPosition(locpos);
569// Int_t phosModule = cpvRecPoint->GetPHOSMod();
570// Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity();
571// Int_t rpMultX, rpMultZ;
572// cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
573// Float_t xrec = locpos.X();
574// Float_t zrec = locpos.Z();
575// Float_t dxmin = 1.e+10;
576// Float_t dzmin = 1.e+10;
577// Float_t r2min = 1.e+10;
578// Float_t r2;
579
580// cpvHits = hitsPerModule[phosModule-1];
581// Int_t nCPVhits = cpvHits->GetEntriesFast();
582// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
583// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
584// xgen = cpvHit->X();
585// zgen = cpvHit->Y();
586// r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
587// if ( r2 < r2min ) {
588// r2min = r2;
589// dxmin = xgen - xrec;
590// dzmin = zgen - zrec;
591// }
592// }
593// hDx ->Fill(dxmin);
594// hDz ->Fill(dzmin);
595// hDr ->Fill(TMath::Sqrt(r2min));
596// hNrp ->Fill(rpMult);
597// hNrpX->Fill(rpMultX);
598// hNrpZ->Fill(rpMultZ);
599// }
600// delete [] hitsPerModule;
601// }
602// // Save histograms
603
604// Text_t outputname[80] ;
605// sprintf(outputname,"%s.analyzed",fRootFile->GetName());
606// TFile output(outputname,"RECREATE");
607// output.cd();
608
609// hDx ->Write() ;
610// hDz ->Write() ;
611// hDr ->Write() ;
612// hNrp ->Write() ;
613// hNrpX->Write() ;
614// hNrpZ->Write() ;
615
616// // Plot histograms
617
618// TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
619// gStyle->SetOptStat(111111);
620// gStyle->SetOptFit(1);
621// gStyle->SetOptDate(1);
622// cpvCanvas->Divide(3,2);
623
624// cpvCanvas->cd(1);
625// gPad->SetFillColor(10);
626// hNrp->SetFillColor(16);
627// hNrp->Draw();
628
629// cpvCanvas->cd(2);
630// gPad->SetFillColor(10);
631// hNrpX->SetFillColor(16);
632// hNrpX->Draw();
633
634// cpvCanvas->cd(3);
635// gPad->SetFillColor(10);
636// hNrpZ->SetFillColor(16);
637// hNrpZ->Draw();
638
639// cpvCanvas->cd(4);
640// gPad->SetFillColor(10);
641// hDx->SetFillColor(16);
642// hDx->Fit("gaus");
643// hDx->Draw();
644
645// cpvCanvas->cd(5);
646// gPad->SetFillColor(10);
647// hDz->SetFillColor(16);
648// hDz->Fit("gaus");
649// hDz->Draw();
650
651// cpvCanvas->cd(6);
652// gPad->SetFillColor(10);
653// hDr->SetFillColor(16);
654// hDr->Draw();
655
656// cpvCanvas->Print("CPV.ps");
c3b9b3f9 657
c3b9b3f9 658}
2bed9e3e 659
c3b9b3f9 660//____________________________________________________________________________
69183710 661 void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
c3b9b3f9 662{
69183710 663 // Calculates Real and Mixed invariant mass distributions
ed4205d8 664
2f04ed65 665 const Int_t knMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
666 Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/knMixedEvents) ;
69183710 667
668 //========== Booking Histograms
669 TH2D * hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ;
670 TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
671 TH2D * hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ;
672 TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
673
674 Int_t ievent;
ed4205d8 675 Int_t eventInMixedLoop ;
69183710 676
2f04ed65 677 Int_t nRecParticles[4];//knMixedEvents] ;
69183710 678
2f04ed65 679 AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", knMixedEvents*1000) ;
69183710 680
ed4205d8 681 for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++ ){
69183710 682 Int_t iRecPhot = 0 ;
c3b9b3f9 683
2f04ed65 684 for ( ievent=0; ievent < knMixedEvents; ievent++){
69183710 685
2f04ed65 686 Int_t absEventNumber = eventInMixedLoop*knMixedEvents + ievent ;
69183710 687
688 //=========== Connects the various Tree's for evt
ed4205d8 689 gAlice->GetEvent(absEventNumber);
690
69183710 691 //========== Creating branches ===================================
ed4205d8 692 fPHOS->SetTreeAddress() ;
69183710 693
ed4205d8 694 gAlice->TreeD()->GetEvent(0) ;
69183710 695 gAlice->TreeR()->GetEvent(0) ;
696
ad8cfaf4 697 TClonesArray * recParticleList = fPHOS->RecParticles() ;
ed4205d8 698
699
700 AliPHOSRecParticle * recParticle ;
69183710 701 Int_t iRecParticle ;
ad8cfaf4 702 for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
69183710 703 {
ad8cfaf4 704 recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
ed4205d8 705 if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
706 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
707 new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
69183710 708 iRecPhot++;
709 }
710 }
711
ed4205d8 712 nRecParticles[ievent] = iRecPhot-1 ;
c3b9b3f9 713 }
69183710 714
69183710 715 //Now calculate invariant mass:
716 Int_t irp1,irp2 ;
ed4205d8 717 Int_t nCurEvent = 0 ;
c3b9b3f9 718
ed4205d8 719 for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
720 AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
c3b9b3f9 721
ed4205d8 722 for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
723 AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
69183710 724
ed4205d8 725 Double_t invMass ;
726 invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
69183710 727 (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
728 (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
729 (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
730
ed4205d8 731 if(invMass> 0)
732 invMass = TMath::Sqrt(invMass);
69183710 733
ed4205d8 734 Double_t pt ;
735 pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
69183710 736
ed4205d8 737 if(irp1 > nRecParticles[nCurEvent])
738 nCurEvent++;
69183710 739
ed4205d8 740 if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
741 hRealEM->Fill(invMass,pt);
69183710 742 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
ed4205d8 743 hRealPhot->Fill(invMass,pt);
69183710 744 }
745 else{
ed4205d8 746 hMixedEM->Fill(invMass,pt);
69183710 747 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
ed4205d8 748 hMixedPhot->Fill(invMass,pt);
69183710 749 } //real-mixed
750
751 } //loop over second rp
752 }//loop over first rp
ed4205d8 753 allRecParticleList->Delete() ;
69183710 754 } //Loop over events
755
ed4205d8 756 delete allRecParticleList ;
69183710 757
758 //writing output
759 TFile output("invmass.root","RECREATE");
760 output.cd();
761
762 hRealEM->Write() ;
763 hRealPhot->Write() ;
764 hMixedEM->Write() ;
765 hMixedPhot->Write() ;
766
767 output.Write();
768 output.Close();
c3b9b3f9 769
770}
771
ed4205d8 772//____________________________________________________________________________
333d1c21 773 void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast)
ed4205d8 774{
037cc66d 775// //
776// // Read and print generated and reconstructed hits in EMC
777// // for events from EvFirst to Nevent.
778// // If only EvFirst is defined, print only this one event.
779// // Author: Yuri Kharlov
780// // 24 November 2000
781// //
782
783// if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
784// Int_t ievent;
785// for (ievent=EvFirst; ievent<=EvLast; ievent++) {
ed4205d8 786
037cc66d 787// //========== Event Number>
788// cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
ed4205d8 789
037cc66d 790// //=========== Connects the various Tree's for evt
791// Int_t ntracks = gAlice->GetEvent(ievent);
792// fPHOS->SetTreeAddress() ;
ed4205d8 793
037cc66d 794// gAlice->TreeD()->GetEvent(0) ;
795// gAlice->TreeR()->GetEvent(0) ;
ed4205d8 796
037cc66d 797// // Loop over reconstructed particles
ed4205d8 798
037cc66d 799// TClonesArray ** recParticleList = fPHOS->RecParticles() ;
800// AliPHOSRecParticle * recParticle ;
801// Int_t iRecParticle ;
802// Int_t *primList;
803// Int_t nPrimary;
804// for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
805// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
806// Float_t recE = recParticle->Energy();
807// primList = recParticle->GetPrimaries(nPrimary);
808// Int_t moduleNumberRec ;
809// Double_t recX, recZ ;
810// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
811// printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
812// moduleNumberRec,recX,recZ,recE,*primList);
813// }
ed4205d8 814
037cc66d 815// // Read and print EMC hits from EMCn branches
ed4205d8 816
037cc66d 817// AliPHOSCPVModule emcModule;
818// TClonesArray *emcHits;
819// Int_t nEMChits;
820// AliPHOSCPVHit *emcHit;
821// TLorentzVector p;
822// Float_t xgen, zgen;
823// Int_t ipart, primary;
824// Int_t nGenHits = 0;
825// for (Int_t itrack=0; itrack<ntracks; itrack++) {
826// //=========== Get the Hits Tree for the Primary track itrack
827// gAlice->ResetHits();
828// gAlice->TreeH()->GetEvent(itrack);
829// Int_t iModule = 0 ;
830// for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
831// emcModule = fPHOS->GetEMCModule(iModule);
832// emcHits = emcModule.Hits();
833// nEMChits = emcHits->GetEntriesFast();
834// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
835// nGenHits++;
836// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
837// p = emcHit->GetMomentum();
838// xgen = emcHit->X();
839// zgen = emcHit->Y();
840// ipart = emcHit->GetIpart();
841// primary= emcHit->GetTrack();
842// printf("EMC hit A: module %d, ",iModule+1);
843// printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
844// p.Px(),p.Py(),p.Pz(),p.Energy());
845// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
846// xgen,zgen,ipart,primary);
847// }
848// }
849// }
ed4205d8 850
037cc66d 851// // // Read and print EMC hits from PHOS branch
852
853// // for (Int_t itrack=0; itrack<ntracks; itrack++) {
854// // //=========== Get the Hits Tree for the Primary track itrack
855// // gAlice->ResetHits();
856// // gAlice->TreeH()->GetEvent(itrack);
857// // TClonesArray *hits = fPHOS->Hits();
858// // AliPHOSHit *hit ;
859// // Int_t ihit;
860// // for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
861// // hit = (AliPHOSHit*)hits->At(ihit) ;
862// // Float_t hitXYZ[3];
863// // hitXYZ[0] = hit->X();
864// // hitXYZ[1] = hit->Y();
865// // hitXYZ[2] = hit->Z();
866// // ipart = hit->GetPid();
867// // primary = hit->GetPrimary();
868// // Int_t absId = hit->GetId();
869// // Int_t relId[4];
870// // fGeom->AbsToRelNumbering(absId, relId) ;
871// // Int_t module = relId[0];
872// // if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
873// // printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
874// // module,hitXYZ[0],hitXYZ[2],ipart,primary);
875// // }
876// // }
ed4205d8 877
037cc66d 878// }
ed4205d8 879}
880
881//____________________________________________________________________________
882 void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
883{
037cc66d 884// //
885// // Read generated and reconstructed hits in EMC for Nevents events.
886// // Plots the coordinate and energy resolution histograms.
887// // Coordinate resolution is a difference between the reconstructed
888// // coordinate and the exact coordinate on the face of the PHOS
889// // Author: Yuri Kharlov
890// // 27 November 2000
891// //
892
893// // Book histograms
894
895// TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.);
896// TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.);
897// TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.);
898
899// TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
900// TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
901// TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
902
903// cout << "Start EMC Analysis"<< endl ;
904// for (Int_t ievent=0; ievent<Nevents; ievent++) {
ed4205d8 905
037cc66d 906// //========== Event Number>
907// if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
908// cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
ed4205d8 909
037cc66d 910// //=========== Connects the various Tree's for evt
911// Int_t ntracks = gAlice->GetEvent(ievent);
ed4205d8 912
037cc66d 913// fPHOS->SetTreeAddress() ;
ed4205d8 914
037cc66d 915// gAlice->TreeD()->GetEvent(0) ;
916// gAlice->TreeR()->GetEvent(0) ;
ed4205d8 917
037cc66d 918// // Create and fill arrays of hits for each EMC module
ed4205d8 919
037cc66d 920// Int_t nOfModules = fGeom->GetNModules();
921// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
922// Int_t iModule;
923// for (iModule=0; iModule < nOfModules; iModule++)
924// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
ed4205d8 925
037cc66d 926// AliPHOSCPVModule emcModule;
927// TClonesArray *emcHits;
928// Int_t nEMChits;
929// AliPHOSCPVHit *emcHit;
ed4205d8 930
037cc66d 931// // First go through all primary tracks and fill the arrays
932// // of hits per each EMC module
933
934// for (Int_t itrack=0; itrack<ntracks; itrack++) {
935// // Get the Hits Tree for the Primary track itrack
936// gAlice->ResetHits();
937// gAlice->TreeH()->GetEvent(itrack);
938// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
939// emcModule = fPHOS->GetEMCModule(iModule);
940// emcHits = emcModule.Hits();
941// nEMChits = emcHits->GetEntriesFast();
942// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
943// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
944// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
945// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
946// }
947// emcModule.Clear();
948// }
949// }
ed4205d8 950
037cc66d 951// // Loop over reconstructed particles
952
953// TClonesArray ** recParticleList = fPHOS->RecParticles() ;
954// AliPHOSRecParticle * recParticle ;
955// Int_t nEMCrecs = (*recParticleList)->GetEntries();
956// if (nEMCrecs == 1) {
957// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
958// Float_t recE = recParticle->Energy();
959// Int_t phosModule;
960// Double_t recX, recZ ;
961// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
962
963// // for this rec.point take the hit list in the same PHOS module
964
965// emcHits = hitsPerModule[phosModule-1];
966// Int_t nEMChits = emcHits->GetEntriesFast();
967// if (nEMChits == 1) {
968// Float_t genX, genZ, genE;
969// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
970// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
971// genX = emcHit->X();
972// genZ = emcHit->Y();
973// genE = emcHit->GetMomentum().E();
974// }
975// Float_t dx = recX - genX;
976// Float_t dz = recZ - genZ;
977// Float_t de = recE - genE;
978// hDx1 ->Fill(dx);
979// hDz1 ->Fill(dz);
980// hDE1 ->Fill(de);
981// hDx2 ->Fill(genE,dx);
982// hDz2 ->Fill(genE,dz);
983// hDE2 ->Fill(genE,recE);
984// }
985// }
986// delete [] hitsPerModule;
987// }
988// // Save histograms
989
990// Text_t outputname[80] ;
991// sprintf(outputname,"%s.analyzed",fRootFile->GetName());
992// TFile output(outputname,"RECREATE");
993// output.cd();
994
995// hDx1 ->Write() ;
996// hDz1 ->Write() ;
997// hDE1 ->Write() ;
998// hDx2 ->Write() ;
999// hDz2 ->Write() ;
1000// hDE2 ->Write() ;
1001
1002// // Plot histograms
1003
1004// TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
1005// gStyle->SetOptStat(111111);
1006// gStyle->SetOptFit(1);
1007// gStyle->SetOptDate(1);
1008// emcCanvas->Divide(3,1);
1009
1010// emcCanvas->cd(1);
1011// gPad->SetFillColor(10);
1012// hDx1->SetFillColor(16);
1013// hDx1->Draw();
1014
1015// emcCanvas->cd(2);
1016// gPad->SetFillColor(10);
1017// hDz1->SetFillColor(16);
1018// hDz1->Draw();
1019
1020// emcCanvas->cd(3);
1021// gPad->SetFillColor(10);
1022// hDE1->SetFillColor(16);
1023// hDE1->Draw();
1024
1025// emcCanvas->Print("EMC.ps");
ed4205d8 1026
1027}
1028
fc879520 1029//____________________________________________________________________________
1030 void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
1031{
1032 // analyzes Nevents events and calculate Energy and Position resolution as well as
1033 // probaility of correct indentifiing of the incident particle
134ce69a 1034
fc879520 1035 //========== Booking Histograms
1036 cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
1037 BookResolutionHistograms();
134ce69a 1038
ed4205d8 1039 Int_t counter[9][5] ;
1040 Int_t i1,i2,totalInd = 0 ;
fc879520 1041 for(i1 = 0; i1<9; i1++)
1042 for(i2 = 0; i2<5; i2++)
ed4205d8 1043 counter[i1][i2] = 0 ;
fc879520 1044
ed4205d8 1045 Int_t totalPrimary = 0 ;
1046 Int_t totalRecPart = 0 ;
1047 Int_t totalRPwithPrim = 0 ;
fc879520 1048 Int_t ievent;
1049
1050 cout << "Start Analysing"<< endl ;
1051 for ( ievent=0; ievent<Nevents; ievent++)
1052 {
1053
1054 //========== Event Number>
69183710 1055 // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
fc879520 1056 cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
1057
1058 //=========== Connects the various Tree's for evt
1059 gAlice->GetEvent(ievent);
134ce69a 1060
69183710 1061 //=========== Gets the Kine TTree
fc879520 1062 gAlice->TreeK()->GetEvent(0) ;
1063
1064 //=========== Gets the list of Primari Particles
fc879520 1065
ed4205d8 1066 TParticle * primary ;
fc879520 1067 Int_t iPrimary ;
2ab0c725 1068 for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
69183710 1069 {
2ab0c725 1070 primary = gAlice->Particle(iPrimary) ;
ed4205d8 1071 Int_t primaryType = primary->GetPdgCode() ;
1072 if( primaryType == 22 ) {
1073 Int_t moduleNumber ;
1074 Double_t primX, primZ ;
1075 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1076 if(moduleNumber){
1077 fhPrimary->Fill(primary->Energy()) ;
1078 if(primary->Energy() > 0.3)
1079 totalPrimary++ ;
69183710 1080 }
1081 }
1082 }
fc879520 1083
ed4205d8 1084 fPHOS->SetTreeAddress() ;
69183710 1085
ed4205d8 1086 gAlice->TreeD()->GetEvent(0) ;
fc879520 1087 gAlice->TreeR()->GetEvent(0) ;
69183710 1088
ad8cfaf4 1089 TClonesArray * recParticleList = fPHOS->RecParticles() ;
ed4205d8 1090
1091 AliPHOSRecParticle * recParticle ;
fc879520 1092 Int_t iRecParticle ;
ad8cfaf4 1093 for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
fc879520 1094 {
ad8cfaf4 1095 recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
ed4205d8 1096 fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
fc879520 1097
ed4205d8 1098 Int_t moduleNumberRec ;
1099 Double_t recX, recZ ;
1100 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
fc879520 1101
ed4205d8 1102 Double_t minDistance = 100. ;
1103 Int_t closestPrimary = -1 ;
fc879520 1104
1105 Int_t numberofprimaries ;
ed4205d8 1106 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
fc879520 1107 Int_t index ;
ed4205d8 1108 TParticle * primary ;
1109 Double_t distance = minDistance ;
1110 Double_t dX, dZ;
b3445972 1111 Double_t dXmin = 0.;
1112 Double_t dZmin = 0. ;
69183710 1113 for ( index = 0 ; index < numberofprimaries ; index++){
2ab0c725 1114 primary = gAlice->Particle(listofprimaries[index]) ;
ed4205d8 1115 Int_t moduleNumber ;
1116 Double_t primX, primZ ;
1117 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1118 if(moduleNumberRec == moduleNumber) {
1119 dX = recX - primX;
1120 dZ = recZ - primZ;
1121 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
1122 if(minDistance > distance) {
1123 minDistance = distance ;
1124 dXmin = dX;
1125 dZmin = dZ;
1126 closestPrimary = listofprimaries[index] ;
69183710 1127 }
ed4205d8 1128 }
69183710 1129 }
ed4205d8 1130 totalRecPart++ ;
134ce69a 1131
ed4205d8 1132 if(closestPrimary >=0 ){
1133 totalRPwithPrim++;
69183710 1134
2ab0c725 1135 Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
1136// TParticlePDG* pDGparticle = gAlice->ParticleAt(closestPrimary)->GetPDG();
2bed9e3e 1137// Double_t charge = PDGparticle->Charge() ;
8f5ada7b 1138// if(charge)
ed4205d8 1139// cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
1140 Int_t primaryCode ;
1141 switch(primaryType)
69183710 1142 {
1143 case 22:
ed4205d8 1144 primaryCode = 0; //Photon
2ab0c725 1145 fhAllEnergy ->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy()) ;
1146 fhAllPosition ->Fill(gAlice->Particle(closestPrimary)->Energy(), minDistance) ;
ed4205d8 1147 fhAllPositionX->Fill(dXmin);
1148 fhAllPositionZ->Fill(dZmin);
69183710 1149 break;
1150 case 11 :
ed4205d8 1151 primaryCode = 1; //Electron
69183710 1152 break;
1153 case -11 :
ed4205d8 1154 primaryCode = 1; //positron
69183710 1155 break;
1156 case 321 :
ed4205d8 1157 primaryCode = 4; //K+
69183710 1158 break;
1159 case -321 :
ed4205d8 1160 primaryCode = 4; //K-
69183710 1161 break;
1162 case 310 :
ed4205d8 1163 primaryCode = 4; //K0s
69183710 1164 break;
1165 case 130 :
ed4205d8 1166 primaryCode = 4; //K0l
69183710 1167 break;
8f5ada7b 1168 case 211 :
ed4205d8 1169 primaryCode = 2; //K0l
8f5ada7b 1170 break;
1171 case -211 :
ed4205d8 1172 primaryCode = 2; //K0l
8f5ada7b 1173 break;
1174 case 2212 :
ed4205d8 1175 primaryCode = 2; //K0l
8f5ada7b 1176 break;
1177 case -2212 :
ed4205d8 1178 primaryCode = 2; //K0l
8f5ada7b 1179 break;
69183710 1180 default:
ed4205d8 1181 primaryCode = 3; //ELSE
69183710 1182 break;
1183 }
1184
ed4205d8 1185 switch(recParticle->GetType())
69183710 1186 {
1187 case AliPHOSFastRecParticle::kGAMMA:
ed4205d8 1188 if(primaryType == 22){
2ab0c725 1189 fhPhotEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1190 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1191 fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
69183710 1192
2ab0c725 1193 fhPhotPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1194 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1195 fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
69183710 1196
ed4205d8 1197 fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1198 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1199 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1200
ed4205d8 1201 fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1202 }
ed4205d8 1203 if(primaryType == 2112){ //neutron
1204 fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1205 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1206 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1207 }
1208
ed4205d8 1209 if(primaryType == -2112){ //neutron ~
1210 fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1211 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1212 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1213
1214 }
ed4205d8 1215 if(primaryCode == 2){
1216 fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1217 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1218 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1219 }
1220
ed4205d8 1221 fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1222 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1223 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1224 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1225 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1226 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1227 counter[0][primaryCode]++;
69183710 1228 break;
1229 case AliPHOSFastRecParticle::kELECTRON:
ed4205d8 1230 if(primaryType == 22){
1231 fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
2ab0c725 1232 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1233 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
ed4205d8 1234 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1235 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1236 }
ed4205d8 1237 if(primaryType == 2112){ //neutron
1238 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1239 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1240 }
1241
ed4205d8 1242 if(primaryType == -2112){ //neutron ~
1243 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1244 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1245
1246 }
ed4205d8 1247 if(primaryCode == 2){
1248 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1249 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1250 }
1251
ed4205d8 1252 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1253 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1254 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1255 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1256 counter[1][primaryCode]++;
69183710 1257 break;
1258 case AliPHOSFastRecParticle::kNEUTRALHA:
ed4205d8 1259 if(primaryType == 22)
1260 fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1261
ed4205d8 1262 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1263 counter[2][primaryCode]++;
69183710 1264 break ;
1265 case AliPHOSFastRecParticle::kNEUTRALEM:
ed4205d8 1266 if(primaryType == 22){
2ab0c725 1267 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(),recParticle->Energy() ) ;
1268 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance ) ;
69183710 1269
ed4205d8 1270 fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1271 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1272 }
ed4205d8 1273 if(primaryType == 2112) //neutron
1274 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1275
ed4205d8 1276 if(primaryType == -2112) //neutron ~
1277 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1278
ed4205d8 1279 if(primaryCode == 2)
1280 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1281
ed4205d8 1282 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1283 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1284 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1285
ed4205d8 1286 counter[3][primaryCode]++;
69183710 1287 break ;
1288 case AliPHOSFastRecParticle::kCHARGEDHA:
ed4205d8 1289 if(primaryType == 22) //photon
1290 fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1291
ed4205d8 1292 counter[4][primaryCode]++ ;
69183710 1293 break ;
1294 case AliPHOSFastRecParticle::kGAMMAHA:
ed4205d8 1295 if(primaryType == 22){ //photon
1296 fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
2ab0c725 1297 fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1298 fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
ed4205d8 1299 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fc879520 1300 }
ed4205d8 1301 if(primaryType == 2112){ //neutron
1302 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fc879520 1303 }
69183710 1304
ed4205d8 1305 if(primaryType == -2112){ //neutron ~
1306 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fc879520 1307 }
ed4205d8 1308 if(primaryCode == 2){
1309 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fc879520 1310 }
69183710 1311
ed4205d8 1312 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1313 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1314 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1315 counter[5][primaryCode]++ ;
fc879520 1316 break ;
69183710 1317 case AliPHOSFastRecParticle::kABSURDEM:
ed4205d8 1318 counter[6][primaryCode]++ ;
1319 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
69183710 1320 break;
1321 case AliPHOSFastRecParticle::kABSURDHA:
ed4205d8 1322 counter[7][primaryCode]++ ;
69183710 1323 break;
1324 default:
ed4205d8 1325 counter[8][primaryCode]++ ;
69183710 1326 break;
1327 }
1328 }
fc879520 1329 }
1330 } // endfor
46b146ca 1331 SaveHistograms();
fc879520 1332 cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
ed4205d8 1333 cout << "Resolutions: Total primary " << totalPrimary << endl ;
1334 cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ;
1335 cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ;
fc879520 1336 cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ;
ed4205d8 1337 cout << " Detected as photon " << counter[0][0] << " " << counter[0][1] << " " << counter[0][2] << " " <<counter[0][3] << " " << counter[0][4] << endl ;
1338 cout << " Detected as electron " << counter[1][0] << " " << counter[1][1] << " " << counter[1][2] << " " <<counter[1][3] << " " << counter[1][4] << endl ;
1339 cout << " Detected as neutral hadron " << counter[2][0] << " " << counter[2][1] << " " << counter[2][2] << " " <<counter[2][3] << " " << counter[2][4] << endl ;
1340 cout << " Detected as neutral EM " << counter[3][0] << " " << counter[3][1] << " " << counter[3][2] << " " <<counter[3][3] << " " << counter[3][4] << endl ;
1341 cout << " Detected as charged hadron " << counter[4][0] << " " << counter[4][1] << " " << counter[4][2] << " " <<counter[4][3] << " " << counter[4][4] << endl ;
1342 cout << " Detected as gamma-hadron " << counter[5][0] << " " << counter[5][1] << " " << counter[5][2] << " " <<counter[5][3] << " " << counter[5][4] << endl ;
1343 cout << " Detected as Absurd EM " << counter[6][0] << " " << counter[6][1] << " " << counter[6][2] << " " <<counter[6][3] << " " << counter[6][4] << endl ;
1344 cout << " Detected as absurd hadron " << counter[7][0] << " " << counter[7][1] << " " << counter[7][2] << " " <<counter[7][3] << " " << counter[7][4] << endl ;
1345 cout << " Detected as undefined " << counter[8][0] << " " << counter[8][1] << " " << counter[8][2] << " " <<counter[8][3] << " " << counter[8][4] << endl ;
fc879520 1346
1347 for(i1 = 0; i1<9; i1++)
1348 for(i2 = 0; i2<5; i2++)
ed4205d8 1349 totalInd+=counter[i1][i2] ;
1350 cout << "Indentified particles " << totalInd << endl ;
fc879520 1351
92862013 1352} // endfunction
1353
1354
1355//____________________________________________________________________________
1356void AliPHOSAnalyze::BookingHistograms()
1357{
b2a60966 1358 // Books the histograms where the results of the analysis are stored (to be changed)
1359
eecb6765 1360 delete fhEmcDigit ;
1361 delete fhVetoDigit ;
1362 delete fhConvertorDigit ;
1363 delete fhEmcCluster ;
1364 delete fhVetoCluster ;
1365 delete fhConvertorCluster ;
1366 delete fhConvertorEmc ;
1367
46b146ca 1368 fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
1369 fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
1370 fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
1371 fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.);
1372 fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
1373 fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
1374 fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
92862013 1375
134ce69a 1376}
1377//____________________________________________________________________________
1378void AliPHOSAnalyze::BookResolutionHistograms()
1379{
1380 // Books the histograms where the results of the Resolution analysis are stored
1381
69183710 1382// if(fhAllEnergy)
1383// delete fhAllEnergy ;
1384// if(fhPhotEnergy)
1385// delete fhPhotEnergy ;
1386// if(fhEMEnergy)
1387// delete fhEMEnergy ;
1388// if(fhPPSDEnergy)
1389// delete fhPPSDEnergy ;
1390
1391
1392 fhAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1393 fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1394 fhEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1395 fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1396
1397// if(fhAllPosition)
1398// delete fhAllPosition ;
1399// if(fhPhotPosition)
1400// delete fhPhotPosition ;
1401// if(fhEMPosition)
1402// delete fhEMPosition ;
1403// if(fhPPSDPosition)
1404// delete fhPPSDPosition ;
1405
1406
1407 fhAllPosition = new TH2F("hAllPosition", "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1408 fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1409 fhEMPosition = new TH2F("hEMPosition", "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1410 fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1411
ed4205d8 1412 fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
1413 fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
1414
69183710 1415// if(fhAllReg)
1416// delete fhAllReg ;
1417// if(fhPhotReg)
1418// delete fhPhotReg ;
1419// if(fhNReg)
1420// delete fhNReg ;
1421// if(fhNBarReg)
1422// delete fhNBarReg ;
1423// if(fhChargedReg)
1424// delete fhChargedReg ;
46b146ca 1425
69183710 1426 fhAllReg = new TH1F("hAllReg", "All primaries registered as photon", 100, 0., 5.);
1427 fhPhotReg = new TH1F("hPhotReg", "Photon registered as photon", 100, 0., 5.);
1428 fhNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.);
1429 fhNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.);
1430 fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
46b146ca 1431
69183710 1432// if(fhAllEM)
1433// delete fhAllEM ;
1434// if(fhPhotEM)
1435// delete fhPhotEM ;
1436// if(fhNEM)
1437// delete fhNEM ;
1438// if(fhNBarEM)
1439// delete fhNBarEM ;
1440// if(fhChargedEM)
1441// delete fhChargedEM ;
46b146ca 1442
69183710 1443 fhAllEM = new TH1F("hAllEM", "All primary registered as EM",100, 0., 5.);
1444 fhPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
1445 fhNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
1446 fhNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
1447 fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
1448
1449// if(fhAllPPSD)
1450// delete fhAllPPSD ;
1451// if(fhPhotPPSD)
1452// delete fhPhotPPSD ;
1453// if(fhNPPSD)
1454// delete fhNPPSD ;
1455// if(fhNBarPPSD)
1456// delete fhNBarPPSD ;
1457// if(fhChargedPPSD)
1458// delete fhChargedPPSD ;
46b146ca 1459
69183710 1460 fhAllPPSD = new TH1F("hAllPPSD", "All primary registered as PPSD",100, 0., 5.);
1461 fhPhotPPSD = new TH1F("hPhotPPSD", "Photon registered as PPSD", 100, 0., 5.);
1462 fhNPPSD = new TH1F("hNPPSD", "N registered as PPSD", 100, 0., 5.);
1463 fhNBarPPSD = new TH1F("hNBarPPSD", "NBar registered as PPSD", 100, 0., 5.);
1464 fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
1465
1466// if(fhPrimary)
1467// delete fhPrimary ;
1468 fhPrimary= new TH1F("hPrimary", "hPrimary", 100, 0., 5.);
1469
1470// if(fhAllRP)
1471// delete fhAllRP ;
1472// if(fhVeto)
1473// delete fhVeto ;
1474// if(fhShape)
1475// delete fhShape ;
1476// if(fhPPSD)
1477// delete fhPPSD ;
1478
1479 fhAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
1480 fhVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.);
1481 fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
1482 fhPPSD = new TH1F("hPPSD", "All PPSD photon particles", 100, 0., 5.);
1483
1484
1485// if(fhPhotPhot)
1486// delete fhPhotPhot ;
1487// if(fhPhotElec)
1488// delete fhPhotElec ;
1489// if(fhPhotNeuH)
1490// delete fhPhotNeuH ;
1491// if(fhPhotNuEM)
1492// delete fhPhotNuEM ;
1493// if(fhPhotChHa)
1494// delete fhPhotChHa ;
1495// if(fhPhotGaHa)
1496// delete fhPhotGaHa ;
1497
1498 fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.); //Photon registered as photon
1499 fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.); //Photon registered as Electron
1500 fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.); //Photon registered as Neutral Hadron
1501 fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.); //Photon registered as Neutral EM
1502 fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.); //Photon registered as Charged Hadron
1503 fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.); //Photon registered as Gamma-Hadron
92862013 1504}
2aad621e 1505
6ad0bfa0 1506//____________________________________________________________________________
1507Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1508{
b2a60966 1509 // Open the root file named "name"
1510
1511 fRootFile = new TFile(name, "update") ;
6ad0bfa0 1512 return fRootFile->IsOpen() ;
1513}
ed4205d8 1514
92862013 1515//____________________________________________________________________________
46b146ca 1516void AliPHOSAnalyze::SaveHistograms()
134ce69a 1517{
1518 // Saves the histograms in a root file named "name.analyzed"
1519
1520 Text_t outputname[80] ;
1521 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1522 TFile output(outputname,"RECREATE");
1523 output.cd();
46b146ca 1524
69183710 1525 if (fhAllEnergy)
1526 fhAllEnergy->Write() ;
1527 if (fhPhotEnergy)
1528 fhPhotEnergy->Write() ;
1529 if(fhEMEnergy)
1530 fhEMEnergy->Write() ;
1531 if(fhPPSDEnergy)
1532 fhPPSDEnergy->Write() ;
1533 if(fhAllPosition)
1534 fhAllPosition->Write() ;
ed4205d8 1535 if(fhAllPositionX)
1536 fhAllPositionX->Write() ;
1537 if(fhAllPositionZ)
1538 fhAllPositionZ->Write() ;
69183710 1539 if(fhPhotPosition)
1540 fhPhotPosition->Write() ;
1541 if(fhEMPosition)
1542 fhEMPosition->Write() ;
1543 if(fhPPSDPosition)
1544 fhPPSDPosition->Write() ;
fc879520 1545 if (fhAllReg)
1546 fhAllReg->Write() ;
69183710 1547 if (fhPhotReg)
1548 fhPhotReg->Write() ;
fc879520 1549 if(fhNReg)
1550 fhNReg->Write() ;
1551 if(fhNBarReg)
1552 fhNBarReg->Write() ;
1553 if(fhChargedReg)
1554 fhChargedReg->Write() ;
fc879520 1555 if (fhAllEM)
1556 fhAllEM->Write() ;
69183710 1557 if (fhPhotEM)
1558 fhPhotEM->Write() ;
fc879520 1559 if(fhNEM)
1560 fhNEM->Write() ;
1561 if(fhNBarEM)
1562 fhNBarEM->Write() ;
1563 if(fhChargedEM)
1564 fhChargedEM->Write() ;
69183710 1565 if (fhAllPPSD)
1566 fhAllPPSD->Write() ;
1567 if (fhPhotPPSD)
1568 fhPhotPPSD->Write() ;
1569 if(fhNPPSD)
1570 fhNPPSD->Write() ;
1571 if(fhNBarPPSD)
1572 fhNBarPPSD->Write() ;
1573 if(fhChargedPPSD)
1574 fhChargedPPSD->Write() ;
fc879520 1575 if(fhPrimary)
1576 fhPrimary->Write() ;
69183710 1577 if(fhAllRP)
1578 fhAllRP->Write() ;
1579 if(fhVeto)
1580 fhVeto->Write() ;
1581 if(fhShape)
1582 fhShape->Write() ;
1583 if(fhPPSD)
1584 fhPPSD->Write() ;
fc879520 1585 if(fhPhotPhot)
1586 fhPhotPhot->Write() ;
1587 if(fhPhotElec)
1588 fhPhotElec->Write() ;
1589 if(fhPhotNeuH)
1590 fhPhotNeuH->Write() ;
1591 if(fhPhotNuEM)
1592 fhPhotNuEM->Write() ;
1593 if(fhPhotNuEM)
1594 fhPhotNuEM->Write() ;
1595 if(fhPhotChHa)
1596 fhPhotChHa->Write() ;
1597 if(fhPhotGaHa)
1598 fhPhotGaHa->Write() ;
46b146ca 1599 if(fhEnergyCorrelations)
1600 fhEnergyCorrelations->Write() ;
1601
92862013 1602 output.Write();
1603 output.Close();
1604}
69183710 1605//____________________________________________________________________________
1606Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
1607{
1608 return ERecPart/0.8783 ;
1609}
1610
46b146ca 1611//____________________________________________________________________________
1612void AliPHOSAnalyze::ResetHistograms()
1613{
1614 fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2)
1615
1616 fhEmcDigit = 0 ; // Histo of digit energies in the Emc
1617 fhVetoDigit = 0 ; // Histo of digit energies in the Veto
1618 fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor
1619 fhEmcCluster = 0 ; // Histo of Cluster energies in Emc
1620 fhVetoCluster = 0 ; // Histo of Cluster energies in Veto
1621 fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor
1622 fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies
1623
69183710 1624 fhAllEnergy = 0 ;
1625 fhPhotEnergy = 0 ; // Total spectrum of detected photons
1626 fhEMEnergy = 0 ; // Spectrum of detected electrons with electron primary
1627 fhPPSDEnergy = 0 ;
1628 fhAllPosition = 0 ;
ed4205d8 1629 fhAllPositionX = 0 ;
1630 fhAllPositionZ = 0 ;
69183710 1631 fhPhotPosition = 0 ;
1632 fhEMPosition = 0 ;
1633 fhPPSDPosition = 0 ;
1634
1635 fhPhotReg = 0 ;
46b146ca 1636 fhAllReg = 0 ;
1637 fhNReg = 0 ;
1638 fhNBarReg = 0 ;
1639 fhChargedReg = 0 ;
69183710 1640 fhPhotEM = 0 ;
46b146ca 1641 fhAllEM = 0 ;
1642 fhNEM = 0 ;
1643 fhNBarEM = 0 ;
1644 fhChargedEM = 0 ;
69183710 1645 fhPhotPPSD = 0 ;
1646 fhAllPPSD = 0 ;
1647 fhNPPSD = 0 ;
1648 fhNBarPPSD = 0 ;
1649 fhChargedPPSD = 0 ;
1650
46b146ca 1651 fhPrimary = 0 ;
1652
1653 fhPhotPhot = 0 ;
1654 fhPhotElec = 0 ;
1655 fhPhotNeuH = 0 ;
1656 fhPhotNuEM = 0 ;
1657 fhPhotChHa = 0 ;
1658 fhPhotGaHa = 0 ;
1659
46b146ca 1660}