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