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