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 PHOS events. In this class we demostrate,
20 // how to handle reconstructed objects with AliPHSOIndexToObject.
21 // As an example we propose sulotions for four most frequently used tasks:
22 // DrawRecon(...) - to draw reconstructed objects in the PHOS plane,
23 // very usefull in the debuging
24 // InvarianMass(...) - to calculate "REAL" and "MIXED" photon pairs
25 // invariant mass distributions
26 // EnergyResoluition(...) -\ Energy and position resolutions of the
27 // PositionResolution(...)-/ reconstructed photons
28 // Contamination(...) - calculates contamination of the photon spectrum and
29 // pobability of reconstruction of several primaries as
30 // kGAMMA,kELECTRON etc.
33 // root [0] AliPHOSAnalyze * a = new AliPHOSAnalyze("galice.root")
34 // // set the file you want to analyse
35 // root [1] a->DrawRecon(1,3)
36 // // plot RecObjects, made in event 1, PHOS module 3
37 // root [2] a->DrawRecon(1,3,"PHOSRP","another PID")
38 // // plot RecObjets made in the event 1, PHOS module 3,
39 // // produced in the another reconstruction pass,
40 // // which produced PHOS RecParticles ("PHOSRP") with
41 // // title "another PID".
42 // root [3] a->InvariantMass()
43 // // Calculates "REAL" and "MIXED" invariant mass
44 // // distributions of kGAMMA and (kGAMMA+kNEUTRALEM)
45 // // and APPENDS this to the file "invmass.root"
46 // root [4] a->PositionResolution()
47 // // calculates two dimentional histos: energy of the primary
48 // // photon vs distance betwin incedence point and reconstructed
49 // // poisition. One can analyse the produced file position.root
50 // // with macro PhotonPosition.C
51 // root [5] a->EnergyResolution()
52 // // calculates two dimentional histos: energy of the primary
53 // // photon vs energy of the reconstructed particle. One can
54 // // analyse the produced file energy.root
55 // // with macro PhotonEnergy.C
56 // root [6] a->Contamination()
57 // // fills spectra of primary photons and several kinds of
58 // // reconstructed particles, so that analyzing them one can
59 // // estimate conatmination, efficiency of registration etc.
61 //*-- Author: Dmitri Peressounko (SUBATECH & RRC Kurchatov Institute)
62 //////////////////////////////////////////////////////////////////////////////
65 // --- ROOT system ---
71 #include "TParticle.h"
72 #include "TClonesArray.h"
78 // --- Standard library ---
82 // --- AliRoot header files ---
85 #include "AliPHOSv1.h"
86 #include "AliPHOSAnalyze.h"
87 #include "AliPHOSDigit.h"
88 #include "AliPHOSSDigitizer.h"
89 #include "AliPHOSTrackSegment.h"
90 #include "AliPHOSRecParticle.h"
91 #include "AliPHOSCpvRecPoint.h"
92 #include "AliPHOSGetter.h"
95 ClassImp(AliPHOSAnalyze)
97 //____________________________________________________________________________
98 AliPHOSAnalyze::AliPHOSAnalyze()
100 // default ctor (useless)
101 fCorrection = 1.2 ; //Value calculated for default parameters of reconstruction
104 //____________________________________________________________________________
105 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName)
107 // ctor: analyze events from root file "name"
108 ffileName = fileName ;
109 fCorrection = 1.05 ; //Value calculated for default parameters of reconstruction
112 //____________________________________________________________________________
113 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
116 ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
119 //____________________________________________________________________________
120 AliPHOSAnalyze::~AliPHOSAnalyze()
125 //____________________________________________________________________________
126 void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,const char* branchTitle){
127 //Draws pimary particles and reconstructed
128 //digits, RecPoints, RecPartices etc
129 //for event Nevent in the module Nmod.
131 //========== Create ObjectGetter
132 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
133 const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
136 Int_t nx = phosgeom->GetNPhi() ;
137 Int_t nz = phosgeom->GetNZ() ;
138 Float_t * cri= phosgeom->GetEMCAGeometry()->GetCrystalHalfSize() ;
139 Float_t x = nx*cri[0] ;
140 Float_t z = nz*cri[2] ;
141 Int_t nxCPV = (Int_t) (nx*phosgeom->GetPadSizePhi()/(2.*cri[0])) ;
142 Int_t nzCPV = (Int_t) (nz*phosgeom->GetPadSizeZ()/(2.*cri[2])) ;
144 TH2F * emcDigits = new TH2F("emcDigits","EMC digits", nx,-x,x,nz,-z,z);
145 TH2F * emcSdigits = new TH2F("emcSdigits","EMC sdigits", nx,-x,x,nz,-z,z);
146 TH2F * emcRecPoints = new TH2F("emcRecPoints","EMC RecPoints",nx,-x,x,nz,-z,z);
147 TH2F * cpvSdigits = new TH2F("cpvSdigits","CPV sdigits", nx,-x,x,nz,-z,z);
148 TH2F * cpvDigits = new TH2F("cpvDigits","CPV digits", nxCPV,-x,x,nzCPV,-z,z) ;
149 TH2F * cpvRecPoints = new TH2F("cpvCl","CPV RecPoints", nxCPV,-x,x,nzCPV,-z,z) ;
150 TH2F * nbar = new TH2F("nbar","Primary nbar", nx,-x,x,nz,-z,z);
151 TH2F * phot = new TH2F("phot","Primary Photon", nx,-x,x,nz,-z,z);
152 TH2F * charg = new TH2F("charg","Primary charged",nx,-x,x,nz,-z,z);
153 TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",nx,-x,x,nz,-z,z);
154 TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", nx,-x,x,nz,-z,z);
157 //Plot Primary Particles
158 const TParticle * primary ;
160 for ( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++)
162 primary = gime->Primary(iPrimary) ;
163 Int_t primaryType = primary->GetPdgCode() ;
164 if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212)
165 ||(primaryType == 11)||(primaryType == -11) ) {
167 Double_t primX, primZ ;
168 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
169 if(moduleNumber==Nmod)
170 charg->Fill(primZ,primX,primary->Energy()) ;
172 if( primaryType == 22 ) {
174 Double_t primX, primZ ;
175 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
176 if(moduleNumber==Nmod)
177 phot->Fill(primZ,primX,primary->Energy()) ;
180 if( primaryType == -2112 ) {
182 Double_t primX, primZ ;
183 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
184 if(moduleNumber==Nmod)
185 nbar->Fill(primZ,primX,primary->Energy()) ;
192 AliPHOSDigit * sdigit ;
193 TClonesArray * sdigits = gime->SDigits() ;
195 for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++)
197 sdigit = (AliPHOSDigit *) sdigits->At(iSDigit) ;
199 phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
201 phosgeom->RelPosInModule(relid,x,z) ;
202 Float_t e = gime->SDigitizer()->Calibrate(sdigit->GetAmp()) ;
204 if(relid[1]==0) //EMC
205 emcSdigits->Fill(x,z,e) ;
207 cpvSdigits->Fill(x,z,e) ;
214 AliPHOSDigit * digit ;
215 TClonesArray * digits = gime->Digits();
217 for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++)
219 digit = (AliPHOSDigit *) digits->At(iDigit) ;
221 phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
223 phosgeom->RelPosInModule(relid,x,z) ;
224 Float_t e = gime->SDigitizer()->Calibrate(digit->GetAmp()) ;
226 if(relid[1]==0) //EMC
227 emcDigits->Fill(x,z,e) ;
229 cpvDigits->Fill(x,z,e) ;
238 TObjArray * emcrp = gime->EmcRecPoints() ;
240 for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){
241 AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
242 if(emc->GetPHOSMod()==Nmod){
243 emc->GetLocalPosition(pos) ;
244 emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy());
249 TObjArray * cpvrp = gime->CpvRecPoints() ;
251 for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){
252 AliPHOSRecPoint * cpv = (AliPHOSCpvRecPoint *) cpvrp->At(irecp) ;
253 if(cpv->GetPHOSMod()==Nmod){
254 cpv->GetLocalPosition(pos) ;
255 cpvRecPoints->Fill(pos.X(),pos.Z(),cpv->GetEnergy());
261 AliPHOSRecParticle * recParticle ;
263 TClonesArray * rp = gime->RecParticles() ;
264 TClonesArray * ts = gime->TrackSegments() ;
265 if(rp && ts && emcrp) {
266 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ )
268 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
269 Int_t moduleNumberRec ;
270 Double_t recX, recZ ;
271 phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
272 if(moduleNumberRec == Nmod){
274 Double_t minDistance = 5. ;
275 Int_t closestPrimary = -1 ;
277 //extract list of primaries: it is stored at EMC RecPoints
278 Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
279 Int_t numberofprimaries ;
280 Int_t * listofprimaries = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
282 const TParticle * primary ;
283 Double_t distance = minDistance ;
285 for ( index = 0 ; index < numberofprimaries ; index++){
286 primary = gime->Primary(listofprimaries[index]) ;
288 Double_t primX, primZ ;
289 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
290 if(moduleNumberRec == moduleNumber)
291 distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
292 if(minDistance > distance)
294 minDistance = distance ;
295 closestPrimary = listofprimaries[index] ;
299 if(closestPrimary >=0 ){
301 Int_t primaryType = gime->Primary(closestPrimary)->GetPdgCode() ;
304 recPhot->Fill(recZ,recX,recParticle->Energy()) ;
306 if(primaryType==-2112)
307 recNbar->Fill(recZ,recX,recParticle->Energy()) ;
314 //Plot made histograms
315 emcSdigits->Draw("box") ;
316 emcDigits->SetLineColor(5) ;
317 emcDigits->Draw("box") ;
318 emcRecPoints->SetLineColor(2) ;
319 emcRecPoints->Draw("boxsame") ;
320 cpvSdigits->SetLineColor(1) ;
321 cpvSdigits->Draw("box") ;
322 charg->SetLineColor(2) ;
323 charg->Draw("boxsame") ;
324 nbar->SetLineColor(6) ;
325 nbar->Draw("boxsame") ;
328 //____________________________________________________________________________
329 void AliPHOSAnalyze::Ls(){
330 //lists branches and titles of PHOS-related branches of TreeR, TreeD, TreeS
332 AliPHOSGetter::GetInstance(ffileName.Data()) ;
335 TObjArray * branches;
337 branches = gAlice->TreeS()->GetListOfBranches() ;
339 cout << "TreeS: " << endl ;
340 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
341 TBranch * branch=(TBranch *) branches->At(ibranch) ;
342 if(strstr(branch->GetName(),"PHOS") )
343 cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ;
347 branches = gAlice->TreeD()->GetListOfBranches() ;
349 cout << "TreeD: " << endl ;
350 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
351 TBranch * branch=(TBranch *) branches->At(ibranch) ;
352 if(strstr(branch->GetName(),"PHOS") )
353 cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ;
358 branches = gAlice->TreeR()->GetListOfBranches() ;
360 cout << "TreeR: " << endl ;
361 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
362 TBranch * branch=(TBranch *) branches->At(ibranch) ;
363 if(strstr(branch->GetName(),"PHOS") )
364 cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ;
370 //____________________________________________________________________________
371 void AliPHOSAnalyze::InvariantMass(const char* branchTitle)
373 // Calculates Real and Mixed invariant mass distributions
374 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
376 Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
379 TFile * mfile = new TFile("invmass.root","update");
381 //========== Reading /Booking Histograms
383 hRealEM = (TH2D*) mfile->Get("hRealEM") ;
385 hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ;
386 TH2D * hRealPhot = 0 ;
388 hRealPhot = (TH2D*)mfile->Get("hRealPhot");
390 hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
392 TH2D * hMixedEM = 0 ;
393 hMixedEM = (TH2D*) mfile->Get("hMixedEM") ;
395 hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ;
397 TH2D * hMixedPhot = 0 ;
398 hMixedPhot = (TH2D*) mfile->Get("hMixedPhot") ;
400 hMixedPhot = new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
403 //reading event and copyng it to TConesArray of all photons
405 TClonesArray * allRecParticleList = new TClonesArray("AliPHOSRecParticle", 1000) ;
406 Int_t * nRecParticles = new Int_t[nMixedEvents] ; // to mark boundaries of each event in the total list
407 for(Int_t index = 0; index < nMixedEvents; index ++)
408 nRecParticles[index] = 0 ;
409 Int_t iRecPhot = 0 ; // number of EM particles in total list
411 //scan over all events
413 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
414 // for(event = 0; event < gime->MaxEvent(); event++ ){
415 for(event = 0; event < maxevent; event++ ){
416 gime->Event(event,"R"); //will read only TreeR
418 //copy EM RecParticles to the "total" list
419 const AliPHOSRecParticle * recParticle ;
421 TClonesArray * rp = gime->RecParticles() ;
423 cout << "AliPHOSAnalyze::InvariantMass --> Can't find RecParticles " << endl ;
427 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ )
429 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
430 if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
431 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW))
432 new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
435 Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle
436 nRecParticles[mevent] = iRecPhot-1 ;
438 //check, if it is time to calculate invariant mass?
439 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
440 if((mevent == 0) && (event +1 == maxevent)){
442 // if((mevent == 0) && (event +1 == gime->MaxEvent())){
444 //calculate invariant mass:
446 Int_t nCurEvent = 0 ;
448 for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
449 AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
451 for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
452 AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
455 invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
456 (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
457 (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
458 (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
461 invMass = TMath::Sqrt(invMass);
464 pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +
465 (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) );
467 if(irp1 > nRecParticles[nCurEvent])
470 if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
471 hRealEM->Fill(invMass,pt);
472 if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
473 (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
474 hRealPhot->Fill(invMass,pt);
477 hMixedEM->Fill(invMass,pt);
478 if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
479 (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
480 hMixedPhot->Fill(invMass,pt);
483 } //loop over second rp
484 }//loop over first rp
486 //Make some cleanings
487 for(Int_t index = 0; index < nMixedEvents; index ++)
488 nRecParticles[index] = 0 ;
490 allRecParticleList->Clear() ;
494 delete allRecParticleList ;
499 hRealEM->Write(0,kOverwrite) ;
500 hRealPhot->Write(0,kOverwrite) ;
501 hMixedEM->Write(0,kOverwrite) ;
502 hMixedPhot->Write(0,kOverwrite) ;
507 delete nRecParticles;
511 //____________________________________________________________________________
512 void AliPHOSAnalyze::EnergyResolution(const char * branchTitle)
514 //fills two dimentional histo: energy of primary vs. energy of reconstructed
516 TH2F * hAllEnergy = 0 ; //all reconstructed with primary photon
517 TH2F * hPhotEnergy= 0 ; //kGamma with primary photon
518 TH2F * hEMEnergy = 0 ; //electromagnetic with primary photon
520 //opening file and reading histograms if any
521 TFile * efile = new TFile("energy.root","update");
523 hAllEnergy = (TH2F*)efile->Get("hAllEnergy") ;
525 hAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
527 hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ;
529 hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
531 hEMEnergy =(TH2F*) efile->Get("hEMEnergy");
533 hEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.);
536 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
537 const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
540 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
541 for ( ievent=0; ievent < maxevent ; ievent++){
543 //read the current event
544 gime->Event(ievent) ;
546 const AliPHOSRecParticle * recParticle ;
548 TClonesArray * rp = gime->RecParticles() ;
550 cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
551 << ", Can't find RecParticles " << endl ;
554 TClonesArray * ts = gime->TrackSegments() ;
556 cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
557 << ", Can't find TrackSegments " << endl ;
560 TObjArray * emcrp = gime->EmcRecPoints() ;
562 cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent
563 << ", Can't find EmcRecPoints " << endl ;
567 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){
568 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
570 //find the closest primary
571 Int_t moduleNumberRec ;
572 Double_t recX, recZ ;
573 phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
575 Double_t minDistance = 100. ;
576 Int_t closestPrimary = -1 ;
578 //extract list of primaries: it is stored at EMC RecPoints
579 Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
580 Int_t numberofprimaries ;
581 Int_t * listofprimaries = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
584 const TParticle * primary ;
585 Double_t distance = minDistance ;
588 Double_t dZmin = 0. ;
589 for ( index = 0 ; index < numberofprimaries ; index++){
590 primary = gime->Primary(listofprimaries[index]) ;
592 Double_t primX, primZ ;
593 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
594 if(moduleNumberRec == moduleNumber) {
597 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
598 if(minDistance > distance) {
599 minDistance = distance ;
602 closestPrimary = listofprimaries[index] ;
607 //if found primary, fill histograms
608 if(closestPrimary >=0 ){
609 const TParticle * primary = gime->Primary(closestPrimary) ;
610 if(primary->GetPdgCode() == 22){
611 hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
612 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
613 hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
614 hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
617 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
618 hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ;
624 //write filled histograms
626 hAllEnergy->Write(0,kOverwrite) ;
627 hPhotEnergy->Write(0,kOverwrite) ;
628 hEMEnergy->Write(0,kOverwrite) ;
634 //____________________________________________________________________________
635 void AliPHOSAnalyze::PositionResolution(const char * branchTitle)
637 //fills two dimentional histo: energy vs. primary - reconstructed distance
641 TH2F * hAllPosition = 0; // Position of any RP with primary photon
642 TH2F * hPhotPosition = 0; // Position of kGAMMA with primary photon
643 TH2F * hEMPosition = 0; // Position of EM with primary photon
645 TH1F * hAllPositionX = 0; // X-Position Resolution of photons with photon primary
646 TH1F * hAllPositionZ = 0; // Z-Position Resolution of photons with photon primary
649 //opening file and reading histograms if any
650 TFile * pfile = new TFile("position.root","update");
652 hAllPosition = (TH2F*)pfile->Get("hAllPosition");
653 if(hAllPosition == 0)
654 hAllPosition = new TH2F("hAllPosition",
655 "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
656 hPhotPosition= (TH2F*)pfile->Get("hPhotPosition");
657 if(hPhotPosition == 0)
658 hPhotPosition = new TH2F("hPhotPosition",
659 "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
660 hEMPosition= (TH2F*)pfile->Get("hEMPosition") ;
662 hEMPosition = new TH2F("hEMPosition",
663 "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.);
664 hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ;
665 if(hAllPositionX == 0)
666 hAllPositionX = new TH1F("hAllPositionX",
667 "Delta X of any RP with primary photon",100, -2., 2.);
668 hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ;
669 if(hAllPositionZ == 0)
670 hAllPositionZ = new TH1F("hAllPositionZ",
671 "Delta X of any RP with primary photon",100, -2., 2.);
674 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
675 const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
678 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
679 for ( ievent=0; ievent < maxevent ; ievent++){
681 //read the current event
682 gime->Event(ievent) ;
683 TClonesArray * rp = gime->RecParticles() ;
685 cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
686 << ", Can't find RecParticles " << endl ;
689 TClonesArray * ts = gime->TrackSegments() ;
691 cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
692 << ", Can't find TrackSegments " << endl ;
695 TObjArray * emcrp = gime->EmcRecPoints() ;
697 cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent
698 << ", Can't find EmcRecPoints " << endl ;
703 const AliPHOSRecParticle * recParticle ;
705 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
706 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
708 //find the closest primary
709 Int_t moduleNumberRec ;
710 Double_t recX, recZ ;
711 phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
713 Double_t minDistance = 100. ;
714 Int_t closestPrimary = -1 ;
716 //extract list of primaries: it is stored at EMC RecPoints
717 Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
718 Int_t numberofprimaries ;
719 Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
722 const TParticle * primary ;
723 Double_t distance = minDistance ;
724 Double_t dX = 1000; // incredible number
725 Double_t dZ = 1000; // for the case if no primary will be found
727 Double_t dZmin = 0. ;
728 for ( index = 0 ; index < numberofprimaries ; index++){
729 primary = gime->Primary(listofprimaries[index]) ;
731 Double_t primX, primZ ;
732 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
733 if(moduleNumberRec == moduleNumber) {
736 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
737 if(minDistance > distance) {
738 minDistance = distance ;
741 closestPrimary = listofprimaries[index] ;
746 //if found primary, fill histograms
747 if(closestPrimary >=0 ){
748 const TParticle * primary = gime->Primary(closestPrimary) ;
749 if(primary->GetPdgCode() == 22){
750 hAllPosition->Fill(primary->Energy(), minDistance) ;
751 hAllPositionX->Fill(primary->Energy(), dX) ;
752 hAllPositionZ->Fill(primary->Energy(), dZ) ;
753 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
754 hPhotPosition->Fill(primary->Energy(), minDistance ) ;
755 hEMPosition->Fill(primary->Energy(), minDistance ) ;
758 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
759 hEMPosition->Fill(primary->Energy(), minDistance ) ;
765 //Write output histgrams
767 hAllPosition->Write(0,kOverwrite) ;
768 hAllPositionX->Write(0,kOverwrite) ;
769 hAllPositionZ->Write(0,kOverwrite) ;
770 hPhotPosition->Write(0,kOverwrite) ;
771 hEMPosition->Write(0,kOverwrite) ;
778 //____________________________________________________________________________
779 void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
780 // fills spectra of primary photons and several kinds of
781 // reconstructed particles, so that analyzing them one can
782 // estimate conatmination, efficiency of registration etc.
784 //define several general histograms
785 TH1F * hPrimary = 0; //spectrum (P_t distribution) of primary photons
786 TH1F * hAllRP = 0; //spectrum of all RecParticles in PHOS
787 TH1F * hPhot = 0; //spectrum of kGAMMA RecParticles
788 TH1F * hShape = 0; //spectrum of all EM RecParticles
789 TH1F * hVeto = 0; //spectrum of all neutral RecParticles
791 //Now separate histograms in accoradance with primary
793 TH1F * hPhotReg = 0; //Registeres as photon
794 TH1F * hPhotEM = 0; //Registered as EM
797 TH1F * hNReg = 0; //Registeres as photon
798 TH1F * hNEM = 0; //Registered as EM
801 TH1F * hNBarReg = 0; //Registeres as photon
802 TH1F * hNBarEM = 0; //Registered as EM
804 //primary - charged hadron (pBar excluded)
805 TH1F * hChargedReg = 0; //Registeres as photon
806 TH1F * hChargedEM = 0; //Registered as EM
809 TH1F * hPbarReg = 0; //Registeres as photon
810 TH1F * hPbarEM = 0; //Registered as EM
813 //Reading histograms from the file
814 TFile * cfile = new TFile("contamination.root","update") ;
816 //read general histograms
817 hPrimary = (TH1F*) cfile->Get("hPrimary") ;
819 hPrimary= new TH1F("hPrimary", "Primary photon spectrum", 100, 0., 5.);
820 hAllRP = (TH1F*)cfile->Get("hAllRP") ;
822 hAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
823 hPhot = (TH1F*)cfile->Get("hPhot") ;
825 hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.);
826 hShape = (TH1F*) cfile->Get("hShape") ;
828 hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.);
829 hVeto= (TH1F*)cfile->Get("hVeto") ;
831 hVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.);
835 hPhotReg = (TH1F*)cfile->Get("hPhotReg");
837 hPhotReg = new TH1F("hPhotReg","Photon registered as photon",100, 0., 5.);
838 hPhotEM =(TH1F*)cfile->Get("hPhotEM");
840 hPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
843 hNReg = (TH1F*)cfile->Get("hNReg");
845 hNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.);
846 hNEM = (TH1F*)cfile->Get("hNEM");
848 hNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
851 hNBarReg =(TH1F*)cfile->Get("hNBarReg");
853 hNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.);
854 hNBarEM =(TH1F*)cfile->Get("hNBarEM");
856 hNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
858 //primary - charged hadron (pBar excluded)
859 hChargedReg = (TH1F*)cfile->Get("hChargedReg");
861 hChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
862 hChargedEM = (TH1F*)cfile->Get("hChargedEM");
864 hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
867 hPbarReg = (TH1F*)cfile->Get("hPbarReg");
869 hPbarReg= new TH1F("hPbarReg", "pBar registered as photon",100, 0., 5.);
870 hPbarEM = (TH1F*)cfile->Get("hPbarEM");
872 hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.);
875 //Now make some initializations
877 Int_t counter[8][5] ; //# of registered particles
879 for(i1 = 0; i1<8; i1++)
880 for(i2 = 0; i2<5; i2++)
881 counter[i1][i2] = 0 ;
885 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),RecPointsTitle) ;
886 const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ;
889 Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
890 for ( ievent=0; ievent < maxevent ; ievent++){
892 gime->Event(ievent) ;
894 TClonesArray * rp = gime->RecParticles() ;
896 cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
897 << ", Can't find RecParticles " << endl ;
900 TClonesArray * ts = gime->TrackSegments() ;
902 cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
903 << ", Can't find TrackSegments " << endl ;
906 TObjArray * emcrp = gime->EmcRecPoints() ;
908 cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent
909 << ", Can't find EmcRecPoints " << endl ;
914 //=========== Make spectrum of the primary photons
915 const TParticle * primary ;
917 for( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++){
918 primary = gime->Primary(iPrimary) ;
919 Int_t primaryType = primary->GetPdgCode() ;
920 if( primaryType == 22 ) {
921 //check, if photons folls onto PHOS
923 Double_t primX, primZ ;
924 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
926 hPrimary->Fill(primary->Energy()) ;
932 //========== Now scan over RecParticles
933 const AliPHOSRecParticle * recParticle ;
935 for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
936 recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
937 //fill histo spectrum of all RecParticles
938 hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
940 //==========find the closest primary
941 Int_t moduleNumberRec ;
942 Double_t recX, recZ ;
943 phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
945 Double_t minDistance = 100. ;
946 Int_t closestPrimary = -1 ;
948 //extract list of primaries: it is stored at EMC RecPoints
949 Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
950 Int_t numberofprimaries ;
951 Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ;
953 const TParticle * primary ;
954 Double_t distance = minDistance ;
957 Double_t dZmin = 0. ;
958 for ( index = 0 ; index < numberofprimaries ; index++){
959 primary = gime->Primary(listofprimaries[index]) ;
961 Double_t primX, primZ ;
962 phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
963 if(moduleNumberRec == moduleNumber) {
966 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
967 if(minDistance > distance) {
968 minDistance = distance ;
971 closestPrimary = listofprimaries[index] ;
976 //===========define the "type" of closest primary
977 if(closestPrimary >=0 ){
978 Int_t primaryCode = -1;
979 const TParticle * primary = gime->Primary(closestPrimary) ;
980 Int_t primaryType = primary->GetPdgCode() ;
981 if(primaryType == 22) // photon ?
984 if(primaryType == 2112) // neutron
987 if(primaryType == -2112) // Anti neutron
990 if(primaryType == -2122) //Anti proton
993 TParticle tempo(*primary) ;
994 if(tempo.GetPDG()->Charge())
998 //==========Now look at the type of RecParticle
999 Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
1000 if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
1001 hPhot->Fill(energy ) ;
1002 switch(primaryCode){
1004 hPhotReg->Fill(energy ) ;
1007 hNReg->Fill(energy ) ;
1010 hNBarReg->Fill(energy ) ;
1013 hChargedReg->Fill(energy ) ;
1016 hPbarReg->Fill(energy ) ;
1022 if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
1023 (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)||
1024 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)||
1025 (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower
1026 hShape->Fill(energy ) ;
1027 switch(primaryCode){
1029 hPhotEM->Fill(energy ) ;
1032 hNEM->Fill(energy ) ;
1035 hNBarEM->Fill(energy ) ;
1038 hChargedEM->Fill(energy ) ;
1041 hPbarEM->Fill(energy ) ;
1048 if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
1049 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) ||
1050 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) ||
1051 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral
1052 hVeto->Fill(energy ) ;
1054 //fill number of primaries identified as ...
1055 if(primaryCode >= 0) // Primary code defined
1056 counter[recParticle->GetType()][primaryCode]++ ;
1060 } // no closest primary found
1064 //=================== SaveHistograms
1066 hPrimary->Write(0,kOverwrite);
1067 hAllRP->Write(0,kOverwrite);
1068 hPhot->Write(0,kOverwrite);
1069 hShape->Write(0,kOverwrite);
1070 hVeto->Write(0,kOverwrite);
1071 hPhotReg->Write(0,kOverwrite);
1072 hPhotEM->Write(0,kOverwrite);
1073 hNReg ->Write(0,kOverwrite);
1074 hNEM ->Write(0,kOverwrite);
1075 hNBarReg ->Write(0,kOverwrite);
1076 hNBarEM ->Write(0,kOverwrite);
1077 hChargedReg ->Write(0,kOverwrite);
1078 hChargedEM ->Write(0,kOverwrite);
1079 hPbarReg ->Write(0,kOverwrite);
1080 hPbarEM ->Write(0,kOverwrite);
1082 cfile->Write(0,kOverwrite);
1088 maxevent = (Int_t)gAlice->TreeE()->GetEntries() ;
1089 // cout << "Resolutions: Analyzed " << gime->MaxEvent() << " event(s)" << endl ;
1091 cout << "Resolutions: Analyzed " << maxevent << " event(s)" << endl ;
1094 cout << " Primary: Photon Neutron Antineutron Charged hadron AntiProton" << endl ;
1095 cout << "--------------------------------------------------------------------------------" << endl ;
1097 << setw(8) << counter[2][0] << setw(9) << counter[2][1] << setw(13) << counter[2][2]
1098 << setw(15)<< counter[2][3] << setw(13) << counter[2][4] << endl ;
1099 cout << " kGAMMAHA: "
1100 << setw(8) << counter[3][0] << setw(9) << counter[3][1] << setw(13) << counter[3][2]
1101 << setw(15)<< counter[3][3] << setw(13) << counter[3][4] << endl ;
1102 cout << " kNEUTRALEM: "
1103 << setw(8) << counter[0][0] << setw(9) << counter[0][1] << setw(13) << counter[0][2]
1104 << setw(15)<< counter[0][3] << setw(13) << counter[0][4] << endl ;
1105 cout << " kNEUTRALHA: "
1106 << setw(8) << counter[1][0] << setw(9) << counter[1][1] << setw(13) << counter[1][2]
1107 << setw(15)<< counter[1][3] << setw(13) << counter[1][4] << endl ;
1108 cout << " kABSURDEM: "
1109 << setw(8) << counter[4][0] << setw(9) << counter[4][1] << setw(13) << counter[4][2]
1110 << setw(15)<< counter[4][3] << setw(13) << counter[4][4] << endl ;
1111 cout << " kABSURDHA: "
1112 << setw(8) << counter[5][0] << setw(9) << counter[5][1] << setw(13) << counter[5][2]
1113 << setw(15)<< counter[5][3] << setw(13) << counter[5][4] << endl ;
1114 cout << " kELECTRON: "
1115 << setw(8) << counter[6][0] << setw(9) << counter[6][1] << setw(13) << counter[6][2]
1116 << setw(15)<< counter[6][3] << setw(13) << counter[6][4] << endl ;
1117 cout << " kCHARGEDHA: "
1118 << setw(8) << counter[7][0] << setw(9) << counter[7][1] << setw(13) << counter[7][2]
1119 << setw(15)<< counter[7][3] << setw(13) << counter[7][4] << endl ;
1120 cout << "--------------------------------------------------------------------------------" << endl ;
1122 Int_t totalInd = 0 ;
1123 for(i1 = 0; i1<8; i1++)
1124 for(i2 = 0; i2<5; i2++)
1125 totalInd+=counter[i1][i2] ;
1126 cout << "Indentified particles: " << totalInd << endl ;