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