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