#include "AliPHOSRecParticle.h"
#include "AliPHOSIndexToObject.h"
#include "AliPHOSHit.h"
-#include "AliPHOSCPVHit.h"
#include "AliPHOSCpvRecPoint.h"
ClassImp(AliPHOSAnalyze)
//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.) ;
//========== Create the Clusterizer
fClu = new AliPHOSClusterizerv1() ;
-
- fClu->SetEmcEnergyThreshold(0.05) ;
- fClu->SetEmcClusteringThreshold(0.20) ;
- fClu->SetPpsdEnergyThreshold (0.0000002) ;
- fClu->SetPpsdClusteringThreshold(0.0000001) ;
- fClu->SetLocalMaxCut(0.03) ;
- fClu->SetCalibrationParameters(0., 0.00000001) ;
gAlice->GetEvent(Nevent);
- TClonesArray * primaryList = gAlice->Particles();
-
TParticle * primary ;
Int_t iPrimary ;
- for ( iPrimary = 0 ; iPrimary < primaryList->GetEntries() ; iPrimary++)
+ for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
{
- primary = (TParticle*)primaryList->At(iPrimary) ;
+ primary = gAlice->Particle(iPrimary) ;
Int_t primaryType = primary->GetPdgCode() ;
if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
Int_t moduleNumber ;
}
}
- fPHOS->SetTreeAddress() ;
+ // fPHOS->SetTreeAddress() ;
- gAlice->TreeD()->GetEvent(0) ;
- gAlice->TreeR()->GetEvent(0) ;
+ gAlice->TreeS()->GetEvent(0) ;
+
+ Int_t iSDigit ;
+ AliPHOSDigit * sdigit ;
- TObjArray ** emcRecPoints = fPHOS->EmcRecPoints() ;
- TObjArray ** ppsdRecPoints = fPHOS->PpsdRecPoints() ;
- TClonesArray ** recParticleList = fPHOS->RecParticles() ;
+ 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 = 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 ;
+ }
- Int_t iDigit ;
- AliPHOSDigit * digit ;
+ gAlice->TreeD()->GetEvent(0) ;
- 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 = 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) ;
+ 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 = 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 ;
- 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());
+ 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());
+ }
}
}
-
- for(irecp = 0; irecp < (*ppsdRecPoints)->GetEntries() ; irecp ++){
- AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)(*ppsdRecPoints)->At(irecp) ;
- 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 EMC rec points read " << endl ;
+ }
+
+ if(ppsdRecPoints ){
+ for(irecp = 0; irecp < ppsdRecPoints->GetEntries() ; irecp ++){
+ AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)ppsdRecPoints->At(irecp) ;
+ 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 ;
- 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 ;
+ if(recParticleList ){
+ for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
+ {
+ recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
- Int_t numberofprimaries ;
- Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
- Int_t index ;
- TParticle * primary ;
- Double_t distance = minDistance ;
+ Int_t moduleNumberRec ;
+ Double_t recX, recZ ;
+ fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
+ if(moduleNumberRec == Nmod){
- for ( index = 0 ; index < numberofprimaries ; index++){
- primary = (TParticle*)primaryList->At(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 ){
+ Double_t minDistance = 5. ;
+ Int_t closestPrimary = -1 ;
- Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
+ Int_t numberofprimaries ;
+ Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
+ Int_t index ;
+ TParticle * primary ;
+ Double_t distance = minDistance ;
- if(primaryType==22)
- recPhot->Fill(recZ,recX,recParticle->Energy()) ;
- else
- if(primaryType==-2112)
- recNbar->Fill(recZ,recX,recParticle->Energy()) ;
- }
+ 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) ;
cout << "Analyze > Starting Reconstructing " << endl ;
//========== Create the Clusterizer
fClu = new AliPHOSClusterizerv1() ;
- fClu->SetEmcEnergyThreshold(0.05) ;
- fClu->SetEmcClusteringThreshold(0.20) ;
- fClu->SetLocalMaxCut(0.03) ;
- if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
- fClu->SetPpsdEnergyThreshold (0.0000002) ;
- fClu->SetPpsdClusteringThreshold(0.0000001) ;
- }
- else if (strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
- fClu->SetLocalMaxCutCPV(0.03) ;
- fClu->SetLogWeightCutCPV(4.0) ;
- fClu->SetCpvEnergyThreshold(0.09) ;
- }
- fClu->SetCalibrationParameters(0., 0.00000001) ;
//========== Creates the track segment maker
fTrs = new AliPHOSTrackSegmentMakerv1() ;
- // fTrs->UnsetUnfoldFlag() ;
+ // fTrs->UnsetUnfoldFlag() ;
- //========== Creates the particle identifier for GPS2 only
- if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
- fPID = new AliPHOSPIDv1() ;
- fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
- }
+ //========== 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 ;
- //=========== Connects the various Tree's for evt
- Int_t tracks = gAlice->GetEvent(ievent);
-
- fPHOS->Hit2Digit(tracks) ;
-
+ fPHOS->Enable() ;
+
+ gAlice->Hits2Digits() ;
+
//=========== Do the reconstruction
fPHOS->Reconstruction(fRec);
-
+
}
if(fClu) {delete fClu ; fClu =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++) {
+// //
+// // 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 ;
+// //========== Event Number>
+// cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
- //=========== Connects the various Tree's for evt
- Int_t ntracks = gAlice->GetEvent(ievent);
+// //=========== 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 ) ;
+// //========== Creating branches ===================================
+// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
+// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
- AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
- gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
+// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
+// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
- // Read and print CPV hits
+// // 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->GetNModules(); 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 = %d\n",
- xgen,zgen,ipart);
- }
- }
- }
-
- // Read and print CPV reconstructed points
-
- //=========== Gets the Reconstruction TTree
- gAlice->TreeR()->GetEvent(0) ;
- TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
- AliPHOSPpsdRecPoint *cpvRecPoint ;
- Int_t nRecPoints = 0;
- while( ( cpvRecPoint = (AliPHOSPpsdRecPoint *)nextRP() ) ) {
- nRecPoints++;
- TVector3 locpos;
- cpvRecPoint->GetLocalPosition(locpos);
- Int_t phosModule = cpvRecPoint->GetPHOSMod();
- printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cm\n",
- phosModule,locpos.X(),locpos.Z());
- }
- printf("This event has %d generated hits and %d reconstructed points\n",
- nGenHits,nRecPoints);
- }
+// 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 = %d\n",
+// xgen,zgen,ipart);
+// }
+// }
+// }
+
+// // Read and print CPV reconstructed points
+
+// //=========== Gets the Reconstruction TTree
+// gAlice->TreeR()->GetEvent(0) ;
+// printf("Recpoints: %d\n",(*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) cm\n",
+// phosModule,locpos.X(),locpos.Z());
+// }
+// printf("This event has %d generated hits and %d reconstructed points\n",
+// 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++) {
+// //
+// // 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 ;
+// //========== 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);
+// //=========== 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 ) ;
+// //========== Creating branches ===================================
+// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
+// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
- AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
- gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
+// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
+// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
- // Create and fill arrays of hits for each CPV module
+// // 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 hits\n",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");
+// 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 hits\n",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");
}
{
// Calculates Real and Mixed invariant mass distributions
- Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
- Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/nMixedEvents) ;
+ 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.) ;
Int_t ievent;
Int_t eventInMixedLoop ;
- Int_t nRecParticles[nMixedEvents] ;
+ Int_t nRecParticles[4];//knMixedEvents] ;
- AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", nMixedEvents*1000) ;
+ AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", knMixedEvents*1000) ;
for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++ ){
Int_t iRecPhot = 0 ;
- for ( ievent=0; ievent < nMixedEvents; ievent++){
+ for ( ievent=0; ievent < knMixedEvents; ievent++){
- Int_t absEventNumber = eventInMixedLoop*nMixedEvents + ievent ;
+ Int_t absEventNumber = eventInMixedLoop*knMixedEvents + ievent ;
//=========== Connects the various Tree's for evt
gAlice->GetEvent(absEventNumber);
gAlice->TreeD()->GetEvent(0) ;
gAlice->TreeR()->GetEvent(0) ;
- TClonesArray ** recParticleList = fPHOS->RecParticles() ;
+ TClonesArray * recParticleList = fPHOS->RecParticles() ;
AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
- for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
+ for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
{
- recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
+ recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
}
//____________________________________________________________________________
- void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst=0, Int_t EvLast=0)
+ 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++) {
+// //
+// // 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 ;
+// //========== 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() ;
+// //=========== Connects the various Tree's for evt
+// Int_t ntracks = gAlice->GetEvent(ievent);
+// fPHOS->SetTreeAddress() ;
- gAlice->TreeD()->GetEvent(0) ;
- gAlice->TreeR()->GetEvent(0) ;
+// gAlice->TreeD()->GetEvent(0) ;
+// gAlice->TreeR()->GetEvent(0) ;
- // Loop over reconstructed particles
+// // 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 = %d\n",
- moduleNumberRec,recX,recZ,recE,*primList);
- }
-
- // Read and print EMC hits from EMCn branches
+// 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 = %d\n",
+// 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 = %d\n",
- 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 = %d\n",
-// module,hitXYZ[0],hitXYZ[2],ipart,primary);
-// }
-// }
-
- }
+// 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 = %d\n",
+// 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 = %d\n",
+// // 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++) {
+// //
+// // 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 ;
+// //========== 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);
+// //=========== Connects the various Tree's for evt
+// Int_t ntracks = gAlice->GetEvent(ievent);
- fPHOS->SetTreeAddress() ;
+// fPHOS->SetTreeAddress() ;
- gAlice->TreeD()->GetEvent(0) ;
- gAlice->TreeR()->GetEvent(0) ;
+// gAlice->TreeD()->GetEvent(0) ;
+// gAlice->TreeR()->GetEvent(0) ;
- // Create and fill arrays of hits for each EMC module
+// // 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
+// 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");
+// 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");
}
gAlice->TreeK()->GetEvent(0) ;
//=========== Gets the list of Primari Particles
- TClonesArray * primaryList = gAlice->Particles();
TParticle * primary ;
Int_t iPrimary ;
- for ( iPrimary = 0 ; iPrimary < primaryList->GetEntries() ; iPrimary++)
+ for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
{
- primary = (TParticle*)primaryList->UncheckedAt(iPrimary) ;
+ primary = gAlice->Particle(iPrimary) ;
Int_t primaryType = primary->GetPdgCode() ;
if( primaryType == 22 ) {
Int_t moduleNumber ;
gAlice->TreeD()->GetEvent(0) ;
gAlice->TreeR()->GetEvent(0) ;
- TClonesArray ** recParticleList = fPHOS->RecParticles() ;
+ TClonesArray * recParticleList = fPHOS->RecParticles() ;
AliPHOSRecParticle * recParticle ;
Int_t iRecParticle ;
- for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ )
+ for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
{
- recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
+ recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
Int_t moduleNumberRec ;
Double_t dXmin = 0.;
Double_t dZmin = 0. ;
for ( index = 0 ; index < numberofprimaries ; index++){
- primary = (TParticle*)primaryList->UncheckedAt(listofprimaries[index]) ;
+ primary = gAlice->Particle(listofprimaries[index]) ;
Int_t moduleNumber ;
Double_t primX, primZ ;
fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
if(closestPrimary >=0 ){
totalRPwithPrim++;
- Int_t primaryType = ((TParticle *)primaryList->At(closestPrimary))->GetPdgCode() ;
-// TParticlePDG* pDGparticle = ((TParticle *)primaryList->At(closestPrimary))->GetPDG();
+ 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 ;
{
case 22:
primaryCode = 0; //Photon
- fhAllEnergy ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy()) ;
- fhAllPosition ->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), minDistance) ;
+ fhAllEnergy ->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy()) ;
+ fhAllPosition ->Fill(gAlice->Particle(closestPrimary)->Energy(), minDistance) ;
fhAllPositionX->Fill(dXmin);
fhAllPositionZ->Fill(dZmin);
break;
{
case AliPHOSFastRecParticle::kGAMMA:
if(primaryType == 22){
- fhPhotEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
- fhEMEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
- fhPPSDEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
+ 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(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
- fhEMPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
- fhPPSDPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
+ 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()) ) ;
case AliPHOSFastRecParticle::kELECTRON:
if(primaryType == 22){
fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
- fhEMEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
- fhEMPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
+ 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()) ) ;
}
break ;
case AliPHOSFastRecParticle::kNEUTRALEM:
if(primaryType == 22){
- fhEMEnergy->Fill(((TParticle *)primaryList->At(closestPrimary))->Energy(),recParticle->Energy() ) ;
- fhEMPosition->Fill(((TParticle *)primaryList->At(closestPrimary))->Energy(),minDistance ) ;
+ 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()) ) ;
case AliPHOSFastRecParticle::kGAMMAHA:
if(primaryType == 22){ //photon
fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
- fhPPSDEnergy->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(), recParticle->Energy() ) ;
- fhPPSDPosition->Fill(((TParticle *) primaryList->At(closestPrimary))->Energy(),minDistance) ;
+ fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
+ fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
}
if(primaryType == 2112){ //neutron
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff