1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //_________________________________________________________________________
19 // Algorythm class to analyze PHOSv1 events:
20 // Construct histograms and displays them.
21 // Use the macro EditorBar.C for best access to the functionnalities
23 //*-- Author: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH)
24 //////////////////////////////////////////////////////////////////////////////
26 // --- ROOT system ---
32 #include "TParticle.h"
33 #include "TClonesArray.h"
39 // --- Standard library ---
44 // --- AliRoot header files ---
47 #include "AliPHOSAnalyze.h"
48 #include "AliPHOSClusterizerv1.h"
49 #include "AliPHOSTrackSegmentMakerv1.h"
50 #include "AliPHOSPIDv1.h"
51 #include "AliPHOSReconstructioner.h"
52 #include "AliPHOSDigit.h"
53 #include "AliPHOSTrackSegment.h"
54 #include "AliPHOSRecParticle.h"
55 #include "AliPHOSIndexToObject.h"
56 #include "AliPHOSHit.h"
57 #include "AliPHOSCPVHit.h"
58 #include "AliPHOSCpvRecPoint.h"
60 ClassImp(AliPHOSAnalyze)
62 //____________________________________________________________________________
63 AliPHOSAnalyze::AliPHOSAnalyze()
65 // default ctor (useless)
70 //____________________________________________________________________________
71 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
73 // ctor: analyze events from root file "name"
75 Bool_t ok = OpenRootFile(name) ;
77 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
80 //========== Get AliRun object from file
81 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
83 //=========== Get the PHOS object and associated geometry from the file
84 fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
85 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
87 //========== Initializes the Index to Object converter
88 fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
89 //========== Current event number
101 //____________________________________________________________________________
102 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
105 ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
108 //____________________________________________________________________________
109 void AliPHOSAnalyze::Copy(TObject & obj)
111 // copy an analysis into an other one
113 // I do nothing more because the copy is silly but the Code checkers requires one
116 //____________________________________________________________________________
117 AliPHOSAnalyze::~AliPHOSAnalyze()
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 ;}
130 //____________________________________________________________________________
131 void 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.
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.);
148 //========== Create the Clusterizer
149 fClu = new AliPHOSClusterizerv1() ;
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) ;
158 gAlice->GetEvent(Nevent);
160 TClonesArray * primaryList = gAlice->Particles();
162 TParticle * primary ;
164 for ( iPrimary = 0 ; iPrimary < primaryList->GetEntries() ; iPrimary++)
166 primary = (TParticle*)primaryList->At(iPrimary) ;
167 Int_t primaryType = primary->GetPdgCode() ;
168 if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
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()) ;
175 if( primaryType == 22 ) {
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()) ;
183 if( primaryType == -2112 ) {
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()) ;
193 fPHOS->SetTreeAddress() ;
195 gAlice->TreeD()->GetEvent(0) ;
196 gAlice->TreeR()->GetEvent(0) ;
198 TObjArray ** emcRecPoints = fPHOS->EmcRecPoints() ;
199 TObjArray ** ppsdRecPoints = fPHOS->PpsdRecPoints() ;
200 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
203 AliPHOSDigit * digit ;
205 for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
207 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
209 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
211 fGeom->RelPosInModule(relid,x,z) ;
212 Float_t e = fClu->Calibrate(digit->GetAmp()) ;
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) ;
219 ppsdLow->Fill(x,z,e) ;
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());
234 for(irecp = 0; irecp < (*ppsdRecPoints)->GetEntries() ; irecp ++){
235 AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)(*ppsdRecPoints)->At(irecp) ;
236 if(ppsd->GetPHOSMod()==Nmod){
237 ppsd->GetLocalPosition(pos) ;
239 ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
241 ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
245 AliPHOSRecParticle * recParticle ;
247 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
249 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
251 Int_t moduleNumberRec ;
252 Double_t recX, recZ ;
253 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
254 if(moduleNumberRec == Nmod){
256 Double_t minDistance = 5. ;
257 Int_t closestPrimary = -1 ;
259 Int_t numberofprimaries ;
260 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
262 TParticle * primary ;
263 Double_t distance = minDistance ;
265 for ( index = 0 ; index < numberofprimaries ; index++){
266 primary = (TParticle*)primaryList->At(listofprimaries[index]) ;
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)
274 minDistance = distance ;
275 closestPrimary = listofprimaries[index] ;
279 if(closestPrimary >=0 ){
281 Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
284 recPhot->Fill(recZ,recX,recParticle->Energy()) ;
286 if(primaryType==-2112)
287 recNbar->Fill(recZ,recX,recParticle->Energy()) ;
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") ;
306 //____________________________________________________________________________
307 void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )
310 // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
311 // for events from FirstEvent to Nevents
314 for ( ievent=firstEvent; ievent<nevents; ievent++) {
315 if (ievent==firstEvent) {
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) ;
322 if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
323 fClu->SetPpsdEnergyThreshold (0.0000002) ;
324 fClu->SetPpsdClusteringThreshold(0.0000001) ;
326 else if (strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
327 fClu->SetLocalMaxCutCPV(0.03) ;
328 fClu->SetLogWeightCutCPV(4.0) ;
329 fClu->SetCpvEnergyThreshold(0.09) ;
331 fClu->SetCalibrationParameters(0., 0.00000001) ;
333 //========== Creates the track segment maker
334 fTrs = new AliPHOSTrackSegmentMakerv1() ;
335 // fTrs->UnsetUnfoldFlag() ;
337 //========== Creates the particle identifier for GPS2 only
338 if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
339 fPID = new AliPHOSPIDv1() ;
340 fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
343 //========== Creates the Reconstructioner
344 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
345 if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
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 ;
352 //=========== Connects the various Tree's for evt
353 Int_t tracks = gAlice->GetEvent(ievent);
355 fPHOS->Hit2Digit(tracks) ;
357 //=========== Do the reconstruction
358 fPHOS->Reconstruction(fRec);
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 ;}
369 //-------------------------------------------------------------------------------------
370 void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast)
373 // Read and print generated and reconstructed hits in CPV
374 // for events from EvFirst to Nevent.
375 // If only EvFirst is defined, print only this one event.
376 // Author: Yuri Kharlov
380 if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
381 for ( Int_t ievent=EvFirst; ievent<=EvLast; ievent++) {
383 //========== Event Number>
384 cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
386 //=========== Connects the various Tree's for evt
387 Int_t ntracks = gAlice->GetEvent(ievent);
389 //========== Creating branches ===================================
390 AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
391 gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
393 AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->CpvRecPoints() ;
394 gAlice->TreeR()->SetBranchAddress( "PHOSCpvRP", cpvRecPoints ) ;
396 // Read and print CPV hits
398 AliPHOSCPVModule cpvModule;
399 TClonesArray *cpvHits;
401 AliPHOSCPVHit *cpvHit;
406 for (Int_t itrack=0; itrack<ntracks; itrack++) {
407 //=========== Get the Hits Tree for the Primary track itrack
409 gAlice->TreeH()->GetEvent(itrack);
411 for (iModule=0; iModule < fGeom->GetNCPVModules(); iModule++) {
412 cpvModule = fPHOS->GetCPVModule(iModule);
413 cpvHits = cpvModule.Hits();
414 nCPVhits = cpvHits->GetEntriesFast();
415 for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
417 cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
418 p = cpvHit->GetMomentum();
421 ipart = cpvHit->GetIpart();
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());
425 printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d\n",
431 // Read and print CPV reconstructed points
433 //=========== Gets the Reconstruction TTree
434 gAlice->TreeR()->GetEvent(0) ;
435 printf("Recpoints: %d\n",(*fPHOS->CpvRecPoints())->GetEntries());
436 TIter nextRP(*fPHOS->CpvRecPoints() ) ;
437 AliPHOSCpvRecPoint *cpvRecPoint ;
438 Int_t nRecPoints = 0;
439 while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
442 cpvRecPoint->GetLocalPosition(locpos);
443 Int_t phosModule = cpvRecPoint->GetPHOSMod();
444 printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cm\n",
445 phosModule,locpos.X(),locpos.Z());
447 printf("This event has %d generated hits and %d reconstructed points\n",
448 nGenHits,nRecPoints);
452 //____________________________________________________________________________
453 void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
456 // Analyzes CPV characteristics
457 // Author: Yuri Kharlov
463 TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.);
464 TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.);
465 TH1F *hDr = new TH1F("hDr" ,"CPV r-resolution@reconstruction",100, 0. , 5.);
466 TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5);
467 TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" , 21,-0.5,20.5);
468 TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" , 21,-0.5,20.5);
470 cout << "Start CPV Analysis"<< endl ;
471 for ( Int_t ievent=0; ievent<Nevents; ievent++) {
473 //========== Event Number>
474 // if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
475 cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
477 //=========== Connects the various Tree's for evt
478 Int_t ntracks = gAlice->GetEvent(ievent);
480 //========== Creating branches ===================================
481 AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
482 gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
484 AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
485 gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
487 // Create and fill arrays of hits for each CPV module
489 Int_t nOfModules = fGeom->GetNModules();
490 TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
492 for (iModule=0; iModule < nOfModules; iModule++)
493 hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
495 AliPHOSCPVModule cpvModule;
496 TClonesArray *cpvHits;
498 AliPHOSCPVHit *cpvHit;
503 // First go through all primary tracks and fill the arrays
504 // of hits per each CPV module
506 for (Int_t itrack=0; itrack<ntracks; itrack++) {
507 // Get the Hits Tree for the Primary track itrack
509 gAlice->TreeH()->GetEvent(itrack);
510 for (Int_t iModule=0; iModule < nOfModules; iModule++) {
511 cpvModule = fPHOS->GetCPVModule(iModule);
512 cpvHits = cpvModule.Hits();
513 nCPVhits = cpvHits->GetEntriesFast();
514 for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
515 cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
516 p = cpvHit->GetMomentum();
517 xzgen[0] = cpvHit->X();
518 xzgen[1] = cpvHit->Y();
519 ipart = cpvHit->GetIpart();
520 TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
521 new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
526 for (iModule=0; iModule < nOfModules; iModule++) {
527 Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
528 printf("Module %d has %d hits\n",iModule,nsum);
531 // Then go through reconstructed points and for each find
533 // The distance from the rec.point to the closest hit
534 // gives the coordinate resolution of the CPV
536 // Get the Reconstruction Tree
537 gAlice->TreeR()->GetEvent(0) ;
538 TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
539 AliPHOSCpvRecPoint *cpvRecPoint ;
541 while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
543 cpvRecPoint->GetLocalPosition(locpos);
544 Int_t phosModule = cpvRecPoint->GetPHOSMod();
545 Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity();
546 Int_t rpMultX, rpMultZ;
547 cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
548 Float_t xrec = locpos.X();
549 Float_t zrec = locpos.Z();
550 Float_t dxmin = 1.e+10;
551 Float_t dzmin = 1.e+10;
552 Float_t r2min = 1.e+10;
555 cpvHits = hitsPerModule[phosModule-1];
556 Int_t nCPVhits = cpvHits->GetEntriesFast();
557 for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
558 cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
561 r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
570 hDr ->Fill(TMath::Sqrt(r2min));
572 hNrpX->Fill(rpMultX);
573 hNrpZ->Fill(rpMultZ);
575 delete [] hitsPerModule;
579 Text_t outputname[80] ;
580 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
581 TFile output(outputname,"RECREATE");
593 TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
594 gStyle->SetOptStat(111111);
595 gStyle->SetOptFit(1);
596 gStyle->SetOptDate(1);
597 cpvCanvas->Divide(3,2);
600 gPad->SetFillColor(10);
601 hNrp->SetFillColor(16);
605 gPad->SetFillColor(10);
606 hNrpX->SetFillColor(16);
610 gPad->SetFillColor(10);
611 hNrpZ->SetFillColor(16);
615 gPad->SetFillColor(10);
616 hDx->SetFillColor(16);
621 gPad->SetFillColor(10);
622 hDz->SetFillColor(16);
627 gPad->SetFillColor(10);
628 hDr->SetFillColor(16);
631 cpvCanvas->Print("CPV.ps");
635 //____________________________________________________________________________
636 void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
638 // Calculates Real and Mixed invariant mass distributions
640 const Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
641 Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/nMixedEvents) ;
643 //========== Booking Histograms
644 TH2D * hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ;
645 TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
646 TH2D * hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ;
647 TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
650 Int_t eventInMixedLoop ;
652 Int_t nRecParticles[4];//nMixedEvents] ;
654 AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", nMixedEvents*1000) ;
656 for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++ ){
659 for ( ievent=0; ievent < nMixedEvents; ievent++){
661 Int_t absEventNumber = eventInMixedLoop*nMixedEvents + ievent ;
663 //=========== Connects the various Tree's for evt
664 gAlice->GetEvent(absEventNumber);
666 //========== Creating branches ===================================
667 fPHOS->SetTreeAddress() ;
669 gAlice->TreeD()->GetEvent(0) ;
670 gAlice->TreeR()->GetEvent(0) ;
672 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
675 AliPHOSRecParticle * recParticle ;
677 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
679 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
680 if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
681 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
682 new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
687 nRecParticles[ievent] = iRecPhot-1 ;
690 //Now calculate invariant mass:
692 Int_t nCurEvent = 0 ;
694 for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
695 AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
697 for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
698 AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
701 invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
702 (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
703 (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
704 (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
707 invMass = TMath::Sqrt(invMass);
710 pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
712 if(irp1 > nRecParticles[nCurEvent])
715 if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
716 hRealEM->Fill(invMass,pt);
717 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
718 hRealPhot->Fill(invMass,pt);
721 hMixedEM->Fill(invMass,pt);
722 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
723 hMixedPhot->Fill(invMass,pt);
726 } //loop over second rp
727 }//loop over first rp
728 allRecParticleList->Delete() ;
731 delete allRecParticleList ;
734 TFile output("invmass.root","RECREATE");
740 hMixedPhot->Write() ;
747 //____________________________________________________________________________
748 void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast)
751 // Read and print generated and reconstructed hits in EMC
752 // for events from EvFirst to Nevent.
753 // If only EvFirst is defined, print only this one event.
754 // Author: Yuri Kharlov
758 if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
760 for (ievent=EvFirst; ievent<=EvLast; ievent++) {
762 //========== Event Number>
763 cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
765 //=========== Connects the various Tree's for evt
766 Int_t ntracks = gAlice->GetEvent(ievent);
767 fPHOS->SetTreeAddress() ;
769 gAlice->TreeD()->GetEvent(0) ;
770 gAlice->TreeR()->GetEvent(0) ;
772 // Loop over reconstructed particles
774 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
775 AliPHOSRecParticle * recParticle ;
779 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
780 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
781 Float_t recE = recParticle->Energy();
782 primList = recParticle->GetPrimaries(nPrimary);
783 Int_t moduleNumberRec ;
784 Double_t recX, recZ ;
785 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
786 printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
787 moduleNumberRec,recX,recZ,recE,*primList);
790 // Read and print EMC hits from EMCn branches
792 AliPHOSCPVModule emcModule;
793 TClonesArray *emcHits;
795 AliPHOSCPVHit *emcHit;
798 Int_t ipart, primary;
800 for (Int_t itrack=0; itrack<ntracks; itrack++) {
801 //=========== Get the Hits Tree for the Primary track itrack
803 gAlice->TreeH()->GetEvent(itrack);
805 for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
806 emcModule = fPHOS->GetEMCModule(iModule);
807 emcHits = emcModule.Hits();
808 nEMChits = emcHits->GetEntriesFast();
809 for (Int_t ihit=0; ihit<nEMChits; ihit++) {
811 emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
812 p = emcHit->GetMomentum();
815 ipart = emcHit->GetIpart();
816 primary= emcHit->GetTrack();
817 printf("EMC hit A: module %d, ",iModule+1);
818 printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
819 p.Px(),p.Py(),p.Pz(),p.Energy());
820 printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
821 xgen,zgen,ipart,primary);
826 // // Read and print EMC hits from PHOS branch
828 // for (Int_t itrack=0; itrack<ntracks; itrack++) {
829 // //=========== Get the Hits Tree for the Primary track itrack
830 // gAlice->ResetHits();
831 // gAlice->TreeH()->GetEvent(itrack);
832 // TClonesArray *hits = fPHOS->Hits();
835 // for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
836 // hit = (AliPHOSHit*)hits->At(ihit) ;
837 // Float_t hitXYZ[3];
838 // hitXYZ[0] = hit->X();
839 // hitXYZ[1] = hit->Y();
840 // hitXYZ[2] = hit->Z();
841 // ipart = hit->GetPid();
842 // primary = hit->GetPrimary();
843 // Int_t absId = hit->GetId();
845 // fGeom->AbsToRelNumbering(absId, relId) ;
846 // Int_t module = relId[0];
847 // if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
848 // printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
849 // module,hitXYZ[0],hitXYZ[2],ipart,primary);
856 //____________________________________________________________________________
857 void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
860 // Read generated and reconstructed hits in EMC for Nevents events.
861 // Plots the coordinate and energy resolution histograms.
862 // Coordinate resolution is a difference between the reconstructed
863 // coordinate and the exact coordinate on the face of the PHOS
864 // Author: Yuri Kharlov
870 TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.);
871 TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.);
872 TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.);
874 TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
875 TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
876 TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
878 cout << "Start EMC Analysis"<< endl ;
879 for (Int_t ievent=0; ievent<Nevents; ievent++) {
881 //========== Event Number>
882 if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
883 cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
885 //=========== Connects the various Tree's for evt
886 Int_t ntracks = gAlice->GetEvent(ievent);
888 fPHOS->SetTreeAddress() ;
890 gAlice->TreeD()->GetEvent(0) ;
891 gAlice->TreeR()->GetEvent(0) ;
893 // Create and fill arrays of hits for each EMC module
895 Int_t nOfModules = fGeom->GetNModules();
896 TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
898 for (iModule=0; iModule < nOfModules; iModule++)
899 hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
901 AliPHOSCPVModule emcModule;
902 TClonesArray *emcHits;
904 AliPHOSCPVHit *emcHit;
906 // First go through all primary tracks and fill the arrays
907 // of hits per each EMC module
909 for (Int_t itrack=0; itrack<ntracks; itrack++) {
910 // Get the Hits Tree for the Primary track itrack
912 gAlice->TreeH()->GetEvent(itrack);
913 for (Int_t iModule=0; iModule < nOfModules; iModule++) {
914 emcModule = fPHOS->GetEMCModule(iModule);
915 emcHits = emcModule.Hits();
916 nEMChits = emcHits->GetEntriesFast();
917 for (Int_t ihit=0; ihit<nEMChits; ihit++) {
918 emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
919 TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
920 new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
926 // Loop over reconstructed particles
928 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
929 AliPHOSRecParticle * recParticle ;
930 Int_t nEMCrecs = (*recParticleList)->GetEntries();
932 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
933 Float_t recE = recParticle->Energy();
935 Double_t recX, recZ ;
936 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
938 // for this rec.point take the hit list in the same PHOS module
940 emcHits = hitsPerModule[phosModule-1];
941 Int_t nEMChits = emcHits->GetEntriesFast();
943 Float_t genX, genZ, genE;
944 for (Int_t ihit=0; ihit<nEMChits; ihit++) {
945 emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
948 genE = emcHit->GetMomentum().E();
950 Float_t dx = recX - genX;
951 Float_t dz = recZ - genZ;
952 Float_t de = recE - genE;
956 hDx2 ->Fill(genE,dx);
957 hDz2 ->Fill(genE,dz);
958 hDE2 ->Fill(genE,recE);
961 delete [] hitsPerModule;
965 Text_t outputname[80] ;
966 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
967 TFile output(outputname,"RECREATE");
979 TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
980 gStyle->SetOptStat(111111);
981 gStyle->SetOptFit(1);
982 gStyle->SetOptDate(1);
983 emcCanvas->Divide(3,1);
986 gPad->SetFillColor(10);
987 hDx1->SetFillColor(16);
991 gPad->SetFillColor(10);
992 hDz1->SetFillColor(16);
996 gPad->SetFillColor(10);
997 hDE1->SetFillColor(16);
1000 emcCanvas->Print("EMC.ps");
1004 //____________________________________________________________________________
1005 void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
1007 // analyzes Nevents events and calculate Energy and Position resolution as well as
1008 // probaility of correct indentifiing of the incident particle
1010 //========== Booking Histograms
1011 cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
1012 BookResolutionHistograms();
1014 Int_t counter[9][5] ;
1015 Int_t i1,i2,totalInd = 0 ;
1016 for(i1 = 0; i1<9; i1++)
1017 for(i2 = 0; i2<5; i2++)
1018 counter[i1][i2] = 0 ;
1020 Int_t totalPrimary = 0 ;
1021 Int_t totalRecPart = 0 ;
1022 Int_t totalRPwithPrim = 0 ;
1025 cout << "Start Analysing"<< endl ;
1026 for ( ievent=0; ievent<Nevents; ievent++)
1029 //========== Event Number>
1030 // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
1031 cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
1033 //=========== Connects the various Tree's for evt
1034 gAlice->GetEvent(ievent);
1036 //=========== Gets the Kine TTree
1037 gAlice->TreeK()->GetEvent(0) ;
1039 //=========== Gets the list of Primari Particles
1040 TClonesArray * primaryList = gAlice->Particles();
1042 TParticle * primary ;
1044 for ( iPrimary = 0 ; iPrimary < primaryList->GetEntries() ; iPrimary++)
1046 primary = (TParticle*)primaryList->UncheckedAt(iPrimary) ;
1047 Int_t primaryType = primary->GetPdgCode() ;
1048 if( primaryType == 22 ) {
1049 Int_t moduleNumber ;
1050 Double_t primX, primZ ;
1051 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1053 fhPrimary->Fill(primary->Energy()) ;
1054 if(primary->Energy() > 0.3)
1060 fPHOS->SetTreeAddress() ;
1062 gAlice->TreeD()->GetEvent(0) ;
1063 gAlice->TreeR()->GetEvent(0) ;
1065 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
1067 AliPHOSRecParticle * recParticle ;
1068 Int_t iRecParticle ;
1069 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
1071 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
1072 fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
1074 Int_t moduleNumberRec ;
1075 Double_t recX, recZ ;
1076 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
1078 Double_t minDistance = 100. ;
1079 Int_t closestPrimary = -1 ;
1081 Int_t numberofprimaries ;
1082 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
1084 TParticle * primary ;
1085 Double_t distance = minDistance ;
1087 Double_t dXmin = 0.;
1088 Double_t dZmin = 0. ;
1089 for ( index = 0 ; index < numberofprimaries ; index++){
1090 primary = (TParticle*)primaryList->UncheckedAt(listofprimaries[index]) ;
1091 Int_t moduleNumber ;
1092 Double_t primX, primZ ;
1093 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1094 if(moduleNumberRec == moduleNumber) {
1097 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
1098 if(minDistance > distance) {
1099 minDistance = distance ;
1102 closestPrimary = listofprimaries[index] ;
1108 if(closestPrimary >=0 ){
1111 Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
1112 // TParticlePDG* pDGparticle = ((TParticle *)primaryList->At(closestPrimary))->GetPDG();
1113 // Double_t charge = PDGparticle->Charge() ;
1115 // cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
1120 primaryCode = 0; //Photon
1121 fhAllEnergy ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy()) ;
1122 fhAllPosition ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), minDistance) ;
1123 fhAllPositionX->Fill(dXmin);
1124 fhAllPositionZ->Fill(dZmin);
1127 primaryCode = 1; //Electron
1130 primaryCode = 1; //positron
1133 primaryCode = 4; //K+
1136 primaryCode = 4; //K-
1139 primaryCode = 4; //K0s
1142 primaryCode = 4; //K0l
1145 primaryCode = 2; //K0l
1148 primaryCode = 2; //K0l
1151 primaryCode = 2; //K0l
1154 primaryCode = 2; //K0l
1157 primaryCode = 3; //ELSE
1161 switch(recParticle->GetType())
1163 case AliPHOSFastRecParticle::kGAMMA:
1164 if(primaryType == 22){
1165 fhPhotEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
1166 fhEMEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
1167 fhPPSDEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
1169 fhPhotPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1170 fhEMPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1171 fhPPSDPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1173 fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1174 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1175 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1177 fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1179 if(primaryType == 2112){ //neutron
1180 fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1181 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1182 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1185 if(primaryType == -2112){ //neutron ~
1186 fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1187 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1188 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1191 if(primaryCode == 2){
1192 fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1193 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1194 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1197 fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1198 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1199 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1200 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1201 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1202 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1203 counter[0][primaryCode]++;
1205 case AliPHOSFastRecParticle::kELECTRON:
1206 if(primaryType == 22){
1207 fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1208 fhEMEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
1209 fhEMPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1210 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1211 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1213 if(primaryType == 2112){ //neutron
1214 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1215 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1218 if(primaryType == -2112){ //neutron ~
1219 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1220 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1223 if(primaryCode == 2){
1224 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1225 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1228 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1229 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1230 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1231 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1232 counter[1][primaryCode]++;
1234 case AliPHOSFastRecParticle::kNEUTRALHA:
1235 if(primaryType == 22)
1236 fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1238 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1239 counter[2][primaryCode]++;
1241 case AliPHOSFastRecParticle::kNEUTRALEM:
1242 if(primaryType == 22){
1243 fhEMEnergy->Fill(((TParticle *)primaryList->At(closestPrimary))->Energy(),recParticle->Energy() ) ;
1244 fhEMPosition->Fill(((TParticle *)primaryList->At(closestPrimary))->Energy(),minDistance ) ;
1246 fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1247 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1249 if(primaryType == 2112) //neutron
1250 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1252 if(primaryType == -2112) //neutron ~
1253 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1255 if(primaryCode == 2)
1256 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1258 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1259 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1260 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1262 counter[3][primaryCode]++;
1264 case AliPHOSFastRecParticle::kCHARGEDHA:
1265 if(primaryType == 22) //photon
1266 fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1268 counter[4][primaryCode]++ ;
1270 case AliPHOSFastRecParticle::kGAMMAHA:
1271 if(primaryType == 22){ //photon
1272 fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1273 fhPPSDEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
1274 fhPPSDPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
1275 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1277 if(primaryType == 2112){ //neutron
1278 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1281 if(primaryType == -2112){ //neutron ~
1282 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1284 if(primaryCode == 2){
1285 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1288 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1289 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1290 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1291 counter[5][primaryCode]++ ;
1293 case AliPHOSFastRecParticle::kABSURDEM:
1294 counter[6][primaryCode]++ ;
1295 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1297 case AliPHOSFastRecParticle::kABSURDHA:
1298 counter[7][primaryCode]++ ;
1301 counter[8][primaryCode]++ ;
1308 cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
1309 cout << "Resolutions: Total primary " << totalPrimary << endl ;
1310 cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ;
1311 cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ;
1312 cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ;
1313 cout << " Detected as photon " << counter[0][0] << " " << counter[0][1] << " " << counter[0][2] << " " <<counter[0][3] << " " << counter[0][4] << endl ;
1314 cout << " Detected as electron " << counter[1][0] << " " << counter[1][1] << " " << counter[1][2] << " " <<counter[1][3] << " " << counter[1][4] << endl ;
1315 cout << " Detected as neutral hadron " << counter[2][0] << " " << counter[2][1] << " " << counter[2][2] << " " <<counter[2][3] << " " << counter[2][4] << endl ;
1316 cout << " Detected as neutral EM " << counter[3][0] << " " << counter[3][1] << " " << counter[3][2] << " " <<counter[3][3] << " " << counter[3][4] << endl ;
1317 cout << " Detected as charged hadron " << counter[4][0] << " " << counter[4][1] << " " << counter[4][2] << " " <<counter[4][3] << " " << counter[4][4] << endl ;
1318 cout << " Detected as gamma-hadron " << counter[5][0] << " " << counter[5][1] << " " << counter[5][2] << " " <<counter[5][3] << " " << counter[5][4] << endl ;
1319 cout << " Detected as Absurd EM " << counter[6][0] << " " << counter[6][1] << " " << counter[6][2] << " " <<counter[6][3] << " " << counter[6][4] << endl ;
1320 cout << " Detected as absurd hadron " << counter[7][0] << " " << counter[7][1] << " " << counter[7][2] << " " <<counter[7][3] << " " << counter[7][4] << endl ;
1321 cout << " Detected as undefined " << counter[8][0] << " " << counter[8][1] << " " << counter[8][2] << " " <<counter[8][3] << " " << counter[8][4] << endl ;
1323 for(i1 = 0; i1<9; i1++)
1324 for(i2 = 0; i2<5; i2++)
1325 totalInd+=counter[i1][i2] ;
1326 cout << "Indentified particles " << totalInd << endl ;
1331 //____________________________________________________________________________
1332 void AliPHOSAnalyze::BookingHistograms()
1334 // Books the histograms where the results of the analysis are stored (to be changed)
1337 delete fhVetoDigit ;
1338 delete fhConvertorDigit ;
1339 delete fhEmcCluster ;
1340 delete fhVetoCluster ;
1341 delete fhConvertorCluster ;
1342 delete fhConvertorEmc ;
1344 fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
1345 fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
1346 fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
1347 fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.);
1348 fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
1349 fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
1350 fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
1353 //____________________________________________________________________________
1354 void AliPHOSAnalyze::BookResolutionHistograms()
1356 // Books the histograms where the results of the Resolution analysis are stored
1359 // delete fhAllEnergy ;
1361 // delete fhPhotEnergy ;
1363 // delete fhEMEnergy ;
1365 // delete fhPPSDEnergy ;
1368 fhAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1369 fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1370 fhEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1371 fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1373 // if(fhAllPosition)
1374 // delete fhAllPosition ;
1375 // if(fhPhotPosition)
1376 // delete fhPhotPosition ;
1378 // delete fhEMPosition ;
1379 // if(fhPPSDPosition)
1380 // delete fhPPSDPosition ;
1383 fhAllPosition = new TH2F("hAllPosition", "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1384 fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1385 fhEMPosition = new TH2F("hEMPosition", "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1386 fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1388 fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
1389 fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
1392 // delete fhAllReg ;
1394 // delete fhPhotReg ;
1398 // delete fhNBarReg ;
1400 // delete fhChargedReg ;
1402 fhAllReg = new TH1F("hAllReg", "All primaries registered as photon", 100, 0., 5.);
1403 fhPhotReg = new TH1F("hPhotReg", "Photon registered as photon", 100, 0., 5.);
1404 fhNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.);
1405 fhNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.);
1406 fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
1411 // delete fhPhotEM ;
1415 // delete fhNBarEM ;
1417 // delete fhChargedEM ;
1419 fhAllEM = new TH1F("hAllEM", "All primary registered as EM",100, 0., 5.);
1420 fhPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
1421 fhNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
1422 fhNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
1423 fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
1426 // delete fhAllPPSD ;
1428 // delete fhPhotPPSD ;
1432 // delete fhNBarPPSD ;
1433 // if(fhChargedPPSD)
1434 // delete fhChargedPPSD ;
1436 fhAllPPSD = new TH1F("hAllPPSD", "All primary registered as PPSD",100, 0., 5.);
1437 fhPhotPPSD = new TH1F("hPhotPPSD", "Photon registered as PPSD", 100, 0., 5.);
1438 fhNPPSD = new TH1F("hNPPSD", "N registered as PPSD", 100, 0., 5.);
1439 fhNBarPPSD = new TH1F("hNBarPPSD", "NBar registered as PPSD", 100, 0., 5.);
1440 fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
1443 // delete fhPrimary ;
1444 fhPrimary= new TH1F("hPrimary", "hPrimary", 100, 0., 5.);
1455 fhAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
1456 fhVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.);
1457 fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
1458 fhPPSD = new TH1F("hPPSD", "All PPSD photon particles", 100, 0., 5.);
1462 // delete fhPhotPhot ;
1464 // delete fhPhotElec ;
1466 // delete fhPhotNeuH ;
1468 // delete fhPhotNuEM ;
1470 // delete fhPhotChHa ;
1472 // delete fhPhotGaHa ;
1474 fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.); //Photon registered as photon
1475 fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.); //Photon registered as Electron
1476 fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.); //Photon registered as Neutral Hadron
1477 fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.); //Photon registered as Neutral EM
1478 fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.); //Photon registered as Charged Hadron
1479 fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.); //Photon registered as Gamma-Hadron
1482 //____________________________________________________________________________
1483 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1485 // Open the root file named "name"
1487 fRootFile = new TFile(name, "update") ;
1488 return fRootFile->IsOpen() ;
1491 //____________________________________________________________________________
1492 void AliPHOSAnalyze::SaveHistograms()
1494 // Saves the histograms in a root file named "name.analyzed"
1496 Text_t outputname[80] ;
1497 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1498 TFile output(outputname,"RECREATE");
1502 fhAllEnergy->Write() ;
1504 fhPhotEnergy->Write() ;
1506 fhEMEnergy->Write() ;
1508 fhPPSDEnergy->Write() ;
1510 fhAllPosition->Write() ;
1512 fhAllPositionX->Write() ;
1514 fhAllPositionZ->Write() ;
1516 fhPhotPosition->Write() ;
1518 fhEMPosition->Write() ;
1520 fhPPSDPosition->Write() ;
1524 fhPhotReg->Write() ;
1528 fhNBarReg->Write() ;
1530 fhChargedReg->Write() ;
1540 fhChargedEM->Write() ;
1542 fhAllPPSD->Write() ;
1544 fhPhotPPSD->Write() ;
1548 fhNBarPPSD->Write() ;
1550 fhChargedPPSD->Write() ;
1552 fhPrimary->Write() ;
1562 fhPhotPhot->Write() ;
1564 fhPhotElec->Write() ;
1566 fhPhotNeuH->Write() ;
1568 fhPhotNuEM->Write() ;
1570 fhPhotNuEM->Write() ;
1572 fhPhotChHa->Write() ;
1574 fhPhotGaHa->Write() ;
1575 if(fhEnergyCorrelations)
1576 fhEnergyCorrelations->Write() ;
1581 //____________________________________________________________________________
1582 Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
1584 return ERecPart/0.8783 ;
1587 //____________________________________________________________________________
1588 void AliPHOSAnalyze::ResetHistograms()
1590 fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2)
1592 fhEmcDigit = 0 ; // Histo of digit energies in the Emc
1593 fhVetoDigit = 0 ; // Histo of digit energies in the Veto
1594 fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor
1595 fhEmcCluster = 0 ; // Histo of Cluster energies in Emc
1596 fhVetoCluster = 0 ; // Histo of Cluster energies in Veto
1597 fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor
1598 fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies
1601 fhPhotEnergy = 0 ; // Total spectrum of detected photons
1602 fhEMEnergy = 0 ; // Spectrum of detected electrons with electron primary
1605 fhAllPositionX = 0 ;
1606 fhAllPositionZ = 0 ;
1607 fhPhotPosition = 0 ;
1609 fhPPSDPosition = 0 ;