- Float_t etot = 0 ;
- Int_t nTrackSegmentsInModule = 0 ;
- for(index = 0; index < nTrackSegments ; index++){
- trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
- etot+= trseg->GetEnergy() ;
- if ( trseg->GetPHOSMod() == module ) {
- nTrackSegmentsInModule++ ;
- trseg->Draw("P");
- }
- }
- Text_t text[80] ;
- sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
- TPaveText * pavetext = new TPaveText(22, 80, 83, 90);
- pavetext->AddText(text) ;
- pavetext->Draw() ;
- trackcanvas->Update() ;
- cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
-
- }
-}
-//____________________________________________________________________________
-Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
-{
- fRootFile = new TFile(name) ;
- return fRootFile->IsOpen() ;
-}
-//____________________________________________________________________________
-void AliPHOSAnalyze::SavingHistograms()
-{
- Text_t outputname[80] ;// = fRootFile->GetName();
- sprintf(outputname,"%s.analyzed",fRootFile->GetName());
- TFile output(outputname,"RECREATE");
- output.cd();
- if (fhEmcDigit )
- fhEmcDigit->Write() ;
- if (fhVetoDigit )
- fhVetoDigit->Write() ;
- if (fhConvertorDigit )
- fhConvertorDigit->Write() ;
- if (fhEmcCluster )
- fhEmcCluster->Write() ;
- if (fhVetoCluster )
- fhVetoCluster->Write() ;
- if (fhConvertorCluster )
- fhConvertorCluster->Write() ;
- if (fhConvertorEmc )
- fhConvertorEmc->Write() ;
- if (fhPhotonEnergy)
- fhPhotonEnergy->Write() ;
- if (fhPhotonPositionX)
- fhPhotonPositionX->Write() ;
- if (fhPhotonPositionY)
- fhPhotonPositionX->Write() ;
- if (fhElectronEnergy)
- fhElectronEnergy->Write() ;
- if (fhElectronPositionX)
- fhElectronPositionX->Write() ;
- if (fhElectronPositionY)
- fhElectronPositionX->Write() ;
- if (fhNeutralHadronEnergy)
- fhNeutralHadronEnergy->Write() ;
- if (fhNeutralHadronPositionX)
- fhNeutralHadronPositionX->Write() ;
- if (fhNeutralHadronPositionY)
- fhNeutralHadronPositionX->Write() ;
- if (fhNeutralEMEnergy)
- fhNeutralEMEnergy->Write() ;
- if (fhNeutralEMPositionX)
- fhNeutralEMPositionX->Write() ;
- if (fhNeutralEMPositionY)
- fhNeutralEMPositionX->Write() ;
- if (fhChargedHadronEnergy)
- fhChargedHadronEnergy->Write() ;
- if (fhChargedHadronPositionX)
- fhChargedHadronPositionX->Write() ;
- if (fhChargedHadronPositionY)
- fhChargedHadronPositionX->Write() ;
- if (fhPhotonHadronEnergy)
- fhPhotonHadronEnergy->Write() ;
- if (fhPhotonHadronPositionX)
- fhPhotonHadronPositionX->Write() ;
- if (fhPhotonHadronPositionY)
- fhPhotonHadronPositionX->Write() ;
-
- output.Write();
- output.Close();
+ Double_t distance = minDistance ;
+ Double_t dX, dZ;
+ Double_t dXmin = 0.;
+ Double_t dZmin = 0. ;
+ for ( index = 0 ; index < numberofprimaries ; index++){
+ primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ;
+ Int_t moduleNumber ;
+ Double_t primX, primZ ;
+ phosgeom->ImpactOnEmc(vtx,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] ;
+ }
+ }
+ }
+
+ //===========define the "type" of closest primary
+ if(closestPrimary >=0 ){
+ Int_t primaryCode = -1;
+ primary = fRunLoader->Stack()->Particle(closestPrimary) ;
+ Int_t primaryType = primary->GetPdgCode() ;
+ if(primaryType == 22) // photon ?
+ primaryCode = 0 ;
+ else
+ if(primaryType == 2112) // neutron
+ primaryCode = 1 ;
+ else
+ if(primaryType == -2112) // Anti neutron
+ primaryCode = 2 ;
+ else
+ if(primaryType == -2122) //Anti proton
+ primaryCode = 4 ;
+ else {
+ TParticle tempo(*primary) ;
+ if(tempo.GetPDG()->Charge())
+ primaryCode = 3 ;
+ }
+
+ //==========Now look at the type of RecParticle
+ Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
+ if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
+ hPhot->Fill(energy ) ;
+ switch(primaryCode){
+ case 0:
+ hPhotReg->Fill(energy ) ;
+ break ;
+ case 1:
+ hNReg->Fill(energy ) ;
+ break ;
+ case 2:
+ hNBarReg->Fill(energy ) ;
+ break ;
+ case 3:
+ hChargedReg->Fill(energy ) ;
+ break ;
+ case 4:
+ hPbarReg->Fill(energy ) ;
+ break ;
+ default:
+ break ;
+ }
+ }
+ if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower
+ hShape->Fill(energy ) ;
+ switch(primaryCode){
+ case 0:
+ hPhotEM->Fill(energy ) ;
+ break ;
+ case 1:
+ hNEM->Fill(energy ) ;
+ break ;
+ case 2:
+ hNBarEM->Fill(energy ) ;
+ break ;
+ case 3:
+ hChargedEM->Fill(energy ) ;
+ break ;
+ case 4:
+ hPbarEM->Fill(energy ) ;
+ break ;
+ default:
+ break ;
+ }
+ }
+
+ if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) ||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) ||
+ (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral
+ hVeto->Fill(energy ) ;
+
+ //fill number of primaries identified as ...
+ if(primaryCode >= 0) // Primary code defined
+ counter[recParticle->GetType()][primaryCode]++ ;
+
+ }
+
+ } // no closest primary found
+ }
+
+
+ //=================== SaveHistograms
+ cfile->cd() ;
+ hPrimary->Write(0,kOverwrite);
+ hAllRP->Write(0,kOverwrite);
+ hPhot->Write(0,kOverwrite);
+ hShape->Write(0,kOverwrite);
+ hVeto->Write(0,kOverwrite);
+ hPhotReg->Write(0,kOverwrite);
+ hPhotEM->Write(0,kOverwrite);
+ hNReg ->Write(0,kOverwrite);
+ hNEM ->Write(0,kOverwrite);
+ hNBarReg ->Write(0,kOverwrite);
+ hNBarEM ->Write(0,kOverwrite);
+ hChargedReg ->Write(0,kOverwrite);
+ hChargedEM ->Write(0,kOverwrite);
+ hPbarReg ->Write(0,kOverwrite);
+ hPbarEM ->Write(0,kOverwrite);
+
+ cfile->Write(0,kOverwrite);
+ cfile->Close();
+ delete cfile ;
+
+
+ //print Final Table
+ maxevent = (Int_t)AliRunLoader::Instance()->TreeE()->GetEntries() ;
+
+ TString message ;
+ message = "Resolutions: Analyzed %d event(s)\n" ;
+
+ message += " Primary: Photon Neutron Antineutron Charged hadron AntiProton\n" ;
+ message += "--------------------------------------------------------------------------------\n" ;
+ message += " kGAMMA: " ;
+ message += "%d %d %d %d %d\n" ;
+ message += " kGAMMAHA: " ;
+ message += "%d %d %d %d %d\n" ;
+ message += " kNEUTRALEM: " ;
+ message += "%d %d %d %d %d\n" ;
+ message += " kNEUTRALHA: " ;
+ message += "%d %d %d %d %d\n" ;
+ message += " kABSURDEM: ";
+ message += "%d %d %d %d %d\n" ;
+ message += " kABSURDHA: " ;
+ message += "%d %d %d %d %d\n" ;
+ message += " kELECTRON: " ;
+ message += "%d %d %d %d %d\n" ;
+ message += " kCHARGEDHA: " ;
+ message += "%d %d %d %d %d\n" ;
+
+ message += "--------------------------------------------------------------------------------" ;
+
+
+ Int_t totalInd = 0 ;
+ for(i1 = 0; i1<8; i1++)
+ for(i2 = 0; i2<5; i2++)
+ totalInd+=counter[i1][i2] ;
+ message += "Indentified particles: %d" ;
+
+ AliInfo(Form(message.Data(), maxevent,
+ counter[2][0], counter[2][1], counter[2][2], counter[2][3], counter[2][4],
+ counter[3][0], counter[3][1], counter[3][2], counter[3][3], counter[3][4],
+ counter[0][0], counter[0][1], counter[0][2], counter[0][3], counter[0][4],
+ counter[1][0], counter[1][1], counter[1][2], counter[1][3], counter[1][4],
+ counter[4][0], counter[4][1], counter[4][2], counter[4][3], counter[4][4],
+ counter[5][0], counter[5][1], counter[5][2], counter[5][3], counter[5][4],
+ counter[6][0], counter[6][1], counter[6][2], counter[6][3], counter[6][4],
+ counter[7][0], counter[7][1], counter[7][2], counter[7][3], counter[7][4],
+ totalInd )) ;
+