/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* Author: The ALICE Off-line Project. *
* Contributors are mentioned in the code where appropriate. *
* *
* Permission to use, copy, modify and distribute this software and its *
* documentation strictly for non-commercial purposes is hereby granted *
* without fee, provided that the above copyright notice appears in all *
* copies and that both the copyright notice and this permission notice *
* appear in the supporting documentation. The authors make no claims *
* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
/* $Id$ */
//_________________________________________________________________________
// Algorythm class to analyze PHOSv1 events:
// Construct histograms and displays them.
// Use the macro EditorBar.C for best access to the functionnalities
//*--
//*-- Author: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH)
//////////////////////////////////////////////////////////////////////////////
// --- ROOT system ---
#include "TFile.h"
#include "TH1.h"
#include "TPad.h"
#include "TH2.h"
#include "TH2.h"
#include "TParticle.h"
#include "TClonesArray.h"
#include "TTree.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TStyle.h"
// --- Standard library ---
#include <iostream.h>
#include <stdio.h>
// --- AliRoot header files ---
#include "AliRun.h"
#include "AliPHOSv1.h"
#include "AliPHOSAnalyze.h"
#include "AliPHOSClusterizerv1.h"
#include "AliPHOSTrackSegmentMakerv1.h"
#include "AliPHOSPIDv1.h"
#include "AliPHOSReconstructioner.h"
#include "AliPHOSDigit.h"
#include "AliPHOSTrackSegment.h"
#include "AliPHOSRecParticle.h"
#include "AliPHOSIndexToObject.h"
#include "AliPHOSHit.h"
#include "AliPHOSCpvRecPoint.h"
#include "AliPHOSPpsdRecPoint.h"
ClassImp(AliPHOSAnalyze)
//____________________________________________________________________________
AliPHOSAnalyze::AliPHOSAnalyze()
{
// default ctor (useless)
fRootFile = 0 ;
}
//____________________________________________________________________________
AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
{
// ctor: analyze events from root file "name"
Bool_t ok = OpenRootFile(name) ;
if ( !ok ) {
cout << " AliPHOSAnalyze > Error opening " << name << endl ;
}
else {
//========== Get AliRun object from file
gAlice = (AliRun*) fRootFile->Get("gAlice") ;
//=========== Get the PHOS object and associated geometry from the file
fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
//========== Initializes the Index to Object converter
fObjGetter = AliPHOSIndexToObject::GetInstance() ; // <--- To be redone
//========== Current event number
fEvt = -999 ;
}
fDebugLevel = 0;
// fClu = 0 ;
// fPID = 0 ;
// fTrs = 0 ;
// fRec = 0 ;
ResetHistograms() ;
}
//____________________________________________________________________________
AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
{
// copy ctor
( (AliPHOSAnalyze &)ana ).Copy(*this) ;
}
//____________________________________________________________________________
void AliPHOSAnalyze::Copy(TObject & obj)
{
// copy an analysis into an other one
TObject::Copy(obj) ;
// I do nothing more because the copy is silly but the Code checkers requires one
}
//____________________________________________________________________________
AliPHOSAnalyze::~AliPHOSAnalyze()
{
// dtor
if(fRootFile->IsOpen()) fRootFile->Close() ;
if(fRootFile) {delete fRootFile ; fRootFile=0 ;}
if(fPHOS) {delete fPHOS ; fPHOS =0 ;}
// if(fClu) {delete fClu ; fClu =0 ;}
// if(fPID) {delete fPID ; fPID =0 ;}
// if(fRec) {delete fRec ; fRec =0 ;}
// if(fTrs) {delete fTrs ; fTrs =0 ;}
}
//____________________________________________________________________________
void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
//Draws pimary particles and reconstructed
//digits, RecPoints, RecPartices etc
//for event Nevent in the module Nmod.
TH2F * digitOccupancy = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.);
TH2F * sdigitOccupancy = new TH2F("sdigitOccupancy","EMC sdigits", 64,-71.,71.,64,-71.,71.);
TH2F * emcOccupancy = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.);
TH2F * ppsdUp = new TH2F("ppsdUp","PPSD Up digits", 128,-71.,71.,128,-71.,71.) ;
TH2F * ppsdUpCl = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ;
TH2F * ppsdLow = new TH2F("ppsdLow","PPSD Low digits", 128,-71.,71.,128,-71.,71.) ;
TH2F * ppsdLowCl = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ;
TH2F * nbar = new TH2F("nbar","Primary nbar", 64,-71.,71.,64,-71.,71.);
TH2F * phot = new TH2F("phot","Primary Photon", 64,-71.,71.,64,-71.,71.);
TH2F * charg = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.);
TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.);
TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", 64,-71.,71.,64,-71.,71.);
//========== Create the Clusterizer
// fClu = new AliPHOSClusterizerv1() ;
gAlice->GetEvent(Nevent);
TParticle * primary ;
Int_t iPrimary ;
for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
{
primary = gAlice->Particle(iPrimary) ;
Int_t primaryType = primary->GetPdgCode() ;
if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
Int_t moduleNumber ;
Double_t primX, primZ ;
fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
if(moduleNumber==Nmod)
charg->Fill(primZ,primX,primary->Energy()) ;
}
if( primaryType == 22 ) {
Int_t moduleNumber ;
Double_t primX, primZ ;
fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
if(moduleNumber==Nmod)
phot->Fill(primZ,primX,primary->Energy()) ;
}
else{
if( primaryType == -2112 ) {
Int_t moduleNumber ;
Double_t primX, primZ ;
fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
if(moduleNumber==Nmod)
nbar->Fill(primZ,primX,primary->Energy()) ;
}
}
}
//Set TreeS here and get AliPHOSSdigitizer
gAlice->TreeS()->GetEvent(0) ;
Int_t iSDigit ;
AliPHOSDigit * sdigit ;
if(fPHOS->SDigits()){
for(iSDigit = 0; iSDigit < fPHOS->SDigits()->GetEntries(); iSDigit++)
{
sdigit = (AliPHOSDigit *) fPHOS->SDigits()->At(iSDigit) ;
Int_t relid[4];
fGeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
Float_t x,z ;
fGeom->RelPosInModule(relid,x,z) ;
Float_t e = 1 ; //<--- fPHOS->Calibrate(sdigit->GetAmp()) ;
if(relid[0]==Nmod){
if(relid[1]==0) //EMC
sdigitOccupancy->Fill(x,z,e) ;
if((relid[1]>0)&&(relid[1]<17))
ppsdUp->Fill(x,z,e) ;
if(relid[1]>16)
ppsdLow->Fill(x,z,e) ;
}
}
}
else{
cout << "No SDigits read " << endl ;
}
gAlice->TreeD()->GetEvent(0) ;
if(fPHOS->Digits()){
Int_t iDigit ;
AliPHOSDigit * digit ;
for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
{
digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
Int_t relid[4];
fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
Float_t x,z ;
fGeom->RelPosInModule(relid,x,z) ;
Float_t e = 1; //<--- fClu->Calibrate(digit->GetAmp()) ;
if(relid[0]==Nmod){
if(relid[1]==0) //EMC
digitOccupancy->Fill(x,z,e) ;
if((relid[1]>0)&&(relid[1]<17))
ppsdUp->Fill(x,z,e) ;
if(relid[1]>16)
ppsdLow->Fill(x,z,e) ;
}
}
}
else{
cout << "No Digits read " << endl ;
}
gAlice->TreeR()->GetEvent(0) ;
TObjArray * emcRecPoints = fPHOS->EmcRecPoints() ;
TObjArray * ppsdRecPoints = fPHOS->PpsdRecPoints() ;
TClonesArray * recParticleList = fPHOS->RecParticles() ;
Int_t irecp ;
TVector3 pos ;
if(emcRecPoints ){
for(irecp = 0; irecp < emcRecPoints->GetEntries() ; irecp ++){
AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)emcRecPoints->At(irecp) ;
if(emc->GetPHOSMod()==Nmod){
emc->GetLocalPosition(pos) ;
emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
}
}
}
else{
cout << "No EMC rec points read " << endl ;
}
if(ppsdRecPoints ){
for(irecp = 0; irecp < ppsdRecPoints->GetEntries() ; irecp ++){
AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)ppsdRecPoints->At(irecp) ;
ppsd->GetLocalPosition(pos) ;
cout << "PPSD " << irecp << " " << ppsd->GetPHOSMod() << " " << pos.X() << " " << pos.Z() << endl ;
if(ppsd->GetPHOSMod()==Nmod){
ppsd->GetLocalPosition(pos) ;
if(ppsd->GetUp())
ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
else
ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
}
}
}
else{
cout << "No PPSD/CPV rec points read " << endl ;
}
AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
if(recParticleList ){
for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
{
recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
Int_t moduleNumberRec ;
Double_t recX, recZ ;
fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
if(moduleNumberRec == Nmod){
Double_t minDistance = 5. ;
Int_t closestPrimary = -1 ;
Int_t numberofprimaries ;
Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
Int_t index ;
TParticle * primary ;
Double_t distance = minDistance ;
for ( index = 0 ; index < numberofprimaries ; index++){
primary = gAlice->Particle(listofprimaries[index]) ;
Int_t moduleNumber ;
Double_t primX, primZ ;
fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
if(moduleNumberRec == moduleNumber)
distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
if(minDistance > distance)
{
minDistance = distance ;
closestPrimary = listofprimaries[index] ;
}
}
if(closestPrimary >=0 ){
Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
if(primaryType==22)
recPhot->Fill(recZ,recX,recParticle->Energy()) ;
else
if(primaryType==-2112)
recNbar->Fill(recZ,recX,recParticle->Energy()) ;
}
}
}
}
else{
cout << "Not Rec Prticles read " << endl ;
}
digitOccupancy->Draw("box") ;
sdigitOccupancy->SetLineColor(5) ;
sdigitOccupancy->Draw("box") ;
emcOccupancy->SetLineColor(2) ;
emcOccupancy->Draw("boxsame") ;
ppsdUp->SetLineColor(3) ;
ppsdUp->Draw("boxsame") ;
ppsdLow->SetLineColor(4) ;
ppsdLow->Draw("boxsame") ;
phot->SetLineColor(8) ;
phot->Draw("boxsame") ;
nbar->SetLineColor(6) ;
nbar->Draw("boxsame") ;
}
// //____________________________________________________________________________
// void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )
// {
// // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
// // for events from FirstEvent to Nevents
// Int_t ievent ;
// for ( ievent=firstEvent; ievent<nevents; ievent++) {
// if (ievent==firstEvent) {
// cout << "Analyze > Starting Reconstructing " << endl ;
// //========== Create the Clusterizer
// fClu = new AliPHOSClusterizerv1() ;
// //========== Creates the track segment maker
// fTrs = new AliPHOSTrackSegmentMakerv1() ;
// // fTrs->UnsetUnfoldFlag() ;
// //========== Creates the particle identifier
// fPID = new AliPHOSPIDv1() ;
// fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
// //========== Creates the Reconstructioner
// fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
// if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
// }
// if (fDebugLevel != 0 ||
// (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
// cout << "======= Analyze ======> Event " << ievent+1 << endl ;
// gAlice->GetEvent(ievent) ;
// gAlice->SetEvent(ievent) ;
// if(gAlice->TreeS() == 0) gAlice->MakeTree("S");
// fPHOS->MakeBranch("S") ;
// fPHOS->Hits2SDigits() ;
// if(gAlice->TreeD() == 0) gAlice->MakeTree("D");
// fPHOS->MakeBranch("D") ;
// fPHOS->SDigits2Digits() ;
// if(gAlice->TreeR() == 0) gAlice->MakeTree("R");
// fPHOS->Reconstruction(fRec);
// gAlice->TreeS()->Fill() ;
// gAlice->TreeS()->Write(0,TObject::kOverwrite);
// gAlice->TreeD()->Fill() ;
// gAlice->TreeD()->Write(0,TObject::kOverwrite);
// }
// if(fClu) {delete fClu ; fClu =0 ;}
// if(fPID) {delete fPID ; fPID =0 ;}
// if(fRec) {delete fRec ; fRec =0 ;}
// if(fTrs) {delete fTrs ; fTrs =0 ;}
// }
//-------------------------------------------------------------------------------------
void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast)
{
// //
// // Read and print generated and reconstructed hits in CPV
// // for events from EvFirst to Nevent.
// // If only EvFirst is defined, print only this one event.
// // Author: Yuri Kharlov
// // 12 October 2000
// //
// if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
// for ( Int_t ievent=EvFirst; ievent<=EvLast; ievent++) {
// //========== Event Number>
// cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
// //=========== Connects the various Tree's for evt
// Int_t ntracks = gAlice->GetEvent(ievent);
// //========== Creating branches ===================================
// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
// // Read and print CPV hits
// AliPHOSCPVModule cpvModule;
// TClonesArray *cpvHits;
// Int_t nCPVhits;
// AliPHOSCPVHit *cpvHit;
// TLorentzVector p;
// Float_t xgen, zgen;
// Int_t ipart;
// Int_t nGenHits = 0;
// for (Int_t itrack=0; itrack<ntracks; itrack++) {
// //=========== Get the Hits Tree for the Primary track itrack
// gAlice->ResetHits();
// gAlice->TreeH()->GetEvent(itrack);
// Int_t iModule = 0 ;
// for (iModule=0; iModule < fGeom->GetNCPVModules(); iModule++) {
// cpvModule = fPHOS->GetCPVModule(iModule);
// cpvHits = cpvModule.Hits();
// nCPVhits = cpvHits->GetEntriesFast();
// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
// nGenHits++;
// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
// p = cpvHit->GetMomentum();
// xgen = cpvHit->X();
// zgen = cpvHit->Y();
// ipart = cpvHit->GetIpart();
// printf("CPV hit in module %d: ",iModule+1);
// printf(" p = (%f, %f, %f, %f) GeV,n",
// p.Px(),p.Py(),p.Pz(),p.Energy());
// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %dn",
// xgen,zgen,ipart);
// }
// }
// }
// // Read and print CPV reconstructed points
// //=========== Gets the Reconstruction TTree
// gAlice->TreeR()->GetEvent(0) ;
// printf("Recpoints: %dn",(*fPHOS->CpvRecPoints())->GetEntries());
// TIter nextRP(*fPHOS->CpvRecPoints() ) ;
// AliPHOSCpvRecPoint *cpvRecPoint ;
// Int_t nRecPoints = 0;
// while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
// nRecPoints++;
// TVector3 locpos;
// cpvRecPoint->GetLocalPosition(locpos);
// Int_t phosModule = cpvRecPoint->GetPHOSMod();
// printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cmn",
// phosModule,locpos.X(),locpos.Z());
// }
// printf("This event has %d generated hits and %d reconstructed pointsn",
// nGenHits,nRecPoints);
// }
}
//____________________________________________________________________________
void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
{
// //
// // Analyzes CPV characteristics
// // Author: Yuri Kharlov
// // 9 October 2000
// //
// // Book histograms
// TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.);
// TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.);
// TH1F *hDr = new TH1F("hDr" ,"CPV r-resolution@reconstruction",100, 0. , 5.);
// TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5);
// TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" , 21,-0.5,20.5);
// TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" , 21,-0.5,20.5);
// cout << "Start CPV Analysis"<< endl ;
// for ( Int_t ievent=0; ievent<Nevents; ievent++) {
// //========== Event Number>
// // if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
// cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
// //=========== Connects the various Tree's for evt
// Int_t ntracks = gAlice->GetEvent(ievent);
// //========== Creating branches ===================================
// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
// // Create and fill arrays of hits for each CPV module
// Int_t nOfModules = fGeom->GetNModules();
// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
// Int_t iModule = 0;
// for (iModule=0; iModule < nOfModules; iModule++)
// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
// AliPHOSCPVModule cpvModule;
// TClonesArray *cpvHits;
// Int_t nCPVhits;
// AliPHOSCPVHit *cpvHit;
// TLorentzVector p;
// Float_t xzgen[2];
// Int_t ipart;
// // First go through all primary tracks and fill the arrays
// // of hits per each CPV module
// for (Int_t itrack=0; itrack<ntracks; itrack++) {
// // Get the Hits Tree for the Primary track itrack
// gAlice->ResetHits();
// gAlice->TreeH()->GetEvent(itrack);
// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
// cpvModule = fPHOS->GetCPVModule(iModule);
// cpvHits = cpvModule.Hits();
// nCPVhits = cpvHits->GetEntriesFast();
// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
// p = cpvHit->GetMomentum();
// xzgen[0] = cpvHit->X();
// xzgen[1] = cpvHit->Y();
// ipart = cpvHit->GetIpart();
// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
// }
// cpvModule.Clear();
// }
// }
// for (iModule=0; iModule < nOfModules; iModule++) {
// Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
// printf("Module %d has %d hitsn",iModule,nsum);
// }
// // Then go through reconstructed points and for each find
// // the closeset hit
// // The distance from the rec.point to the closest hit
// // gives the coordinate resolution of the CPV
// // Get the Reconstruction Tree
// gAlice->TreeR()->GetEvent(0) ;
// TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
// AliPHOSCpvRecPoint *cpvRecPoint ;
// Float_t xgen, zgen;
// while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
// TVector3 locpos;
// cpvRecPoint->GetLocalPosition(locpos);
// Int_t phosModule = cpvRecPoint->GetPHOSMod();
// Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity();
// Int_t rpMultX, rpMultZ;
// cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
// Float_t xrec = locpos.X();
// Float_t zrec = locpos.Z();
// Float_t dxmin = 1.e+10;
// Float_t dzmin = 1.e+10;
// Float_t r2min = 1.e+10;
// Float_t r2;
// cpvHits = hitsPerModule[phosModule-1];
// Int_t nCPVhits = cpvHits->GetEntriesFast();
// for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
// cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
// xgen = cpvHit->X();
// zgen = cpvHit->Y();
// r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
// if ( r2 < r2min ) {
// r2min = r2;
// dxmin = xgen - xrec;
// dzmin = zgen - zrec;
// }
// }
// hDx ->Fill(dxmin);
// hDz ->Fill(dzmin);
// hDr ->Fill(TMath::Sqrt(r2min));
// hNrp ->Fill(rpMult);
// hNrpX->Fill(rpMultX);
// hNrpZ->Fill(rpMultZ);
// }
// delete [] hitsPerModule;
// }
// // Save histograms
// Text_t outputname[80] ;
// sprintf(outputname,"%s.analyzed",fRootFile->GetName());
// TFile output(outputname,"RECREATE");
// output.cd();
// hDx ->Write() ;
// hDz ->Write() ;
// hDr ->Write() ;
// hNrp ->Write() ;
// hNrpX->Write() ;
// hNrpZ->Write() ;
// // Plot histograms
// TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
// gStyle->SetOptStat(111111);
// gStyle->SetOptFit(1);
// gStyle->SetOptDate(1);
// cpvCanvas->Divide(3,2);
// cpvCanvas->cd(1);
// gPad->SetFillColor(10);
// hNrp->SetFillColor(16);
// hNrp->Draw();
// cpvCanvas->cd(2);
// gPad->SetFillColor(10);
// hNrpX->SetFillColor(16);
// hNrpX->Draw();
// cpvCanvas->cd(3);
// gPad->SetFillColor(10);
// hNrpZ->SetFillColor(16);
// hNrpZ->Draw();
// cpvCanvas->cd(4);
// gPad->SetFillColor(10);
// hDx->SetFillColor(16);
// hDx->Fit("gaus");
// hDx->Draw();
// cpvCanvas->cd(5);
// gPad->SetFillColor(10);
// hDz->SetFillColor(16);
// hDz->Fit("gaus");
// hDz->Draw();
// cpvCanvas->cd(6);
// gPad->SetFillColor(10);
// hDr->SetFillColor(16);
// hDr->Draw();
// cpvCanvas->Print("CPV.ps");
}
//____________________________________________________________________________
void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
{
// Calculates Real and Mixed invariant mass distributions
const Int_t knMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/knMixedEvents) ;
//========== Booking Histograms
TH2D * hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ;
TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
TH2D * hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ;
TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
Int_t ievent;
Int_t eventInMixedLoop ;
Int_t nRecParticles[4];//knMixedEvents] ;
AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", knMixedEvents*1000) ;
for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++ ){
Int_t iRecPhot = 0 ;
for ( ievent=0; ievent < knMixedEvents; ievent++){
Int_t absEventNumber = eventInMixedLoop*knMixedEvents + ievent ;
//=========== Connects the various Tree's for evt
gAlice->GetEvent(absEventNumber);
//========== Creating branches ===================================
fPHOS->SetTreeAddress() ;
gAlice->TreeD()->GetEvent(0) ;
gAlice->TreeR()->GetEvent(0) ;
TClonesArray * recParticleList = fPHOS->RecParticles() ;
AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
{
recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
iRecPhot++;
}
}
nRecParticles[ievent] = iRecPhot-1 ;
}
//Now calculate invariant mass:
Int_t irp1,irp2 ;
Int_t nCurEvent = 0 ;
for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
Double_t invMass ;
invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
(rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
(rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
if(invMass> 0)
invMass = TMath::Sqrt(invMass);
Double_t pt ;
pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
if(irp1 > nRecParticles[nCurEvent])
nCurEvent++;
if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
hRealEM->Fill(invMass,pt);
if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
hRealPhot->Fill(invMass,pt);
}
else{
hMixedEM->Fill(invMass,pt);
if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
hMixedPhot->Fill(invMass,pt);
} //real-mixed
} //loop over second rp
}//loop over first rp
allRecParticleList->Delete() ;
} //Loop over events
delete allRecParticleList ;
//writing output
TFile output("invmass.root","RECREATE");
output.cd();
hRealEM->Write() ;
hRealPhot->Write() ;
hMixedEM->Write() ;
hMixedPhot->Write() ;
output.Write();
output.Close();
}
//____________________________________________________________________________
void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast)
{
// //
// // Read and print generated and reconstructed hits in EMC
// // for events from EvFirst to Nevent.
// // If only EvFirst is defined, print only this one event.
// // Author: Yuri Kharlov
// // 24 November 2000
// //
// if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
// Int_t ievent;
// for (ievent=EvFirst; ievent<=EvLast; ievent++) {
// //========== Event Number>
// cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
// //=========== Connects the various Tree's for evt
// Int_t ntracks = gAlice->GetEvent(ievent);
// fPHOS->SetTreeAddress() ;
// gAlice->TreeD()->GetEvent(0) ;
// gAlice->TreeR()->GetEvent(0) ;
// // Loop over reconstructed particles
// TClonesArray ** recParticleList = fPHOS->RecParticles() ;
// AliPHOSRecParticle * recParticle ;
// Int_t iRecParticle ;
// Int_t *primList;
// Int_t nPrimary;
// for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
// Float_t recE = recParticle->Energy();
// primList = recParticle->GetPrimaries(nPrimary);
// Int_t moduleNumberRec ;
// Double_t recX, recZ ;
// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
// printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %dn",
// moduleNumberRec,recX,recZ,recE,*primList);
// }
// // Read and print EMC hits from EMCn branches
// AliPHOSCPVModule emcModule;
// TClonesArray *emcHits;
// Int_t nEMChits;
// AliPHOSCPVHit *emcHit;
// TLorentzVector p;
// Float_t xgen, zgen;
// Int_t ipart, primary;
// Int_t nGenHits = 0;
// for (Int_t itrack=0; itrack<ntracks; itrack++) {
// //=========== Get the Hits Tree for the Primary track itrack
// gAlice->ResetHits();
// gAlice->TreeH()->GetEvent(itrack);
// Int_t iModule = 0 ;
// for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
// emcModule = fPHOS->GetEMCModule(iModule);
// emcHits = emcModule.Hits();
// nEMChits = emcHits->GetEntriesFast();
// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
// nGenHits++;
// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
// p = emcHit->GetMomentum();
// xgen = emcHit->X();
// zgen = emcHit->Y();
// ipart = emcHit->GetIpart();
// primary= emcHit->GetTrack();
// printf("EMC hit A: module %d, ",iModule+1);
// printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,n",
// p.Px(),p.Py(),p.Pz(),p.Energy());
// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %dn",
// xgen,zgen,ipart,primary);
// }
// }
// }
// // // Read and print EMC hits from PHOS branch
// // for (Int_t itrack=0; itrack<ntracks; itrack++) {
// // //=========== Get the Hits Tree for the Primary track itrack
// // gAlice->ResetHits();
// // gAlice->TreeH()->GetEvent(itrack);
// // TClonesArray *hits = fPHOS->Hits();
// // AliPHOSHit *hit ;
// // Int_t ihit;
// // for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
// // hit = (AliPHOSHit*)hits->At(ihit) ;
// // Float_t hitXYZ[3];
// // hitXYZ[0] = hit->X();
// // hitXYZ[1] = hit->Y();
// // hitXYZ[2] = hit->Z();
// // ipart = hit->GetPid();
// // primary = hit->GetPrimary();
// // Int_t absId = hit->GetId();
// // Int_t relId[4];
// // fGeom->AbsToRelNumbering(absId, relId) ;
// // Int_t module = relId[0];
// // if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
// // printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %dn",
// // module,hitXYZ[0],hitXYZ[2],ipart,primary);
// // }
// // }
// }
}
//____________________________________________________________________________
void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
{
// //
// // Read generated and reconstructed hits in EMC for Nevents events.
// // Plots the coordinate and energy resolution histograms.
// // Coordinate resolution is a difference between the reconstructed
// // coordinate and the exact coordinate on the face of the PHOS
// // Author: Yuri Kharlov
// // 27 November 2000
// //
// // Book histograms
// TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.);
// TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.);
// TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.);
// TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
// TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
// TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
// cout << "Start EMC Analysis"<< endl ;
// for (Int_t ievent=0; ievent<Nevents; ievent++) {
// //========== Event Number>
// if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
// cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
// //=========== Connects the various Tree's for evt
// Int_t ntracks = gAlice->GetEvent(ievent);
// fPHOS->SetTreeAddress() ;
// gAlice->TreeD()->GetEvent(0) ;
// gAlice->TreeR()->GetEvent(0) ;
// // Create and fill arrays of hits for each EMC module
// Int_t nOfModules = fGeom->GetNModules();
// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
// Int_t iModule;
// for (iModule=0; iModule < nOfModules; iModule++)
// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
// AliPHOSCPVModule emcModule;
// TClonesArray *emcHits;
// Int_t nEMChits;
// AliPHOSCPVHit *emcHit;
// // First go through all primary tracks and fill the arrays
// // of hits per each EMC module
// for (Int_t itrack=0; itrack<ntracks; itrack++) {
// // Get the Hits Tree for the Primary track itrack
// gAlice->ResetHits();
// gAlice->TreeH()->GetEvent(itrack);
// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
// emcModule = fPHOS->GetEMCModule(iModule);
// emcHits = emcModule.Hits();
// nEMChits = emcHits->GetEntriesFast();
// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
// }
// emcModule.Clear();
// }
// }
// // Loop over reconstructed particles
// TClonesArray ** recParticleList = fPHOS->RecParticles() ;
// AliPHOSRecParticle * recParticle ;
// Int_t nEMCrecs = (*recParticleList)->GetEntries();
// if (nEMCrecs == 1) {
// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
// Float_t recE = recParticle->Energy();
// Int_t phosModule;
// Double_t recX, recZ ;
// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
// // for this rec.point take the hit list in the same PHOS module
// emcHits = hitsPerModule[phosModule-1];
// Int_t nEMChits = emcHits->GetEntriesFast();
// if (nEMChits == 1) {
// Float_t genX, genZ, genE;
// for (Int_t ihit=0; ihit<nEMChits; ihit++) {
// emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
// genX = emcHit->X();
// genZ = emcHit->Y();
// genE = emcHit->GetMomentum().E();
// }
// Float_t dx = recX - genX;
// Float_t dz = recZ - genZ;
// Float_t de = recE - genE;
// hDx1 ->Fill(dx);
// hDz1 ->Fill(dz);
// hDE1 ->Fill(de);
// hDx2 ->Fill(genE,dx);
// hDz2 ->Fill(genE,dz);
// hDE2 ->Fill(genE,recE);
// }
// }
// delete [] hitsPerModule;
// }
// // Save histograms
// Text_t outputname[80] ;
// sprintf(outputname,"%s.analyzed",fRootFile->GetName());
// TFile output(outputname,"RECREATE");
// output.cd();
// hDx1 ->Write() ;
// hDz1 ->Write() ;
// hDE1 ->Write() ;
// hDx2 ->Write() ;
// hDz2 ->Write() ;
// hDE2 ->Write() ;
// // Plot histograms
// TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
// gStyle->SetOptStat(111111);
// gStyle->SetOptFit(1);
// gStyle->SetOptDate(1);
// emcCanvas->Divide(3,1);
// emcCanvas->cd(1);
// gPad->SetFillColor(10);
// hDx1->SetFillColor(16);
// hDx1->Draw();
// emcCanvas->cd(2);
// gPad->SetFillColor(10);
// hDz1->SetFillColor(16);
// hDz1->Draw();
// emcCanvas->cd(3);
// gPad->SetFillColor(10);
// hDE1->SetFillColor(16);
// hDE1->Draw();
// emcCanvas->Print("EMC.ps");
}
//____________________________________________________________________________
void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
{
// analyzes Nevents events and calculate Energy and Position resolution as well as
// probaility of correct indentifiing of the incident particle
//========== Booking Histograms
cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
BookResolutionHistograms();
Int_t counter[9][5] ;
Int_t i1,i2,totalInd = 0 ;
for(i1 = 0; i1<9; i1++)
for(i2 = 0; i2<5; i2++)
counter[i1][i2] = 0 ;
Int_t totalPrimary = 0 ;
Int_t totalRecPart = 0 ;
Int_t totalRPwithPrim = 0 ;
Int_t ievent;
cout << "Start Analysing"<< endl ;
for ( ievent=0; ievent<Nevents; ievent++)
{
//========== Event Number>
// if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
//=========== Connects the various Tree's for evt
gAlice->GetEvent(ievent);
//=========== Gets the Kine TTree
gAlice->TreeK()->GetEvent(0) ;
//=========== Gets the list of Primari Particles
TParticle * primary ;
Int_t iPrimary ;
for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
{
primary = gAlice->Particle(iPrimary) ;
Int_t primaryType = primary->GetPdgCode() ;
if( primaryType == 22 ) {
Int_t moduleNumber ;
Double_t primX, primZ ;
fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
if(moduleNumber){
fhPrimary->Fill(primary->Energy()) ;
if(primary->Energy() > 0.3)
totalPrimary++ ;
}
}
}
fPHOS->SetTreeAddress() ;
gAlice->TreeD()->GetEvent(0) ;
gAlice->TreeR()->GetEvent(0) ;
TClonesArray * recParticleList = fPHOS->RecParticles() ;
AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
{
recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
Int_t moduleNumberRec ;
Double_t recX, recZ ;
fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
Double_t minDistance = 100. ;
Int_t closestPrimary = -1 ;
Int_t numberofprimaries ;
Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
Int_t index ;
TParticle * primary ;
Double_t distance = minDistance ;
Double_t dX, dZ;
Double_t dXmin = 0.;
Double_t dZmin = 0. ;
for ( index = 0 ; index < numberofprimaries ; index++){
primary = gAlice->Particle(listofprimaries[index]) ;
Int_t moduleNumber ;
Double_t primX, primZ ;
fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
if(moduleNumberRec == moduleNumber) {
dX = recX - primX;
dZ = recZ - primZ;
distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
if(minDistance > distance) {
minDistance = distance ;
dXmin = dX;
dZmin = dZ;
closestPrimary = listofprimaries[index] ;
}
}
}
totalRecPart++ ;
if(closestPrimary >=0 ){
totalRPwithPrim++;
Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
// TParticlePDG* pDGparticle = gAlice->ParticleAt(closestPrimary)->GetPDG();
// Double_t charge = PDGparticle->Charge() ;
// if(charge)
// cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
Int_t primaryCode ;
switch(primaryType)
{
case 22:
primaryCode = 0; //Photon
fhAllEnergy ->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy()) ;
fhAllPosition ->Fill(gAlice->Particle(closestPrimary)->Energy(), minDistance) ;
fhAllPositionX->Fill(dXmin);
fhAllPositionZ->Fill(dZmin);
break;
case 11 :
primaryCode = 1; //Electron
break;
case -11 :
primaryCode = 1; //positron
break;
case 321 :
primaryCode = 4; //K+
break;
case -321 :
primaryCode = 4; //K-
break;
case 310 :
primaryCode = 4; //K0s
break;
case 130 :
primaryCode = 4; //K0l
break;
case 211 :
primaryCode = 2; //K0l
break;
case -211 :
primaryCode = 2; //K0l
break;
case 2212 :
primaryCode = 2; //K0l
break;
case -2212 :
primaryCode = 2; //K0l
break;
default:
primaryCode = 3; //ELSE
break;
}
switch(recParticle->GetType())
{
case AliPHOSFastRecParticle::kGAMMA:
if(primaryType == 22){
fhPhotEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
fhPhotPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryType == 2112){ //neutron
fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryType == -2112){ //neutron ~
fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryCode == 2){
fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
counter[0][primaryCode]++;
break;
case AliPHOSFastRecParticle::kELECTRON:
if(primaryType == 22){
fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryType == 2112){ //neutron
fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryType == -2112){ //neutron ~
fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryCode == 2){
fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
counter[1][primaryCode]++;
break;
case AliPHOSFastRecParticle::kNEUTRALHA:
if(primaryType == 22)
fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
counter[2][primaryCode]++;
break ;
case AliPHOSFastRecParticle::kNEUTRALEM:
if(primaryType == 22){
fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(),recParticle->Energy() ) ;
fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance ) ;
fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryType == 2112) //neutron
fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
if(primaryType == -2112) //neutron ~
fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
if(primaryCode == 2)
fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
counter[3][primaryCode]++;
break ;
case AliPHOSFastRecParticle::kCHARGEDHA:
if(primaryType == 22) //photon
fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
counter[4][primaryCode]++ ;
break ;
case AliPHOSFastRecParticle::kGAMMAHA:
if(primaryType == 22){ //photon
fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryType == 2112){ //neutron
fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryType == -2112){ //neutron ~
fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryCode == 2){
fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
counter[5][primaryCode]++ ;
break ;
case AliPHOSFastRecParticle::kABSURDEM:
counter[6][primaryCode]++ ;
fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
break;
case AliPHOSFastRecParticle::kABSURDHA:
counter[7][primaryCode]++ ;
break;
default:
counter[8][primaryCode]++ ;
break;
}
}
}
} // endfor
SaveHistograms();
cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
cout << "Resolutions: Total primary " << totalPrimary << endl ;
cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ;
cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ;
cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ;
cout << " Detected as photon " << counter[0][0] << " " << counter[0][1] << " " << counter[0][2] << " " <<counter[0][3] << " " << counter[0][4] << endl ;
cout << " Detected as electron " << counter[1][0] << " " << counter[1][1] << " " << counter[1][2] << " " <<counter[1][3] << " " << counter[1][4] << endl ;
cout << " Detected as neutral hadron " << counter[2][0] << " " << counter[2][1] << " " << counter[2][2] << " " <<counter[2][3] << " " << counter[2][4] << endl ;
cout << " Detected as neutral EM " << counter[3][0] << " " << counter[3][1] << " " << counter[3][2] << " " <<counter[3][3] << " " << counter[3][4] << endl ;
cout << " Detected as charged hadron " << counter[4][0] << " " << counter[4][1] << " " << counter[4][2] << " " <<counter[4][3] << " " << counter[4][4] << endl ;
cout << " Detected as gamma-hadron " << counter[5][0] << " " << counter[5][1] << " " << counter[5][2] << " " <<counter[5][3] << " " << counter[5][4] << endl ;
cout << " Detected as Absurd EM " << counter[6][0] << " " << counter[6][1] << " " << counter[6][2] << " " <<counter[6][3] << " " << counter[6][4] << endl ;
cout << " Detected as absurd hadron " << counter[7][0] << " " << counter[7][1] << " " << counter[7][2] << " " <<counter[7][3] << " " << counter[7][4] << endl ;
cout << " Detected as undefined " << counter[8][0] << " " << counter[8][1] << " " << counter[8][2] << " " <<counter[8][3] << " " << counter[8][4] << endl ;
for(i1 = 0; i1<9; i1++)
for(i2 = 0; i2<5; i2++)
totalInd+=counter[i1][i2] ;
cout << "Indentified particles " << totalInd << endl ;
} // endfunction
//____________________________________________________________________________
void AliPHOSAnalyze::BookingHistograms()
{
// Books the histograms where the results of the analysis are stored (to be changed)
delete fhEmcDigit ;
delete fhVetoDigit ;
delete fhConvertorDigit ;
delete fhEmcCluster ;
delete fhVetoCluster ;
delete fhConvertorCluster ;
delete fhConvertorEmc ;
fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.);
fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
}
//____________________________________________________________________________
void AliPHOSAnalyze::BookResolutionHistograms()
{
// Books the histograms where the results of the Resolution analysis are stored
// if(fhAllEnergy)
// delete fhAllEnergy ;
// if(fhPhotEnergy)
// delete fhPhotEnergy ;
// if(fhEMEnergy)
// delete fhEMEnergy ;
// if(fhPPSDEnergy)
// delete fhPPSDEnergy ;
fhAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
fhEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.);
fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
// if(fhAllPosition)
// delete fhAllPosition ;
// if(fhPhotPosition)
// delete fhPhotPosition ;
// if(fhEMPosition)
// delete fhEMPosition ;
// if(fhPPSDPosition)
// delete fhPPSDPosition ;
fhAllPosition = new TH2F("hAllPosition", "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
fhEMPosition = new TH2F("hEMPosition", "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.);
fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
// if(fhAllReg)
// delete fhAllReg ;
// if(fhPhotReg)
// delete fhPhotReg ;
// if(fhNReg)
// delete fhNReg ;
// if(fhNBarReg)
// delete fhNBarReg ;
// if(fhChargedReg)
// delete fhChargedReg ;
fhAllReg = new TH1F("hAllReg", "All primaries registered as photon", 100, 0., 5.);
fhPhotReg = new TH1F("hPhotReg", "Photon registered as photon", 100, 0., 5.);
fhNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.);
fhNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.);
fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
// if(fhAllEM)
// delete fhAllEM ;
// if(fhPhotEM)
// delete fhPhotEM ;
// if(fhNEM)
// delete fhNEM ;
// if(fhNBarEM)
// delete fhNBarEM ;
// if(fhChargedEM)
// delete fhChargedEM ;
fhAllEM = new TH1F("hAllEM", "All primary registered as EM",100, 0., 5.);
fhPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
fhNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
fhNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
// if(fhAllPPSD)
// delete fhAllPPSD ;
// if(fhPhotPPSD)
// delete fhPhotPPSD ;
// if(fhNPPSD)
// delete fhNPPSD ;
// if(fhNBarPPSD)
// delete fhNBarPPSD ;
// if(fhChargedPPSD)
// delete fhChargedPPSD ;
fhAllPPSD = new TH1F("hAllPPSD", "All primary registered as PPSD",100, 0., 5.);
fhPhotPPSD = new TH1F("hPhotPPSD", "Photon registered as PPSD", 100, 0., 5.);
fhNPPSD = new TH1F("hNPPSD", "N registered as PPSD", 100, 0., 5.);
fhNBarPPSD = new TH1F("hNBarPPSD", "NBar registered as PPSD", 100, 0., 5.);
fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
// if(fhPrimary)
// delete fhPrimary ;
fhPrimary= new TH1F("hPrimary", "hPrimary", 100, 0., 5.);
// if(fhAllRP)
// delete fhAllRP ;
// if(fhVeto)
// delete fhVeto ;
// if(fhShape)
// delete fhShape ;
// if(fhPPSD)
// delete fhPPSD ;
fhAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
fhVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.);
fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
fhPPSD = new TH1F("hPPSD", "All PPSD photon particles", 100, 0., 5.);
// if(fhPhotPhot)
// delete fhPhotPhot ;
// if(fhPhotElec)
// delete fhPhotElec ;
// if(fhPhotNeuH)
// delete fhPhotNeuH ;
// if(fhPhotNuEM)
// delete fhPhotNuEM ;
// if(fhPhotChHa)
// delete fhPhotChHa ;
// if(fhPhotGaHa)
// delete fhPhotGaHa ;
fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.); //Photon registered as photon
fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.); //Photon registered as Electron
fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.); //Photon registered as Neutral Hadron
fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.); //Photon registered as Neutral EM
fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.); //Photon registered as Charged Hadron
fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.); //Photon registered as Gamma-Hadron
}
//____________________________________________________________________________
Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
{
// Open the root file named "name"
fRootFile = new TFile(name, "update") ;
return fRootFile->IsOpen() ;
}
//____________________________________________________________________________
void AliPHOSAnalyze::SaveHistograms()
{
// Saves the histograms in a root file named "name.analyzed"
Text_t outputname[80] ;
sprintf(outputname,"%s.analyzed",fRootFile->GetName());
TFile output(outputname,"RECREATE");
output.cd();
if (fhAllEnergy)
fhAllEnergy->Write() ;
if (fhPhotEnergy)
fhPhotEnergy->Write() ;
if(fhEMEnergy)
fhEMEnergy->Write() ;
if(fhPPSDEnergy)
fhPPSDEnergy->Write() ;
if(fhAllPosition)
fhAllPosition->Write() ;
if(fhAllPositionX)
fhAllPositionX->Write() ;
if(fhAllPositionZ)
fhAllPositionZ->Write() ;
if(fhPhotPosition)
fhPhotPosition->Write() ;
if(fhEMPosition)
fhEMPosition->Write() ;
if(fhPPSDPosition)
fhPPSDPosition->Write() ;
if (fhAllReg)
fhAllReg->Write() ;
if (fhPhotReg)
fhPhotReg->Write() ;
if(fhNReg)
fhNReg->Write() ;
if(fhNBarReg)
fhNBarReg->Write() ;
if(fhChargedReg)
fhChargedReg->Write() ;
if (fhAllEM)
fhAllEM->Write() ;
if (fhPhotEM)
fhPhotEM->Write() ;
if(fhNEM)
fhNEM->Write() ;
if(fhNBarEM)
fhNBarEM->Write() ;
if(fhChargedEM)
fhChargedEM->Write() ;
if (fhAllPPSD)
fhAllPPSD->Write() ;
if (fhPhotPPSD)
fhPhotPPSD->Write() ;
if(fhNPPSD)
fhNPPSD->Write() ;
if(fhNBarPPSD)
fhNBarPPSD->Write() ;
if(fhChargedPPSD)
fhChargedPPSD->Write() ;
if(fhPrimary)
fhPrimary->Write() ;
if(fhAllRP)
fhAllRP->Write() ;
if(fhVeto)
fhVeto->Write() ;
if(fhShape)
fhShape->Write() ;
if(fhPPSD)
fhPPSD->Write() ;
if(fhPhotPhot)
fhPhotPhot->Write() ;
if(fhPhotElec)
fhPhotElec->Write() ;
if(fhPhotNeuH)
fhPhotNeuH->Write() ;
if(fhPhotNuEM)
fhPhotNuEM->Write() ;
if(fhPhotNuEM)
fhPhotNuEM->Write() ;
if(fhPhotChHa)
fhPhotChHa->Write() ;
if(fhPhotGaHa)
fhPhotGaHa->Write() ;
if(fhEnergyCorrelations)
fhEnergyCorrelations->Write() ;
output.Write();
output.Close();
}
//____________________________________________________________________________
Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
{
return ERecPart/0.8783 ;
}
//____________________________________________________________________________
void AliPHOSAnalyze::ResetHistograms()
{
fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2)
fhEmcDigit = 0 ; // Histo of digit energies in the Emc
fhVetoDigit = 0 ; // Histo of digit energies in the Veto
fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor
fhEmcCluster = 0 ; // Histo of Cluster energies in Emc
fhVetoCluster = 0 ; // Histo of Cluster energies in Veto
fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor
fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies
fhAllEnergy = 0 ;
fhPhotEnergy = 0 ; // Total spectrum of detected photons
fhEMEnergy = 0 ; // Spectrum of detected electrons with electron primary
fhPPSDEnergy = 0 ;
fhAllPosition = 0 ;
fhAllPositionX = 0 ;
fhAllPositionZ = 0 ;
fhPhotPosition = 0 ;
fhEMPosition = 0 ;
fhPPSDPosition = 0 ;
fhPhotReg = 0 ;
fhAllReg = 0 ;
fhNReg = 0 ;
fhNBarReg = 0 ;
fhChargedReg = 0 ;
fhPhotEM = 0 ;
fhAllEM = 0 ;
fhNEM = 0 ;
fhNBarEM = 0 ;
fhChargedEM = 0 ;
fhPhotPPSD = 0 ;
fhAllPPSD = 0 ;
fhNPPSD = 0 ;
fhNBarPPSD = 0 ;
fhChargedPPSD = 0 ;
fhPrimary = 0 ;
fhPhotPhot = 0 ;
fhPhotElec = 0 ;
fhPhotNeuH = 0 ;
fhPhotNuEM = 0 ;
fhPhotChHa = 0 ;
fhPhotGaHa = 0 ;
}
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.