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 gAlice->GetEvent(Nevent);
153 TParticle * primary ;
155 for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
157 primary = gAlice->Particle(iPrimary) ;
158 Int_t primaryType = primary->GetPdgCode() ;
159 if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
161 Double_t primX, primZ ;
162 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
163 if(moduleNumber==Nmod)
164 charg->Fill(primZ,primX,primary->Energy()) ;
166 if( primaryType == 22 ) {
168 Double_t primX, primZ ;
169 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
170 if(moduleNumber==Nmod)
171 phot->Fill(primZ,primX,primary->Energy()) ;
174 if( primaryType == -2112 ) {
176 Double_t primX, primZ ;
177 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
178 if(moduleNumber==Nmod)
179 nbar->Fill(primZ,primX,primary->Energy()) ;
184 fPHOS->SetTreeAddress() ;
186 gAlice->TreeD()->GetEvent(0) ;
187 gAlice->TreeR()->GetEvent(0) ;
189 TObjArray ** emcRecPoints = fPHOS->EmcRecPoints() ;
190 TObjArray ** ppsdRecPoints = fPHOS->PpsdRecPoints() ;
191 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
194 AliPHOSDigit * digit ;
196 for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
198 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
200 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
202 fGeom->RelPosInModule(relid,x,z) ;
203 Float_t e = fClu->Calibrate(digit->GetAmp()) ;
205 if(relid[1]==0) //EMC
206 digitOccupancy->Fill(x,z,e) ;
207 if((relid[1]>0)&&(relid[1]<17))
208 ppsdUp->Fill(x,z,e) ;
210 ppsdLow->Fill(x,z,e) ;
217 for(irecp = 0; irecp < (*emcRecPoints)->GetEntries() ; irecp ++){
218 AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)(*emcRecPoints)->At(irecp) ;
219 if(emc->GetPHOSMod()==Nmod){
220 emc->GetLocalPosition(pos) ;
221 emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
225 for(irecp = 0; irecp < (*ppsdRecPoints)->GetEntries() ; irecp ++){
226 AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)(*ppsdRecPoints)->At(irecp) ;
227 if(ppsd->GetPHOSMod()==Nmod){
228 ppsd->GetLocalPosition(pos) ;
230 ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
232 ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
236 AliPHOSRecParticle * recParticle ;
238 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
240 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
242 Int_t moduleNumberRec ;
243 Double_t recX, recZ ;
244 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
245 if(moduleNumberRec == Nmod){
247 Double_t minDistance = 5. ;
248 Int_t closestPrimary = -1 ;
250 Int_t numberofprimaries ;
251 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
253 TParticle * primary ;
254 Double_t distance = minDistance ;
256 for ( index = 0 ; index < numberofprimaries ; index++){
257 primary = gAlice->Particle(listofprimaries[index]) ;
259 Double_t primX, primZ ;
260 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
261 if(moduleNumberRec == moduleNumber)
262 distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
263 if(minDistance > distance)
265 minDistance = distance ;
266 closestPrimary = listofprimaries[index] ;
270 if(closestPrimary >=0 ){
272 Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
275 recPhot->Fill(recZ,recX,recParticle->Energy()) ;
277 if(primaryType==-2112)
278 recNbar->Fill(recZ,recX,recParticle->Energy()) ;
284 digitOccupancy->Draw("box") ;
285 emcOccupancy->SetLineColor(2) ;
286 emcOccupancy->Draw("boxsame") ;
287 ppsdUp->SetLineColor(3) ;
288 ppsdUp->Draw("boxsame") ;
289 ppsdLow->SetLineColor(4) ;
290 ppsdLow->Draw("boxsame") ;
291 phot->SetLineColor(8) ;
292 phot->Draw("boxsame") ;
293 nbar->SetLineColor(6) ;
294 nbar->Draw("boxsame") ;
297 //____________________________________________________________________________
298 void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )
301 // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
302 // for events from FirstEvent to Nevents
305 for ( ievent=firstEvent; ievent<nevents; ievent++) {
306 if (ievent==firstEvent) {
307 cout << "Analyze > Starting Reconstructing " << endl ;
308 //========== Create the Clusterizer
309 fClu = new AliPHOSClusterizerv1() ;
311 //========== Creates the track segment maker
312 fTrs = new AliPHOSTrackSegmentMakerv1() ;
313 // fTrs->UnsetUnfoldFlag() ;
315 //========== Creates the particle identifier
316 fPID = new AliPHOSPIDv1() ;
317 fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
319 //========== Creates the Reconstructioner
320 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
321 if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
324 if (fDebugLevel != 0 ||
325 (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
326 cout << "======= Analyze ======> Event " << ievent+1 << endl ;
328 //=========== Connects the various Tree's for evt
329 Int_t tracks = gAlice->GetEvent(ievent);
331 fPHOS->Hit2Digit(tracks) ;
333 //=========== Do the reconstruction
334 fPHOS->Reconstruction(fRec);
338 if(fClu) {delete fClu ; fClu =0 ;}
339 if(fPID) {delete fPID ; fPID =0 ;}
340 if(fRec) {delete fRec ; fRec =0 ;}
341 if(fTrs) {delete fTrs ; fTrs =0 ;}
345 //-------------------------------------------------------------------------------------
346 void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast)
349 // // Read and print generated and reconstructed hits in CPV
350 // // for events from EvFirst to Nevent.
351 // // If only EvFirst is defined, print only this one event.
352 // // Author: Yuri Kharlov
353 // // 12 October 2000
356 // if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
357 // for ( Int_t ievent=EvFirst; ievent<=EvLast; ievent++) {
359 // //========== Event Number>
360 // cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
362 // //=========== Connects the various Tree's for evt
363 // Int_t ntracks = gAlice->GetEvent(ievent);
365 // //========== Creating branches ===================================
366 // AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
367 // gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
369 // AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
370 // gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
372 // // Read and print CPV hits
374 // AliPHOSCPVModule cpvModule;
375 // TClonesArray *cpvHits;
377 // AliPHOSCPVHit *cpvHit;
379 // Float_t xgen, zgen;
381 // Int_t nGenHits = 0;
382 // for (Int_t itrack=0; itrack<ntracks; itrack++) {
383 // //=========== Get the Hits Tree for the Primary track itrack
384 // gAlice->ResetHits();
385 // gAlice->TreeH()->GetEvent(itrack);
386 // Int_t iModule = 0 ;
387 // for (iModule=0; iModule < fGeom->GetNCPVModules(); iModule++) {
388 // cpvModule = fPHOS->GetCPVModule(iModule);
389 // cpvHits = cpvModule.Hits();
390 // nCPVhits = cpvHits->GetEntriesFast();
391 // for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
393 // cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
394 // p = cpvHit->GetMomentum();
395 // xgen = cpvHit->X();
396 // zgen = cpvHit->Y();
397 // ipart = cpvHit->GetIpart();
398 // printf("CPV hit in module %d: ",iModule+1);
399 // printf(" p = (%f, %f, %f, %f) GeV,\n",
400 // p.Px(),p.Py(),p.Pz(),p.Energy());
401 // printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d\n",
407 // // Read and print CPV reconstructed points
409 // //=========== Gets the Reconstruction TTree
410 // gAlice->TreeR()->GetEvent(0) ;
411 // printf("Recpoints: %d\n",(*fPHOS->CpvRecPoints())->GetEntries());
412 // TIter nextRP(*fPHOS->CpvRecPoints() ) ;
413 // AliPHOSCpvRecPoint *cpvRecPoint ;
414 // Int_t nRecPoints = 0;
415 // while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
418 // cpvRecPoint->GetLocalPosition(locpos);
419 // Int_t phosModule = cpvRecPoint->GetPHOSMod();
420 // printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cm\n",
421 // phosModule,locpos.X(),locpos.Z());
423 // printf("This event has %d generated hits and %d reconstructed points\n",
424 // nGenHits,nRecPoints);
428 //____________________________________________________________________________
429 void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
432 // Analyzes CPV characteristics
433 // Author: Yuri Kharlov
439 TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.);
440 TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.);
441 TH1F *hDr = new TH1F("hDr" ,"CPV r-resolution@reconstruction",100, 0. , 5.);
442 TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5);
443 TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" , 21,-0.5,20.5);
444 TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" , 21,-0.5,20.5);
446 cout << "Start CPV Analysis"<< endl ;
447 for ( Int_t ievent=0; ievent<Nevents; ievent++) {
449 //========== Event Number>
450 // if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
451 cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
453 //=========== Connects the various Tree's for evt
454 Int_t ntracks = gAlice->GetEvent(ievent);
456 //========== Creating branches ===================================
457 AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
458 gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
460 AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
461 gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
463 // Create and fill arrays of hits for each CPV module
465 Int_t nOfModules = fGeom->GetNModules();
466 TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
468 for (iModule=0; iModule < nOfModules; iModule++)
469 hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
471 AliPHOSCPVModule cpvModule;
472 TClonesArray *cpvHits;
474 AliPHOSCPVHit *cpvHit;
479 // First go through all primary tracks and fill the arrays
480 // of hits per each CPV module
482 for (Int_t itrack=0; itrack<ntracks; itrack++) {
483 // Get the Hits Tree for the Primary track itrack
485 gAlice->TreeH()->GetEvent(itrack);
486 for (Int_t iModule=0; iModule < nOfModules; iModule++) {
487 cpvModule = fPHOS->GetCPVModule(iModule);
488 cpvHits = cpvModule.Hits();
489 nCPVhits = cpvHits->GetEntriesFast();
490 for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
491 cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
492 p = cpvHit->GetMomentum();
493 xzgen[0] = cpvHit->X();
494 xzgen[1] = cpvHit->Y();
495 ipart = cpvHit->GetIpart();
496 TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
497 new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
502 for (iModule=0; iModule < nOfModules; iModule++) {
503 Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
504 printf("Module %d has %d hits\n",iModule,nsum);
507 // Then go through reconstructed points and for each find
509 // The distance from the rec.point to the closest hit
510 // gives the coordinate resolution of the CPV
512 // Get the Reconstruction Tree
513 gAlice->TreeR()->GetEvent(0) ;
514 TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
515 AliPHOSCpvRecPoint *cpvRecPoint ;
517 while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
519 cpvRecPoint->GetLocalPosition(locpos);
520 Int_t phosModule = cpvRecPoint->GetPHOSMod();
521 Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity();
522 Int_t rpMultX, rpMultZ;
523 cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
524 Float_t xrec = locpos.X();
525 Float_t zrec = locpos.Z();
526 Float_t dxmin = 1.e+10;
527 Float_t dzmin = 1.e+10;
528 Float_t r2min = 1.e+10;
531 cpvHits = hitsPerModule[phosModule-1];
532 Int_t nCPVhits = cpvHits->GetEntriesFast();
533 for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
534 cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
537 r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
546 hDr ->Fill(TMath::Sqrt(r2min));
548 hNrpX->Fill(rpMultX);
549 hNrpZ->Fill(rpMultZ);
551 delete [] hitsPerModule;
555 Text_t outputname[80] ;
556 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
557 TFile output(outputname,"RECREATE");
569 TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
570 gStyle->SetOptStat(111111);
571 gStyle->SetOptFit(1);
572 gStyle->SetOptDate(1);
573 cpvCanvas->Divide(3,2);
576 gPad->SetFillColor(10);
577 hNrp->SetFillColor(16);
581 gPad->SetFillColor(10);
582 hNrpX->SetFillColor(16);
586 gPad->SetFillColor(10);
587 hNrpZ->SetFillColor(16);
591 gPad->SetFillColor(10);
592 hDx->SetFillColor(16);
597 gPad->SetFillColor(10);
598 hDz->SetFillColor(16);
603 gPad->SetFillColor(10);
604 hDr->SetFillColor(16);
607 cpvCanvas->Print("CPV.ps");
611 //____________________________________________________________________________
612 void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
614 // Calculates Real and Mixed invariant mass distributions
616 const Int_t knMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
617 Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/knMixedEvents) ;
619 //========== Booking Histograms
620 TH2D * hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ;
621 TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
622 TH2D * hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ;
623 TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
626 Int_t eventInMixedLoop ;
628 Int_t nRecParticles[4];//knMixedEvents] ;
630 AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", knMixedEvents*1000) ;
632 for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++ ){
635 for ( ievent=0; ievent < knMixedEvents; ievent++){
637 Int_t absEventNumber = eventInMixedLoop*knMixedEvents + ievent ;
639 //=========== Connects the various Tree's for evt
640 gAlice->GetEvent(absEventNumber);
642 //========== Creating branches ===================================
643 fPHOS->SetTreeAddress() ;
645 gAlice->TreeD()->GetEvent(0) ;
646 gAlice->TreeR()->GetEvent(0) ;
648 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
651 AliPHOSRecParticle * recParticle ;
653 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
655 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
656 if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
657 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
658 new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
663 nRecParticles[ievent] = iRecPhot-1 ;
666 //Now calculate invariant mass:
668 Int_t nCurEvent = 0 ;
670 for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
671 AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
673 for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
674 AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
677 invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
678 (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
679 (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
680 (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
683 invMass = TMath::Sqrt(invMass);
686 pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
688 if(irp1 > nRecParticles[nCurEvent])
691 if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
692 hRealEM->Fill(invMass,pt);
693 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
694 hRealPhot->Fill(invMass,pt);
697 hMixedEM->Fill(invMass,pt);
698 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
699 hMixedPhot->Fill(invMass,pt);
702 } //loop over second rp
703 }//loop over first rp
704 allRecParticleList->Delete() ;
707 delete allRecParticleList ;
710 TFile output("invmass.root","RECREATE");
716 hMixedPhot->Write() ;
723 //____________________________________________________________________________
724 void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast)
727 // Read and print generated and reconstructed hits in EMC
728 // for events from EvFirst to Nevent.
729 // If only EvFirst is defined, print only this one event.
730 // Author: Yuri Kharlov
734 if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
736 for (ievent=EvFirst; ievent<=EvLast; ievent++) {
738 //========== Event Number>
739 cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
741 //=========== Connects the various Tree's for evt
742 Int_t ntracks = gAlice->GetEvent(ievent);
743 fPHOS->SetTreeAddress() ;
745 gAlice->TreeD()->GetEvent(0) ;
746 gAlice->TreeR()->GetEvent(0) ;
748 // Loop over reconstructed particles
750 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
751 AliPHOSRecParticle * recParticle ;
755 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
756 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
757 Float_t recE = recParticle->Energy();
758 primList = recParticle->GetPrimaries(nPrimary);
759 Int_t moduleNumberRec ;
760 Double_t recX, recZ ;
761 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
762 printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
763 moduleNumberRec,recX,recZ,recE,*primList);
766 // Read and print EMC hits from EMCn branches
768 AliPHOSCPVModule emcModule;
769 TClonesArray *emcHits;
771 AliPHOSCPVHit *emcHit;
774 Int_t ipart, primary;
776 for (Int_t itrack=0; itrack<ntracks; itrack++) {
777 //=========== Get the Hits Tree for the Primary track itrack
779 gAlice->TreeH()->GetEvent(itrack);
781 for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
782 emcModule = fPHOS->GetEMCModule(iModule);
783 emcHits = emcModule.Hits();
784 nEMChits = emcHits->GetEntriesFast();
785 for (Int_t ihit=0; ihit<nEMChits; ihit++) {
787 emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
788 p = emcHit->GetMomentum();
791 ipart = emcHit->GetIpart();
792 primary= emcHit->GetTrack();
793 printf("EMC hit A: module %d, ",iModule+1);
794 printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
795 p.Px(),p.Py(),p.Pz(),p.Energy());
796 printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
797 xgen,zgen,ipart,primary);
802 // // Read and print EMC hits from PHOS branch
804 // for (Int_t itrack=0; itrack<ntracks; itrack++) {
805 // //=========== Get the Hits Tree for the Primary track itrack
806 // gAlice->ResetHits();
807 // gAlice->TreeH()->GetEvent(itrack);
808 // TClonesArray *hits = fPHOS->Hits();
811 // for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
812 // hit = (AliPHOSHit*)hits->At(ihit) ;
813 // Float_t hitXYZ[3];
814 // hitXYZ[0] = hit->X();
815 // hitXYZ[1] = hit->Y();
816 // hitXYZ[2] = hit->Z();
817 // ipart = hit->GetPid();
818 // primary = hit->GetPrimary();
819 // Int_t absId = hit->GetId();
821 // fGeom->AbsToRelNumbering(absId, relId) ;
822 // Int_t module = relId[0];
823 // if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
824 // printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
825 // module,hitXYZ[0],hitXYZ[2],ipart,primary);
832 //____________________________________________________________________________
833 void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
836 // Read generated and reconstructed hits in EMC for Nevents events.
837 // Plots the coordinate and energy resolution histograms.
838 // Coordinate resolution is a difference between the reconstructed
839 // coordinate and the exact coordinate on the face of the PHOS
840 // Author: Yuri Kharlov
846 TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.);
847 TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.);
848 TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.);
850 TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
851 TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
852 TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
854 cout << "Start EMC Analysis"<< endl ;
855 for (Int_t ievent=0; ievent<Nevents; ievent++) {
857 //========== Event Number>
858 if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
859 cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
861 //=========== Connects the various Tree's for evt
862 Int_t ntracks = gAlice->GetEvent(ievent);
864 fPHOS->SetTreeAddress() ;
866 gAlice->TreeD()->GetEvent(0) ;
867 gAlice->TreeR()->GetEvent(0) ;
869 // Create and fill arrays of hits for each EMC module
871 Int_t nOfModules = fGeom->GetNModules();
872 TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
874 for (iModule=0; iModule < nOfModules; iModule++)
875 hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
877 AliPHOSCPVModule emcModule;
878 TClonesArray *emcHits;
880 AliPHOSCPVHit *emcHit;
882 // First go through all primary tracks and fill the arrays
883 // of hits per each EMC module
885 for (Int_t itrack=0; itrack<ntracks; itrack++) {
886 // Get the Hits Tree for the Primary track itrack
888 gAlice->TreeH()->GetEvent(itrack);
889 for (Int_t iModule=0; iModule < nOfModules; iModule++) {
890 emcModule = fPHOS->GetEMCModule(iModule);
891 emcHits = emcModule.Hits();
892 nEMChits = emcHits->GetEntriesFast();
893 for (Int_t ihit=0; ihit<nEMChits; ihit++) {
894 emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
895 TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
896 new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
902 // Loop over reconstructed particles
904 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
905 AliPHOSRecParticle * recParticle ;
906 Int_t nEMCrecs = (*recParticleList)->GetEntries();
908 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
909 Float_t recE = recParticle->Energy();
911 Double_t recX, recZ ;
912 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
914 // for this rec.point take the hit list in the same PHOS module
916 emcHits = hitsPerModule[phosModule-1];
917 Int_t nEMChits = emcHits->GetEntriesFast();
919 Float_t genX, genZ, genE;
920 for (Int_t ihit=0; ihit<nEMChits; ihit++) {
921 emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
924 genE = emcHit->GetMomentum().E();
926 Float_t dx = recX - genX;
927 Float_t dz = recZ - genZ;
928 Float_t de = recE - genE;
932 hDx2 ->Fill(genE,dx);
933 hDz2 ->Fill(genE,dz);
934 hDE2 ->Fill(genE,recE);
937 delete [] hitsPerModule;
941 Text_t outputname[80] ;
942 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
943 TFile output(outputname,"RECREATE");
955 TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
956 gStyle->SetOptStat(111111);
957 gStyle->SetOptFit(1);
958 gStyle->SetOptDate(1);
959 emcCanvas->Divide(3,1);
962 gPad->SetFillColor(10);
963 hDx1->SetFillColor(16);
967 gPad->SetFillColor(10);
968 hDz1->SetFillColor(16);
972 gPad->SetFillColor(10);
973 hDE1->SetFillColor(16);
976 emcCanvas->Print("EMC.ps");
980 //____________________________________________________________________________
981 void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
983 // analyzes Nevents events and calculate Energy and Position resolution as well as
984 // probaility of correct indentifiing of the incident particle
986 //========== Booking Histograms
987 cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
988 BookResolutionHistograms();
990 Int_t counter[9][5] ;
991 Int_t i1,i2,totalInd = 0 ;
992 for(i1 = 0; i1<9; i1++)
993 for(i2 = 0; i2<5; i2++)
994 counter[i1][i2] = 0 ;
996 Int_t totalPrimary = 0 ;
997 Int_t totalRecPart = 0 ;
998 Int_t totalRPwithPrim = 0 ;
1001 cout << "Start Analysing"<< endl ;
1002 for ( ievent=0; ievent<Nevents; ievent++)
1005 //========== Event Number>
1006 // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
1007 cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
1009 //=========== Connects the various Tree's for evt
1010 gAlice->GetEvent(ievent);
1012 //=========== Gets the Kine TTree
1013 gAlice->TreeK()->GetEvent(0) ;
1015 //=========== Gets the list of Primari Particles
1017 TParticle * primary ;
1019 for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
1021 primary = gAlice->Particle(iPrimary) ;
1022 Int_t primaryType = primary->GetPdgCode() ;
1023 if( primaryType == 22 ) {
1024 Int_t moduleNumber ;
1025 Double_t primX, primZ ;
1026 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1028 fhPrimary->Fill(primary->Energy()) ;
1029 if(primary->Energy() > 0.3)
1035 fPHOS->SetTreeAddress() ;
1037 gAlice->TreeD()->GetEvent(0) ;
1038 gAlice->TreeR()->GetEvent(0) ;
1040 TClonesArray ** recParticleList = fPHOS->RecParticles() ;
1042 AliPHOSRecParticle * recParticle ;
1043 Int_t iRecParticle ;
1044 for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
1046 recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
1047 fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
1049 Int_t moduleNumberRec ;
1050 Double_t recX, recZ ;
1051 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
1053 Double_t minDistance = 100. ;
1054 Int_t closestPrimary = -1 ;
1056 Int_t numberofprimaries ;
1057 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
1059 TParticle * primary ;
1060 Double_t distance = minDistance ;
1062 Double_t dXmin = 0.;
1063 Double_t dZmin = 0. ;
1064 for ( index = 0 ; index < numberofprimaries ; index++){
1065 primary = gAlice->Particle(listofprimaries[index]) ;
1066 Int_t moduleNumber ;
1067 Double_t primX, primZ ;
1068 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1069 if(moduleNumberRec == moduleNumber) {
1072 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
1073 if(minDistance > distance) {
1074 minDistance = distance ;
1077 closestPrimary = listofprimaries[index] ;
1083 if(closestPrimary >=0 ){
1086 Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
1087 // TParticlePDG* pDGparticle = gAlice->ParticleAt(closestPrimary)->GetPDG();
1088 // Double_t charge = PDGparticle->Charge() ;
1090 // cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
1095 primaryCode = 0; //Photon
1096 fhAllEnergy ->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy()) ;
1097 fhAllPosition ->Fill(gAlice->Particle(closestPrimary)->Energy(), minDistance) ;
1098 fhAllPositionX->Fill(dXmin);
1099 fhAllPositionZ->Fill(dZmin);
1102 primaryCode = 1; //Electron
1105 primaryCode = 1; //positron
1108 primaryCode = 4; //K+
1111 primaryCode = 4; //K-
1114 primaryCode = 4; //K0s
1117 primaryCode = 4; //K0l
1120 primaryCode = 2; //K0l
1123 primaryCode = 2; //K0l
1126 primaryCode = 2; //K0l
1129 primaryCode = 2; //K0l
1132 primaryCode = 3; //ELSE
1136 switch(recParticle->GetType())
1138 case AliPHOSFastRecParticle::kGAMMA:
1139 if(primaryType == 22){
1140 fhPhotEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1141 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1142 fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1144 fhPhotPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1145 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1146 fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1148 fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1149 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1150 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1152 fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1154 if(primaryType == 2112){ //neutron
1155 fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1156 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1157 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1160 if(primaryType == -2112){ //neutron ~
1161 fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1162 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1163 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1166 if(primaryCode == 2){
1167 fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1168 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1169 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1172 fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1173 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1174 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1175 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1176 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1177 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1178 counter[0][primaryCode]++;
1180 case AliPHOSFastRecParticle::kELECTRON:
1181 if(primaryType == 22){
1182 fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1183 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1184 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1185 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1186 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1188 if(primaryType == 2112){ //neutron
1189 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1190 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1193 if(primaryType == -2112){ //neutron ~
1194 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1195 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1198 if(primaryCode == 2){
1199 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1200 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1203 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1204 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1205 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1206 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1207 counter[1][primaryCode]++;
1209 case AliPHOSFastRecParticle::kNEUTRALHA:
1210 if(primaryType == 22)
1211 fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1213 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1214 counter[2][primaryCode]++;
1216 case AliPHOSFastRecParticle::kNEUTRALEM:
1217 if(primaryType == 22){
1218 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(),recParticle->Energy() ) ;
1219 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance ) ;
1221 fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1222 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1224 if(primaryType == 2112) //neutron
1225 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1227 if(primaryType == -2112) //neutron ~
1228 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1230 if(primaryCode == 2)
1231 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1233 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1234 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1235 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1237 counter[3][primaryCode]++;
1239 case AliPHOSFastRecParticle::kCHARGEDHA:
1240 if(primaryType == 22) //photon
1241 fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1243 counter[4][primaryCode]++ ;
1245 case AliPHOSFastRecParticle::kGAMMAHA:
1246 if(primaryType == 22){ //photon
1247 fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1248 fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1249 fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1250 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1252 if(primaryType == 2112){ //neutron
1253 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1256 if(primaryType == -2112){ //neutron ~
1257 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1259 if(primaryCode == 2){
1260 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1263 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1264 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1265 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1266 counter[5][primaryCode]++ ;
1268 case AliPHOSFastRecParticle::kABSURDEM:
1269 counter[6][primaryCode]++ ;
1270 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1272 case AliPHOSFastRecParticle::kABSURDHA:
1273 counter[7][primaryCode]++ ;
1276 counter[8][primaryCode]++ ;
1283 cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
1284 cout << "Resolutions: Total primary " << totalPrimary << endl ;
1285 cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ;
1286 cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ;
1287 cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ;
1288 cout << " Detected as photon " << counter[0][0] << " " << counter[0][1] << " " << counter[0][2] << " " <<counter[0][3] << " " << counter[0][4] << endl ;
1289 cout << " Detected as electron " << counter[1][0] << " " << counter[1][1] << " " << counter[1][2] << " " <<counter[1][3] << " " << counter[1][4] << endl ;
1290 cout << " Detected as neutral hadron " << counter[2][0] << " " << counter[2][1] << " " << counter[2][2] << " " <<counter[2][3] << " " << counter[2][4] << endl ;
1291 cout << " Detected as neutral EM " << counter[3][0] << " " << counter[3][1] << " " << counter[3][2] << " " <<counter[3][3] << " " << counter[3][4] << endl ;
1292 cout << " Detected as charged hadron " << counter[4][0] << " " << counter[4][1] << " " << counter[4][2] << " " <<counter[4][3] << " " << counter[4][4] << endl ;
1293 cout << " Detected as gamma-hadron " << counter[5][0] << " " << counter[5][1] << " " << counter[5][2] << " " <<counter[5][3] << " " << counter[5][4] << endl ;
1294 cout << " Detected as Absurd EM " << counter[6][0] << " " << counter[6][1] << " " << counter[6][2] << " " <<counter[6][3] << " " << counter[6][4] << endl ;
1295 cout << " Detected as absurd hadron " << counter[7][0] << " " << counter[7][1] << " " << counter[7][2] << " " <<counter[7][3] << " " << counter[7][4] << endl ;
1296 cout << " Detected as undefined " << counter[8][0] << " " << counter[8][1] << " " << counter[8][2] << " " <<counter[8][3] << " " << counter[8][4] << endl ;
1298 for(i1 = 0; i1<9; i1++)
1299 for(i2 = 0; i2<5; i2++)
1300 totalInd+=counter[i1][i2] ;
1301 cout << "Indentified particles " << totalInd << endl ;
1306 //____________________________________________________________________________
1307 void AliPHOSAnalyze::BookingHistograms()
1309 // Books the histograms where the results of the analysis are stored (to be changed)
1312 delete fhVetoDigit ;
1313 delete fhConvertorDigit ;
1314 delete fhEmcCluster ;
1315 delete fhVetoCluster ;
1316 delete fhConvertorCluster ;
1317 delete fhConvertorEmc ;
1319 fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
1320 fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
1321 fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
1322 fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.);
1323 fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
1324 fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
1325 fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
1328 //____________________________________________________________________________
1329 void AliPHOSAnalyze::BookResolutionHistograms()
1331 // Books the histograms where the results of the Resolution analysis are stored
1334 // delete fhAllEnergy ;
1336 // delete fhPhotEnergy ;
1338 // delete fhEMEnergy ;
1340 // delete fhPPSDEnergy ;
1343 fhAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1344 fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1345 fhEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1346 fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1348 // if(fhAllPosition)
1349 // delete fhAllPosition ;
1350 // if(fhPhotPosition)
1351 // delete fhPhotPosition ;
1353 // delete fhEMPosition ;
1354 // if(fhPPSDPosition)
1355 // delete fhPPSDPosition ;
1358 fhAllPosition = new TH2F("hAllPosition", "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1359 fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1360 fhEMPosition = new TH2F("hEMPosition", "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1361 fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1363 fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
1364 fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
1367 // delete fhAllReg ;
1369 // delete fhPhotReg ;
1373 // delete fhNBarReg ;
1375 // delete fhChargedReg ;
1377 fhAllReg = new TH1F("hAllReg", "All primaries registered as photon", 100, 0., 5.);
1378 fhPhotReg = new TH1F("hPhotReg", "Photon registered as photon", 100, 0., 5.);
1379 fhNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.);
1380 fhNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.);
1381 fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
1386 // delete fhPhotEM ;
1390 // delete fhNBarEM ;
1392 // delete fhChargedEM ;
1394 fhAllEM = new TH1F("hAllEM", "All primary registered as EM",100, 0., 5.);
1395 fhPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
1396 fhNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
1397 fhNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
1398 fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
1401 // delete fhAllPPSD ;
1403 // delete fhPhotPPSD ;
1407 // delete fhNBarPPSD ;
1408 // if(fhChargedPPSD)
1409 // delete fhChargedPPSD ;
1411 fhAllPPSD = new TH1F("hAllPPSD", "All primary registered as PPSD",100, 0., 5.);
1412 fhPhotPPSD = new TH1F("hPhotPPSD", "Photon registered as PPSD", 100, 0., 5.);
1413 fhNPPSD = new TH1F("hNPPSD", "N registered as PPSD", 100, 0., 5.);
1414 fhNBarPPSD = new TH1F("hNBarPPSD", "NBar registered as PPSD", 100, 0., 5.);
1415 fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
1418 // delete fhPrimary ;
1419 fhPrimary= new TH1F("hPrimary", "hPrimary", 100, 0., 5.);
1430 fhAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
1431 fhVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.);
1432 fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
1433 fhPPSD = new TH1F("hPPSD", "All PPSD photon particles", 100, 0., 5.);
1437 // delete fhPhotPhot ;
1439 // delete fhPhotElec ;
1441 // delete fhPhotNeuH ;
1443 // delete fhPhotNuEM ;
1445 // delete fhPhotChHa ;
1447 // delete fhPhotGaHa ;
1449 fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.); //Photon registered as photon
1450 fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.); //Photon registered as Electron
1451 fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.); //Photon registered as Neutral Hadron
1452 fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.); //Photon registered as Neutral EM
1453 fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.); //Photon registered as Charged Hadron
1454 fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.); //Photon registered as Gamma-Hadron
1457 //____________________________________________________________________________
1458 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1460 // Open the root file named "name"
1462 fRootFile = new TFile(name, "update") ;
1463 return fRootFile->IsOpen() ;
1466 //____________________________________________________________________________
1467 void AliPHOSAnalyze::SaveHistograms()
1469 // Saves the histograms in a root file named "name.analyzed"
1471 Text_t outputname[80] ;
1472 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1473 TFile output(outputname,"RECREATE");
1477 fhAllEnergy->Write() ;
1479 fhPhotEnergy->Write() ;
1481 fhEMEnergy->Write() ;
1483 fhPPSDEnergy->Write() ;
1485 fhAllPosition->Write() ;
1487 fhAllPositionX->Write() ;
1489 fhAllPositionZ->Write() ;
1491 fhPhotPosition->Write() ;
1493 fhEMPosition->Write() ;
1495 fhPPSDPosition->Write() ;
1499 fhPhotReg->Write() ;
1503 fhNBarReg->Write() ;
1505 fhChargedReg->Write() ;
1515 fhChargedEM->Write() ;
1517 fhAllPPSD->Write() ;
1519 fhPhotPPSD->Write() ;
1523 fhNBarPPSD->Write() ;
1525 fhChargedPPSD->Write() ;
1527 fhPrimary->Write() ;
1537 fhPhotPhot->Write() ;
1539 fhPhotElec->Write() ;
1541 fhPhotNeuH->Write() ;
1543 fhPhotNuEM->Write() ;
1545 fhPhotNuEM->Write() ;
1547 fhPhotChHa->Write() ;
1549 fhPhotGaHa->Write() ;
1550 if(fhEnergyCorrelations)
1551 fhEnergyCorrelations->Write() ;
1556 //____________________________________________________________________________
1557 Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
1559 return ERecPart/0.8783 ;
1562 //____________________________________________________________________________
1563 void AliPHOSAnalyze::ResetHistograms()
1565 fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2)
1567 fhEmcDigit = 0 ; // Histo of digit energies in the Emc
1568 fhVetoDigit = 0 ; // Histo of digit energies in the Veto
1569 fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor
1570 fhEmcCluster = 0 ; // Histo of Cluster energies in Emc
1571 fhVetoCluster = 0 ; // Histo of Cluster energies in Veto
1572 fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor
1573 fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies
1576 fhPhotEnergy = 0 ; // Total spectrum of detected photons
1577 fhEMEnergy = 0 ; // Spectrum of detected electrons with electron primary
1580 fhAllPositionX = 0 ;
1581 fhAllPositionZ = 0 ;
1582 fhPhotPosition = 0 ;
1584 fhPPSDPosition = 0 ;