#include "TChain.h" #include "TTree.h" #include "TObjArray.h" #include "TF1.h" #include "TFile.h" #include "TH1F.h" #include "TH2F.h" #include "TH2I.h" #include "TH3F.h" #include "TParticle.h" #include "TCanvas.h" #include "TStyle.h" #include "TRandom.h" #include "TROOT.h" #include "THashList.h" #include "AliAnalysisManager.h" #include "AliMCEventHandler.h" #include "AliMCEvent.h" #include "AliStack.h" #include "AliAnalysisTaskSE.h" #include "AliPHOSHijingEfficiency.h" #include "AliCaloPhoton.h" #include "AliPHOSGeometry.h" #include "AliPHOSEsdCluster.h" #include "AliPHOSCalibData.h" #include "AliESDEvent.h" #include "AliESDCaloCells.h" #include "AliESDVertex.h" #include "AliESDtrackCuts.h" #include "AliLog.h" #include "AliPID.h" #include "AliCDBManager.h" #include "AliCentrality.h" #include "AliESDtrackCuts.h" #include "AliEventplane.h" #include "TProfile.h" #include "AliOADBContainer.h" // Calculates pi0 and photon efficiencies using Hijing events // Also calculates pi0 contamination: non-vertex pi0s // And photon contamination: // - off-vertex photons // - misidentified photons // Authors: Dmitri Peressounko // Date : 05.02.2013 ClassImp(AliPHOSHijingEfficiency) //________________________________________________________________________ Double_t rnlin(Double_t x) { //This is standard PHOS non-linearity used in reconstruction //and re-calibration by 2.2% return 1.022*(0.0241+1.0504*x+0.000249*x*x) ; } //________________________________________________________________________ AliPHOSHijingEfficiency::AliPHOSHijingEfficiency(const char *name) : AliAnalysisTaskSE(name), fStack(0x0), fOutputContainer(0x0), fPHOSEvent(0x0), fPHOSCalibData(0x0), fRP(0.),fRPV0A(0.),fRPV0C(0.),fHaveTPCRP(0), fRunNumber(0), fCentrality(0.), fCenBin(0), fPHOSGeo(0), fEventCounter(0) { // Constructor for(Int_t i=0;i<1;i++){ for(Int_t j=0;j<10;j++) for(Int_t k=0;k<11;k++) fPHOSEvents[i][j][k]=0 ; } // Output slots #0 write into a TH1 container DefineOutput(1,TList::Class()); // Set bad channel map for(Int_t i=0; i<6; i++){ fPHOSBadMap[i]=new TH2I(Form("PHOS_BadMap_mod%d",i),"Bad Modules map",64,0.,64.,56,0.,56.) ; } } //________________________________________________________________________ void AliPHOSHijingEfficiency::UserCreateOutputObjects() { // Create histograms // Called once const Int_t nRuns=200 ; // ESD histograms if(fOutputContainer != NULL){ delete fOutputContainer; } fOutputContainer = new THashList(); fOutputContainer->SetOwner(kTRUE); PostData(1, fOutputContainer); //========QA histograms======= //Event selection fOutputContainer->Add(new TH2F("hSelEvents","Event selection", 10,0.,10.,nRuns,0.,float(nRuns))) ; fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", 10,0.,10.)) ; //vertex distribution fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ; //Centrality fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ; fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ; fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ; fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ; fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,20.,100,0.,100.)) ; fOutputContainer->Add(new TH2F("hCluEvsCluM","ClusterMult vs E",200,0.,20.,100,0.,20.)) ; fOutputContainer->Add(new TH3F("hCPVr","CPV radius",100,0.,20.,100,0.,20.,10,0.,100.)); //Bad Map fOutputContainer->Add(new TH2F("hCluLowM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluLowM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluLowM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluHighM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluHighM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluHighM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluVetoM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluVetoM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluVetoM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluDispM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluDispM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5)); fOutputContainer->Add(new TH2F("hCluDispM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5)); //Single photon and pi0 spectrum Int_t nPtPhot = 300 ; Double_t ptPhotMax = 30 ; Int_t nM = 500; Double_t mMin = 0.0; Double_t mMax = 1.0; //PHOS calibration QA fOutputContainer->Add(new TH2F("hPi0M11","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F("hPi0M12","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F("hPi0M13","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F("hPi0M22","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F("hPi0M23","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F("hPi0M33","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); char key[55] ; for(Int_t cent=0; cent<6; cent++){ fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad2_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad4_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad6_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotAll_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotAllcore_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotAllwou_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotDisp_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotDisp2_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotDispcore_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotDisp2core_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotDispwou_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotCPV_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotCPVcore_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotCPV2_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotBoth_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotBothcore_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotBoth2_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH1F(Form("hPhotBoth2core_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0All_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Allcore_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Allwou_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Disp_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Disp2_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Dispcore_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Disp2core_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Dispwou_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0CPV_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0CPVcore_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0CPV2_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Both_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Bothcore_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Both2_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); fOutputContainer->Add(new TH2F(Form("hPi0Both2core_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hPi0All_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hPi0Disp_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hPi0CPV_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hPi0CPV2_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hPi0Both_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hPi0Both2_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleAll_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleAllcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleAllwou_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleDisp_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleDisp2_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleDispcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleDisp2core_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleDispwou_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleCPV_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleCPVcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleCPV2_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleBoth_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleBothcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleBoth2_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hSingleBoth2core_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0All_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Allcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Allwou_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Disp_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Disp2_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Dispwou_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Dispcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Disp2core_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0CPV_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0CPVcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0CPV2_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Both_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Bothcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Both2_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Both2core_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0All_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Disp_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0CPV_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0CPV2_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Both_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiPi0Both2_a07_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleAll_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleAllwou_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleAllcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleDisp_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleDisp2_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleDispwou_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleDispcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleDisp2core_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleCPV_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleCPVcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleCPV2_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleBoth_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleBothcore_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleBoth2_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); snprintf(key,55,"hMiSingleBoth2core_cen%d",cent) ; fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax)); } //MC for(Int_t cent=0; cent<6; cent++){ snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ; snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ; snprintf(key,55,"hMC_rap_eta_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Rapidity eta",200,-1.,1.)) ; snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ; snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ; snprintf(key,55,"hMC_phi_eta_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ; snprintf(key,55,"hMC_all_gamma_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ; snprintf(key,55,"hMC_all_pi0_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ; snprintf(key,55,"hMC_all_eta_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ; snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ; snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ; snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ; fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ; } fOutputContainer->Add(new TH2F("hMC_gamma_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ; fOutputContainer->Add(new TH2F("hMC_pi0_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ; fOutputContainer->Add(new TH2F("hMC_eta_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ; Int_t nPt = 200; Double_t ptMin = 0; Double_t ptMax = 20; fOutputContainer->Add(new TH2F("Vertex","Pi0 creation vertex",nPt,ptMin,ptMax,5000,0.,500.)); fOutputContainer->Add(new TH3F("hSecondPi0RphiZ","Secondary pi0 vertex",450,0.,450.,100,0.,TMath::TwoPi(),200,-100.,100.)); fOutputContainer->Add(new TH2F("hSecondPi0RE","Secondary pi0 vertex",450,0.,450.,200,0.,20.)); fOutputContainer->Add(new TH3F("hMass_R","Mass vs radius any parent",50,0.,0.25,100,0.,10.,300,0.,600.)); fOutputContainer->Add(new TH3F("Real_pi_R","All clusters",50,0.,0.25,100,0.,10.,250,0.,500.)); fOutputContainer->Add(new TH3F("Real_pi_Z","All clusters",50,0.,0.25,100,0.,10.,100,-100.,100.)); // fOutputContainer->Add(new TH2F(Form("Real_npnp_RZ"),"All clusters",250,0.,500.,100,-100.,100.)); // fOutputContainer->Add(new TH3F(Form("Real_mass_R"),"All clusters",50,0.,0.25,100,0.,10.,300,0.,600.)); for(Int_t cen=0; cen<6; cen++){ fOutputContainer->Add(new TH1F(Form("hPrimPhot_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimEl_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimPi0_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimEta_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimPipm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimP_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimPbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimN_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimNbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimK0S_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimK0L_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimKpm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); fOutputContainer->Add(new TH1F(Form("hPrimOther_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax)); //pairs from common parents fOutputContainer->Add(new TH2F(Form("hParentAll_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentK0s_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentGamma_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentEl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentOther_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentDirPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); //common parent - pi0 fOutputContainer->Add(new TH2F(Form("hParentPi0NoPrim_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentPi0Eta_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentPi0Omega_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentPi0Pipm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentPi0Kpm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentPi0Ks_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentPi0Kl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentPi0pn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); fOutputContainer->Add(new TH2F(Form("hParentPi0antipn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax)); } //Photon contaminations fOutputContainer->Add(new TH2F("hPipmGammaConvR","Conversion radius" ,200,0.,20.,1000,0.,500.)); fOutputContainer->Add(new TH2F("hPipmElConvR","Conversion radius" ,200,0.,20.,1000,0.,500.)); fOutputContainer->Add(new TH2F("hPipmNConvR","Conversion radius" ,200,0.,20.,1000,0.,500.)); fOutputContainer->Add(new TH2F("hPipmOtherConvR","Conversion radius" ,200,0.,20.,1000,0.,500.)); fOutputContainer->Add(new TH2F("hPipmGammaConvRZ","Conversion radius" ,400,-200.,200.,1000,0.,500.)); const Int_t nTypes=24 ; char partTypes[nTypes][55] ; snprintf(partTypes[0],55,"hGammaNoPrim") ; // snprintf(partTypes[1],55,"hGammaPhot") ; // snprintf(partTypes[2],55,"hGammaEl") ; // snprintf(partTypes[3],55,"hGammaPi0") ; // snprintf(partTypes[4],55,"hGammaEta") ; // snprintf(partTypes[5],55,"hhGammaOmega") ; // snprintf(partTypes[6],55,"hGammaPipm") ; // snprintf(partTypes[7],55,"hGammaP") ; // snprintf(partTypes[8],55,"hGammaPbar") ; // snprintf(partTypes[9],55,"hGammaN") ; // snprintf(partTypes[10],55,"hGammaNbar") ; // snprintf(partTypes[11],55,"hGammaK0S") ; // snprintf(partTypes[12],55,"hGammaK0L") ; // snprintf(partTypes[13],55,"hGammaKpm") ; // snprintf(partTypes[14],55,"hGammaKstar") ; // snprintf(partTypes[15],55,"hGammaDelta") ; // snprintf(partTypes[16],55,"hGammaOtherCharged") ; // snprintf(partTypes[17],55,"hGammaOtherNeutral") ; // snprintf(partTypes[18],55,"hGammaPipmGamma") ; // snprintf(partTypes[19],55,"hGammaPipmEl") ; // snprintf(partTypes[20],55,"hGammaPipmOther") ; // snprintf(partTypes[21],55,"hGammaPipmDirect") ; // snprintf(partTypes[22],55,"hGammaPipmp") ; // snprintf(partTypes[23],55,"hGammaPipmn") ; // const Int_t nPID=12 ; char cPID[12][25] ; snprintf(cPID[0],25,"All") ; snprintf(cPID[1],25,"Allcore") ; snprintf(cPID[2],25,"CPV") ; snprintf(cPID[3],25,"CPVcore") ; snprintf(cPID[4],25,"CPV2") ; snprintf(cPID[5],25,"CPV2core") ; snprintf(cPID[6],25,"Disp") ; snprintf(cPID[7],25,"Dispcore") ; snprintf(cPID[8],25,"Disp2") ; snprintf(cPID[9],25,"Disp2core") ; snprintf(cPID[10],25,"Both") ; snprintf(cPID[11],25,"Bothcore") ; for(Int_t itype=0; itypeAdd(new TH1F(Form("%s_%s_cen%d",partTypes[itype],cPID[iPID],cen),"Cluster parents",nPt,ptMin,ptMax)); } } } PostData(1, fOutputContainer); printf("Init done \n") ; } //________________________________________________________________________ void AliPHOSHijingEfficiency::UserExec(Option_t *) { // Main loop, called for each event // Analyze ESD/AOD fStack=0 ; if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){ if(static_cast(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()) fStack = static_cast(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack(); } AliESDEvent *event = dynamic_cast(InputEvent()); if (!event) { Printf("ERROR: Could not retrieve event"); PostData(1, fOutputContainer); return; } fRunNumber=ConvertRunNumber(event->GetRunNumber()) ; FillHistogram("hSelEvents",0.5,fRunNumber-0.5) ; FillHistogram("hTotSelEvents",0.5) ; // Get PHOS rotation matrices from ESD and set them to the PHOS geometry if(fEventCounter == 0) { // Initialize the PHOS geometry fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ; //We have to apply re-calibration for pass1 LCH10h // Initialize decalibration factors in the form of the OCDB object const Float_t wDecalib = 0.055; fPHOSCalibData = new AliPHOSCalibData(); fPHOSCalibData->SetName("PhosCalibData"); for(Int_t module=1; module<=5; module++) { for(Int_t column=1; column<=56; column++) { for(Int_t row=1; row<=64; row++) { Float_t adcChannelEmc = gRandom->Gaus(1.,wDecalib); fPHOSCalibData->SetADCchannelEmc(module,column,row,adcChannelEmc); } } } for(Int_t mod=0; mod<5; mod++) { if(!event->GetPHOSMatrix(mod)) continue; fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ; Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod); } fEventCounter++ ; } // Checks if we have a primary vertex // Get primary vertices form ESD const AliESDVertex *esdVertex5 = event->GetPrimaryVertex(); Double_t vtx0[3] = {0,0,0}; // don't rely on ESD vertex, assume (0,0,0) Double_t vtx5[3] ={0.,0.,0.}; vtx5[0] = esdVertex5->GetX(); vtx5[1] = esdVertex5->GetY(); vtx5[2] = esdVertex5->GetZ(); FillHistogram("hZvertex",esdVertex5->GetZ(),fRunNumber-0.5); if (TMath::Abs(esdVertex5->GetZ()) > 10. ){ PostData(1, fOutputContainer); return; } FillHistogram("hSelEvents",1.5,fRunNumber-0.5) ; FillHistogram("hTotSelEvents",1.5) ; /* if(event->IsPileupFromSPD()){ PostData(1, fOutputContainer); return; } */ FillHistogram("hSelEvents",2.5,fRunNumber-0.5) ; FillHistogram("hTotSelEvents",2.5) ; /* //Vtx class z-bin Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ; if(zvtx<0)zvtx=0 ; if(zvtx>9)zvtx=9 ; */ //No dependence on zVtx observed, save memory Int_t zvtx=0 ; AliCentrality *centrality = event->GetCentrality(); fCentrality=centrality->GetCentralityPercentile("V0M"); if( fCentrality <= 0. || fCentrality>80. ){ PostData(1, fOutputContainer); return; } FillHistogram("hSelEvents",3.5,fRunNumber-0.5) ; FillHistogram("hTotSelEvents",3.5) ; if(fCentrality<5.) fCenBin=0 ; else if(fCentrality<10.) fCenBin=1 ; else if(fCentrality<20.) fCenBin=2 ; else if(fCentrality<40.) fCenBin=3 ; else if(fCentrality<60.) fCenBin=4 ; else if(fCentrality<80.) fCenBin=5 ; else{ PostData(1, fOutputContainer); return; } fRP=0. ; FillHistogram("hSelEvents",4.5,fRunNumber-0.5) ; FillHistogram("hTotSelEvents",4.5) ; //All event selections done FillHistogram("hCentrality",fCentrality,fRunNumber-0.5) ; //Reaction plane is defined in the range (0;pi) //We have 10 bins Double_t averageRP = 0.;//fRPV0A+fRPV0C+fRP ; Int_t irp=Int_t(10.*(averageRP)/TMath::Pi()); if(irp>9)irp=9 ; if(!fPHOSEvents[zvtx][fCenBin][irp]) fPHOSEvents[zvtx][fCenBin][irp]=new TList() ; TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ; ProcessMC() ; if(fPHOSEvent) fPHOSEvent->Clear() ; else fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ; //For re-calibration const Double_t logWeight=4.5 ; TVector3 vertex(vtx5); Int_t multClust = event->GetNumberOfCaloClusters(); Int_t inPHOS=0,inMod1=0, inMod2=0, inMod3=0 ; Double_t avrgEm1=0,avrgEm2=0,avrgEm3=0; //average cluster energy AliESDCaloCells * cells = event->GetPHOSCells() ; FillHistogram("hCenPHOSCells",fCentrality,cells->GetNumberOfCells()) ; FillHistogram("hCenTrack",fCentrality,event->GetNumberOfTracks()) ; TVector3 localPos ; for (Int_t i=0; iGetCaloCluster(i); if ( !clu->IsPHOS() || clu->E()<0.3) continue; Float_t position[3]; clu->GetPosition(position); TVector3 global(position) ; Int_t relId[4] ; fPHOSGeo->GlobalPos2RelId(global,relId) ; Int_t mod = relId[0] ; Int_t cellX = relId[2]; Int_t cellZ = relId[3] ; if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) continue ; FillHistogram("hCluEvsClu",clu->E(),clu->GetNCells()) ; FillHistogram("hCluEvsCluM",clu->E(),clu->GetM02()) ; //Apply re-Calibreation AliPHOSEsdCluster cluPHOS1(*clu); cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies Reclusterize(&cluPHOS1) ; cluPHOS1.EvalAll(logWeight,vertex); // recalculate the cluster parameters cluPHOS1.SetE(rnlin(cluPHOS1.E()));// Users's nonlinearity if(cluPHOS1.E()<0.3) continue; if(clu->GetNCells()<3) continue ; if(clu->GetM02()<0.2) continue ; if(mod==1){ avrgEm1+=cluPHOS1.E() ; inMod1++ ; } if(mod==2){ avrgEm2+=cluPHOS1.E() ; inMod2++ ; } if(mod==3){ avrgEm3+=cluPHOS1.E() ; inMod3++ ; } TLorentzVector pv1 ; cluPHOS1.GetMomentum(pv1 ,vtx0); Double_t ecore=CoreEnergy(&cluPHOS1) ; FillHistogram(Form("hCluLowM%d",mod),cellX,cellZ,1.); if(cluPHOS1.E()>1.5){ FillHistogram(Form("hCluHighM%d",mod),cellX,cellZ,1.); } if(inPHOS>=fPHOSEvent->GetSize()){ fPHOSEvent->Expand(inPHOS+50) ; } new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ; AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ; ph->SetModule(mod) ; pv1*= ecore/pv1.E() ; ph->SetMomV2(&pv1) ; ph->SetNCells(clu->GetNCells()); Double_t m02=0.,m20=0.; EvalLambdas(&cluPHOS1,0,m02, m20); ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ; EvalLambdas(&cluPHOS1,1,m02, m20); ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ; Double_t distBC=clu->GetDistanceToBadChannel(); if(distBC>2.){ FillHistogram(Form("hPhotAll_DistBad2_cen%d",fCenBin),ph->Pt()) ; if(distBC>4.){ FillHistogram(Form("hPhotAll_DistBad4_cen%d",fCenBin),ph->Pt()) ; if(distBC>6.) FillHistogram(Form("hPhotAll_DistBad6_cen%d",fCenBin),ph->Pt()) ; } } if(ph->IsDispOK()){ FillHistogram(Form("hCluDispM%d",mod),cellX,cellZ,1.); } //Track matching Double_t dx=clu->GetTrackDx() ; Double_t dz=clu->GetTrackDz() ; Bool_t cpvBit=kTRUE ; //No track matched by default Bool_t cpvBit2=kTRUE ; //More Strict criterion TArrayI * itracks = clu->GetTracksMatched() ; if(itracks->GetSize()>0){ Int_t iTr = itracks->At(0); if(iTr>=0 && iTrGetNumberOfTracks()){ AliESDtrack *track = event->GetTrack(iTr) ; Double_t pt = track->Pt() ; Short_t charge = track->Charge() ; Double_t r=TestCPV(dx, dz, pt,charge) ; cpvBit=(r>2.) ; cpvBit2=(r>4.) ; FillHistogram("hCPVr",r,pv1.E(),fCentrality); } } ph->SetCPVBit(cpvBit) ; ph->SetCPV2Bit(cpvBit2) ; if(cpvBit){ FillHistogram(Form("hCluVetoM%d",mod),cellX,cellZ,1.); } ph->SetEMCx(float(cellX)) ; ph->SetEMCz(float(cellZ)) ; ph->SetLambdas(cluPHOS1.GetM20(),cluPHOS1.GetM02()) ; ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolded Bool_t sure= kTRUE; Int_t primary=FindPrimary(clu,sure) ; //номер праймари частицы в стеке ph->SetPrimary(primary) ; ph->SetWeight(PrimaryWeight(primary)) ; inPHOS++ ; } FillHistogram("hCenPHOS",fCentrality,inPHOS) ; if(inPHOS==0){ PostData(1, fOutputContainer); fEventCounter++; return ; } for (Int_t i1=0; i1At(i1) ; Double_t w=ph1->GetWeight() ; Double_t pt = ph1->Pt() ; Double_t ptv =ph1->GetMomV2()->Pt() ; FillHistogram(Form("hPhotAll_cen%d",fCenBin),pt,w) ; FillHistogram(Form("hPhotAllcore_cen%d",fCenBin),ptv,w) ; if(ph1->IsntUnfolded()){ FillHistogram(Form("hPhotAllwou_cen%d",fCenBin),pt,w) ; } if(ph1->IsCPVOK()){ FillHistogram(Form("hPhotCPV_cen%d",fCenBin),pt,w) ; FillHistogram(Form("hPhotCPVcore_cen%d",fCenBin),ptv,w) ; } if(ph1->IsCPV2OK()){ FillHistogram(Form("hPhotCPV2_cen%d",fCenBin),pt,w) ; } if(ph1->IsDispOK()){ FillHistogram(Form("hPhotDisp_cen%d",fCenBin),pt,w) ; FillHistogram(Form("hPhotDispcore_cen%d",fCenBin),ptv,w) ; if(ph1->IsntUnfolded()){ FillHistogram(Form("hPhotDispwou_cen%d",fCenBin),pt,w) ; } if(ph1->IsCPVOK()){ FillHistogram(Form("hPhotBoth_cen%d",fCenBin),pt,w) ; FillHistogram(Form("hPhotBothcore_cen%d",fCenBin),ptv,w) ; } } if(ph1->IsDisp2OK()){ FillHistogram(Form("hPhotDisp2_cen%d",fCenBin),pt,w) ; FillHistogram(Form("hPhotDisp2core_cen%d",fCenBin),ptv,w) ; if(ph1->IsCPVOK()){ FillHistogram(Form("hPhotBoth2_cen%d",fCenBin),pt,w) ; FillHistogram(Form("hPhotBoth2core_cen%d",fCenBin),ptv,w) ; } } } const Double_t alphaCut=0.7 ; //pi0 for (Int_t i1=0; i1At(i1) ; Double_t w1=ph1->GetWeight() ; for (Int_t i2=i1+1; i2At(i2) ; Double_t w2=ph2->GetWeight() ; Double_t w = w1*w2 ; if(FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary())!=-1) w=w1 ; TLorentzVector p12 = *ph1 + *ph2; TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2()); Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ; Double_t m = p12.M() ; Double_t pt = p12.Pt() ; Double_t mv = p12.M() ; Double_t ptv = p12.Pt() ; FillHistogram(Form("hPi0All_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hPi0Allcore_cen%d",fCenBin),mv,ptv,w) ; if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){ FillHistogram(Form("hPi0Allwou_cen%d",fCenBin),m,pt,w) ; } Double_t pt1 = ph1->Pt() ; Double_t ptv1 = ph1->GetMomV2()->Pt() ; Double_t pt2 = ph2->Pt() ; Double_t ptv2 = ph2->GetMomV2()->Pt() ; FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv1,w) ; FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv2,w) ; if(ph1->IsntUnfolded()) FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt1,w) ; if(ph2->IsntUnfolded()) FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt2,w) ; if(ph1->IsCPVOK()){ FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsCPVOK()){ FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv2,w) ; } if(ph1->IsCPV2OK()){ FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt1,w) ; } if(ph2->IsCPV2OK()){ FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt2,w) ; } if(ph1->IsDispOK()){ FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt1,w) ; if(ph1->IsntUnfolded()){ FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt1,w) ; } FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsDispOK()){ FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt2,w) ; if(ph1->IsntUnfolded()){ FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt2,w) ; } FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv2,w) ; } if(ph1->IsDisp2OK()){ FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsDisp2OK()){ FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv2,w) ; } if(ph1->IsDispOK() && ph1->IsCPVOK()){ FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsDispOK() && ph2->IsCPVOK()){ FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv2,w) ; } if(ph1->IsDisp2OK() && ph1->IsCPVOK()){ FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsDisp2OK() && ph2->IsCPVOK()){ FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv2,w) ; } if(aIsCPVOK() && ph2->IsCPVOK()){ //pi0 CPV FillHistogram(Form("hPi0CPV_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hPi0CPVcore_cen%d",fCenBin),mv,ptv,w) ; if(aIsCPV2OK() && ph2->IsCPV2OK()){ FillHistogram(Form("hPi0CPV2_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; if(aIsDisp2OK() && ph2->IsDisp2OK()){ FillHistogram(Form("hPi0Disp2_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hPi0Disp2core_cen%d",fCenBin),mv,ptv,w) ; if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 Both FillHistogram(Form("hPi0Both2_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hPi0Both2core_cen%d",fCenBin),mv,ptv,w) ; if(aIsDispOK() && ph2->IsDispOK()){ //pi0 Disp FillHistogram(Form("hPi0Disp_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hPi0Dispcore_cen%d",fCenBin),mv,ptv,w) ; if(aIsCPVOK() && ph2->IsCPVOK()){ //pi0 Both FillHistogram(Form("hPi0Both_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hPi0Bothcore_cen%d",fCenBin),mv,ptv,w) ; if(aAt(i1) ; Double_t w1=ph1->GetWeight() ; for(Int_t ev=0; evGetSize();ev++){ TClonesArray * mixPHOS = static_cast(prevPHOS->At(ev)) ; for(Int_t i2=0; i2GetEntriesFast();i2++){ AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ; Double_t w2=ph2->GetWeight() ; TLorentzVector p12 = *ph1 + *ph2; TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2()); Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ; Double_t m = p12.M() ; Double_t pt = p12.Pt() ; Double_t mv = p12.M() ; Double_t ptv = p12.Pt() ; Double_t w=w1*w2 ; FillHistogram(Form("hMiPi0All_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hMiPi0Allcore_cen%d",fCenBin),mv,ptv,w) ; if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){ FillHistogram(Form("hMiPi0Allwou_cen%d",fCenBin),m,pt,w) ; } Double_t pt1 = ph1->Pt() ; Double_t ptv1 = ph1->GetMomV2()->Pt() ; Double_t pt2 = ph2->Pt() ; Double_t ptv2 = ph2->GetMomV2()->Pt() ; FillHistogram(Form("hMiSingleAll_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hMiSingleAll_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv1,w) ; FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv2,w) ; if(ph1->IsntUnfolded()) FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt1,w) ; if(ph2->IsntUnfolded()) FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt2,w) ; if(ph1->IsCPVOK()){ FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsCPVOK()){ FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv2,w) ; } if(ph1->IsCPV2OK()){ FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt1,w) ; } if(ph2->IsCPV2OK()){ FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt2,w) ; } if(ph1->IsDispOK()){ FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt1,w) ; if(ph1->IsntUnfolded()){ FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt1,w) ; } FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsDispOK()){ FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt2,w) ; if(ph2->IsntUnfolded()){ FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt2,w) ; } FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv2,w) ; } if(ph1->IsDisp2OK()){ FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsDisp2OK()){ FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv2,w) ; } if(ph1->IsDispOK() && ph1->IsCPVOK()){ FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsDispOK() && ph2->IsCPVOK()){ FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv2,w) ; } if(ph1->IsDisp2OK() && ph1->IsCPVOK()){ FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt1,w) ; FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv1,w) ; } if(ph2->IsDisp2OK() && ph2->IsCPVOK()){ FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt2,w) ; FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv2,w) ; } if(aIsCPVOK() && ph2->IsCPVOK()){ FillHistogram(Form("hMiPi0CPV_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hMiPi0CPVcore_cen%d",fCenBin),mv,ptv,w) ; if(aIsCPV2OK() && ph2->IsCPV2OK()){ FillHistogram(Form("hMiPi0CPV2_cen%d",fCenBin),m,pt,w) ; if(aIsDisp2OK() && ph2->IsDisp2OK()){ FillHistogram(Form("hMiPi0Disp2_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hMiPi0Disp2core_cen%d",fCenBin),mv,ptv,w) ; if(ph1->IsCPVOK() && ph2->IsCPVOK()){ FillHistogram(Form("hMiPi0Both2_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hMiPi0Both2core_cen%d",fCenBin),mv,ptv,w) ; } } if(ph1->IsDispOK() && ph2->IsDispOK()){ FillHistogram(Form("hMiPi0Disp_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hMiPi0Dispcore_cen%d",fCenBin),mv,ptv,w) ; if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){ FillHistogram(Form("hMiPi0Dispwou_cen%d",fCenBin),m,pt,w) ; } if(aIsCPVOK() && ph2->IsCPVOK()){ FillHistogram(Form("hMiPi0Both_cen%d",fCenBin),m,pt,w) ; FillHistogram(Form("hMiPi0Bothcore_cen%d",fCenBin),mv,ptv,w) ; if(aGetEntriesFast()>0){ prevPHOS->AddFirst(fPHOSEvent) ; fPHOSEvent=0; if(prevPHOS->GetSize()>kMixEvents[fCenBin]){//Remove redundant events TClonesArray * tmp = static_cast(prevPHOS->Last()) ; prevPHOS->RemoveLast() ; delete tmp ; } } // Post output data. PostData(1, fOutputContainer); fEventCounter++; } //________________________________________________________________________ void AliPHOSHijingEfficiency::FillSecondaries(){ //Sort secondaires const Double_t rcut = 1. ; //Fill spectra of primary particles //with proper weight for(Int_t i=0; iGetNtrack(); i++){ TParticle * p = fStack->Particle(i) ; if(p->R()>rcut) continue ; if(TMath::Abs(p->Y())>0.5) continue ; Double_t w = PrimaryParticleWeight(p) ; Int_t primPdgCode=p->GetPdgCode() ; switch(primPdgCode){ case 22: FillHistogram(Form("hPrimPhot_cen%d",fCenBin),p->Pt(),w); break ; case 11: case -11: FillHistogram(Form("hPrimEl_cen%d",fCenBin),p->Pt(),w); break ; case 111: FillHistogram(Form("hPrimPi0_cen%d",fCenBin),p->Pt(),w); break ; case 221: FillHistogram(Form("hPrimEta_cen%d",fCenBin),p->Pt(),w); break ; case 211: case -211: FillHistogram(Form("hPrimPipm_cen%d",fCenBin),p->Pt(),w); break ; case 2212: //p FillHistogram(Form("hPrimP_cen%d",fCenBin),p->Pt(),w); break ; case -2212: //pbar FillHistogram(Form("hPrimPbar_cen%d",fCenBin),p->Pt(),w); break ; case 2112: //n FillHistogram(Form("hPrimN_cen%d",fCenBin),p->Pt(),w); break ; case -2112: //nbar FillHistogram(Form("hPrimNbar_cen%d",fCenBin),p->Pt(),w); break ; case 310: //nbar FillHistogram(Form("hPrimK0S_cen%d",fCenBin),p->Pt(),w); break ; case 130: //nbar FillHistogram(Form("hPrimK0L_cen%d",fCenBin),p->Pt(),w); break ; case 321: //K+ case -321: //K- FillHistogram(Form("hPrimKpm_cen%d",fCenBin),p->Pt(),w); break ; default: //other FillHistogram(Form("hPrimOther_cen%d",fCenBin),p->Pt(),w); } } //Origins of secondary pi0s for(Int_t i=0; iGetNtrack(); i++){ TParticle * p = fStack->Particle(i) ; if(p->GetPdgCode()!=111) continue ; FillHistogram("Vertex",p->Pt(),p->R()); if(p->R()Phi() ; while(phi<0.)phi+=TMath::TwoPi() ; while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ; FillHistogram("hSecondPi0RphiZ",p->R(),phi,p->Vz()) ; Double_t w = PrimaryParticleWeight(p) ; FillHistogram("hSecondPi0RE",p->R(),p->Pt(),w) ; } TLorentzVector p1; Int_t inPHOS=fPHOSEvent->GetEntries() ; for (Int_t i1=0; i1At(i1) ; Double_t w1=ph1->GetWeight() ; for (Int_t i2=i1+1; i2At(i2) ; TLorentzVector p12 = *ph1 + *ph2; Double_t w2=ph2->GetWeight() ; Double_t w = TMath::Sqrt(w1*w2) ; FillHistogram(Form("hParentAll_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; Int_t prim=FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary()) ; if(prim>-1){ TParticle * particle = (TParticle *)fStack->Particle(prim); FillHistogram("hMass_R",p12.M(),p12.Pt(),TMath::Sqrt(particle->R()*particle->R()+particle->Vz()*particle->Vz())) ; Int_t pdgCode=particle->GetPdgCode() ; if(pdgCode!=111){ //common parent - not pi0 if(pdgCode==22) FillHistogram(Form("hParentGamma_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; else{ if(pdgCode==11 || pdgCode==-11){ FillHistogram(Form("hParentEl_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; } else{ if(InPi0mass(p12.M() ,p12.Pt())){ printf("Common parent: %d \n",pdgCode) ; } FillHistogram(Form("hParentOther_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; } }//Not photons }//Common parent not pi0 else{ //common parent - pi0 FillHistogram(Form("hParentPi0_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; FillHistogram(Form("Real_pi_R"),p12.M(),p12.Pt(),particle->R(),w) ; FillHistogram(Form("Real_pi_Z"),p12.M(),p12.Pt(),particle->Vz(),w) ; if(particle->R()Vz())<10.){ FillHistogram(Form("hParentDirPi0_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; continue ; } //Common particle pi0, created off-vertex Int_t primPi0=particle->GetFirstMother(); if(primPi0==-1){ FillHistogram(Form("hParentPi0NoPrim_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; } else{ Int_t primPdgCode=fStack->Particle(primPi0)->GetPdgCode() ; switch(primPdgCode){ case 221: FillHistogram(Form("hParentPi0Eta_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //eta break ; case 223: FillHistogram(Form("hParentPi0Omega_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //omega break ; case 211: //pi+- case -211: FillHistogram(Form("hParentPi0Pipm_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // break ; case 321: //K+- case -321: FillHistogram(Form("hParentPi0Kpm_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // break ; case 310: FillHistogram(Form("hParentPi0Ks_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // K0S break ; case 130: FillHistogram(Form("hParentPi0Kl_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // K0L break ; case 2212: //p case 2112: //n FillHistogram(Form("hParentPi0pn_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // pn break ; case -2212: //pbar case -2112: //nbar FillHistogram(Form("hParentPi0antipn_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // pn break ; default: //other FillHistogram(Form("hParentPi0Other_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // }//switch }//pi0 with primary }//common parent - pi0 }//there is common primary }//seond photon loop }//first photon loop //Now look at photon contaiminations for (Int_t i1=0; i1At(i1) ; Int_t iprim = ph1->GetPrimary() ; if(iprim<0) FillAllHistograms(Form("hGammaNoPrim_cen%d",fCenBin),ph1) ; // else{ //Find primary at vertex TParticle * primPHOS = fStack->Particle(iprim) ; Int_t iprimV=primPHOS->GetFirstMother(); TParticle * primVtx = primPHOS ; while((iprimV>-1) && primVtx->R()>rcut){ primVtx = fStack->Particle(iprimV) ; iprimV=primVtx->GetFirstMother(); } //photon Int_t primPdgCode=primVtx->GetPdgCode() ; switch(primPdgCode){ case 22: FillAllHistograms("hGammaPhot",ph1); break ; case 11: case -11: FillAllHistograms("hGammaEl",ph1); break ; case 111: FillAllHistograms("hGammaPi0",ph1); break ; case 221: FillAllHistograms("hGammaEta",ph1); break ; case 223: FillAllHistograms("hGammaOmega",ph1) ; //omega break ; case 211: case -211: FillAllHistograms("hGammaPipm",ph1); //Find particle entered PHOS if(primVtx == primPHOS) FillAllHistograms("hGammaPipmDirect",ph1); else{ Int_t primPdgPHOS=primPHOS->GetPdgCode() ; if(primPdgPHOS==22){ FillAllHistograms("hGammaPipmGamma",ph1); FillHistogram("hPipmGammaConvR",ph1->Pt(),primPHOS->R()); FillHistogram("hPipmGammaConvRZ",primPHOS->Vz(),primPHOS->R()); break ; } if(TMath::Abs(primPdgPHOS)==11){ FillAllHistograms("hGammaPipmEl",ph1); FillHistogram("hPipmElConvR",ph1->Pt(),primPHOS->R()); break ; } if(TMath::Abs(primPdgPHOS)==2212){ FillAllHistograms("hGammaPipmp",ph1); FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R()); break ; } if(TMath::Abs(primPdgPHOS)==2112){ FillAllHistograms("hGammaPipmn",ph1); FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R()); break ; } FillAllHistograms("hGammaPipmOther",ph1); FillHistogram("hPipmOtherConvR",ph1->Pt(),primPHOS->R()); } break ; case 2212: //p FillAllHistograms("hGammaP",ph1); break ; case -2212: //pbar FillAllHistograms("hGammaPbar",ph1); break ; case 2112: //n FillAllHistograms("hGammaN",ph1); break ; case -2112: //nbar FillAllHistograms("hGammaNbar",ph1) ; // pn break ; case 310: //nbar FillAllHistograms("hGammaK0S",ph1) ; // pn break ; case 130: //nbar FillAllHistograms("hGammaK0L",ph1) ; // pn break ; case 321: //K+ case -321: //K- FillAllHistograms("hGammaKpm",ph1) ; // pn break ; case -323: case 323: case -313: case 313: FillAllHistograms("hGammaKstar",ph1) ; // K*(892) break ; case -2224 : //Deltas case 2224 : //Deltas case -2214 : //Deltas case 2214 : //Deltas case -2114 : //Deltas case 2114 : //Deltas case -1114 : //Deltas case 1114 : //Deltas FillAllHistograms("hGammaDelta",ph1) ; // pn break ; default: //other if(primVtx->GetPDG()->Charge()) FillAllHistograms("hGammaOtherCharged",ph1) ; // else FillAllHistograms("hGammaOtherNeutral",ph1) ; // } } }//single photons } //________________________________________________________________________ void AliPHOSHijingEfficiency::Terminate(Option_t *) { // Draw result to the screen // Called once at the end of the query } //________________________________________________________________________ Bool_t AliPHOSHijingEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz) { //Check if this channel belogs to the good ones if(strcmp(det,"PHOS")==0){ if(mod>5 || mod<1){ AliError(Form("No bad map for PHOS module %d ",mod)) ; return kTRUE ; } if(!fPHOSBadMap[mod]){ AliError(Form("No Bad map for PHOS module %d",mod)) ; return kTRUE ; } if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0) return kFALSE ; else return kTRUE ; } else{ AliError(Form("Can not find bad channels for detector %s ",det)) ; } return kTRUE ; } //_____________________________________________________________________________ void AliPHOSHijingEfficiency::FillAllHistograms(const char * key,AliCaloPhoton * ph)const{ //Fill All PID histograms Double_t w=ph->GetWeight() ; Double_t pt = ph->Pt() ; Double_t ptC=ph->GetMomV2()->Pt() ; FillHistogram(Form("%s_All_cen%d",key,fCenBin),pt,w) ; FillHistogram(Form("%s_Allcore_cen%d",key,fCenBin),ptC,w) ; if(ph->IsCPVOK()){ FillHistogram(Form("%s_CPV_cen%d",key,fCenBin),pt,w) ; FillHistogram(Form("%s_CPVcore_cen%d",key,fCenBin),ptC,w) ; } if(ph->IsCPV2OK()){ FillHistogram(Form("%s_CPV2_cen%d",key,fCenBin),pt,w) ; FillHistogram(Form("%s_CPV2core_cen%d",key,fCenBin),ptC,w) ; } if(ph->IsDispOK()){ FillHistogram(Form("%s_Disp_cen%d",key,fCenBin),pt,w) ; FillHistogram(Form("%s_Dispcore_cen%d",key,fCenBin),ptC,w) ; if(ph->IsDisp2OK()){ FillHistogram(Form("%s_Disp2_cen%d",key,fCenBin),pt,w) ; FillHistogram(Form("%s_Disp2core_cen%d",key,fCenBin),ptC,w) ; } if(ph->IsCPVOK()){ FillHistogram(Form("%s_Both_cen%d",key,fCenBin),pt,w) ; FillHistogram(Form("%s_Bothcore_cen%d",key,fCenBin),ptC,w) ; } } } //_____________________________________________________________________________ void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x)const{ //FillHistogram TObject * tmp = fOutputContainer->FindObject(key) ; if(!tmp){ AliInfo(Form("can not find histogram <%s> ",key)) ; return ; } if(tmp->IsA() == TClass::GetClass("TH1I")){ ((TH1I*)tmp)->Fill(x) ; return ; } if(tmp->IsA() == TClass::GetClass("TH1F")){ ((TH1F*)tmp)->Fill(x) ; return ; } if(tmp->IsA() == TClass::GetClass("TH1D")){ ((TH1D*)tmp)->Fill(x) ; return ; } AliInfo(Form("can not find 1D histogram <%s> ",key)) ; } //_____________________________________________________________________________ void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{ //FillHistogram TObject * tmp = fOutputContainer->FindObject(key) ; if(!tmp){ AliInfo(Form("can not find histogram <%s> ",key)) ; return ; } if(tmp->IsA() == TClass::GetClass("TH1F")){ ((TH1F*)tmp)->Fill(x,y) ; return ; } if(tmp->IsA() == TClass::GetClass("TH2F")){ ((TH2F*)tmp)->Fill(x,y) ; return ; } AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ; } //_____________________________________________________________________________ void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{ //Fills 1D histograms with key TObject * tmp = fOutputContainer->FindObject(key) ; if(!tmp){ AliInfo(Form("can not find histogram <%s> ",key)) ; return ; } if(tmp->IsA() == TClass::GetClass("TH2F")){ ((TH2F*)tmp)->Fill(x,y,z) ; return ; } if(tmp->IsA() == TClass::GetClass("TH3F")){ ((TH3F*)tmp)->Fill(x,y,z) ; return ; } } //_____________________________________________________________________________ void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z, Double_t w) const{ //Fills 1D histograms with key TObject * tmp = fOutputContainer->FindObject(key) ; if(!tmp){ AliInfo(Form("can not find histogram <%s> ",key)) ; return ; } if(tmp->IsA() == TClass::GetClass("TH3F")){ ((TH3F*)tmp)->Fill(x,y,z,w) ; return ; } } //___________________________________________________________________________ Int_t AliPHOSHijingEfficiency::ConvertRunNumber(Int_t run){ switch(run){ case 139517 : return 137; case 139514 : return 136; case 139513 : return 135; case 139511 : return 134; case 139510 : return 133; case 139507 : return 132; case 139505 : return 131; case 139504 : return 130; case 139503 : return 129; case 139470 : return 128; case 139467 : return 127; case 139466 : return 126; case 139465 : return 125; case 139440 : return 124; case 139439 : return 123; case 139438 : return 122; case 139437 : return 121; case 139360 : return 120; case 139329 : return 119; case 139328 : return 118; case 139314 : return 117; case 139311 : return 116; case 139310 : return 115; case 139309 : return 114; case 139308 : return 113; case 139173 : return 112; case 139172 : return 111; case 139110 : return 110; case 139107 : return 109; case 139105 : return 108; case 139104 : return 107; case 139042 : return 106; case 139038 : return 105; case 139037 : return 104; case 139036 : return 103; case 139029 : return 102; case 139028 : return 101; case 138983 : return 100; case 138982 : return 99; case 138980 : return 98; case 138979 : return 97; case 138978 : return 96; case 138977 : return 95; case 138976 : return 94; case 138973 : return 93; case 138972 : return 92; case 138965 : return 91; case 138924 : return 90; case 138872 : return 89; case 138871 : return 88; case 138870 : return 87; case 138837 : return 86; case 138830 : return 85; case 138828 : return 84; case 138826 : return 83; case 138796 : return 82; case 138795 : return 81; case 138742 : return 80; case 138732 : return 79; case 138730 : return 78; case 138666 : return 77; case 138662 : return 76; case 138653 : return 75; case 138652 : return 74; case 138638 : return 73; case 138624 : return 72; case 138621 : return 71; case 138583 : return 70; case 138582 : return 69; case 138579 : return 68; case 138578 : return 67; case 138534 : return 66; case 138469 : return 65; case 138442 : return 64; case 138439 : return 63; case 138438 : return 62; case 138396 : return 61; case 138364 : return 60; case 138359 : return 59; case 138275 : return 58; case 138225 : return 57; case 138201 : return 56; case 138200 : return 55; case 138197 : return 54; case 138192 : return 53; case 138190 : return 52; case 138154 : return 51; case 138153 : return 50; case 138151 : return 49; case 138150 : return 48; case 138126 : return 47; case 138125 : return 46; case 137848 : return 45; case 137847 : return 44; case 137844 : return 43; case 137843 : return 42; case 137752 : return 41; case 137751 : return 40; case 137748 : return 39; case 137724 : return 38; case 137722 : return 37; case 137718 : return 36; case 137704 : return 35; case 137693 : return 34; case 137692 : return 33; case 137691 : return 32; case 137689 : return 31; case 137686 : return 30; case 137685 : return 29; case 137639 : return 28; case 137638 : return 27; case 137608 : return 26; case 137595 : return 25; case 137549 : return 24; case 137546 : return 23; case 137544 : return 22; case 137541 : return 21; case 137539 : return 20; case 137531 : return 19; case 137530 : return 18; case 137443 : return 17; case 137441 : return 16; case 137440 : return 15; case 137439 : return 14; case 137434 : return 13; case 137432 : return 12; case 137431 : return 11; case 137430 : return 10; case 137366 : return 9; case 137243 : return 8; case 137236 : return 7; case 137235 : return 6; case 137232 : return 5; case 137231 : return 4; case 137165 : return 3; case 137162 : return 2; case 137161 : return 1; default : return 199; } } //_____________________________________________________________________________ Bool_t AliPHOSHijingEfficiency::TestLambda(Double_t pt,Double_t l1,Double_t l2){ //For R=3.5 Double_t l1Mean = 1.170014 -0.059465/(1.+0.019343*pt+0.147544*pt*pt) ; Double_t l2Mean = 1.626270 + 0.761554*exp(-1.213839*pt)-0.020027*pt ; Double_t l1Sigma = 0.133409 + 0.261307*exp(-0.636874*pt)-0.002849*pt ; Double_t l2Sigma = 0.289698 + 0.459400*exp(-1.214242*pt)-0.012578*pt ; Double_t c=-0.124103 ; /* Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ; Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ; Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt; Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt; Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ; */ Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma + 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma + 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ; return (R2<2.5*2.5) ; } //_____________________________________________________________________________ Bool_t AliPHOSHijingEfficiency::TestLambda2(Double_t pt,Double_t l1,Double_t l2){ //For R=4.5 Double_t l1Mean = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ; Double_t l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ; Double_t l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ; Double_t l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ; Double_t c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ; /* Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ; Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ; Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt; Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt; Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ; */ Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma + 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma + 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ; return (R2<2.5*2.5) ; } //____________________________________________________________________________ Double_t AliPHOSHijingEfficiency::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){ //Parameterization of LHC10h period //_true if neutral_ Double_t meanX=0; Double_t meanZ=0.; Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+ 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219); Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ; AliESDEvent *event = dynamic_cast(InputEvent()); if(!event) //not possible, but required by Coverity return 999.; Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for -- if(mf<0.){ //field -- meanZ = -0.468318 ; if(charge>0) meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ; else meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ; } else{ //Field ++ meanZ= -0.468318; if(charge>0) meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ; else meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ; } Double_t rz=(dz-meanZ)/sz ; Double_t rx=(dx-meanX)/sx ; return TMath::Sqrt(rx*rx+rz*rz) ; } //____________________________________________________________________________ Bool_t AliPHOSHijingEfficiency::TestTOF(Double_t t, Double_t e){ Double_t sigma = TMath::Sqrt(2.23183e-09*2.23183e-09 +2.24611e-09*2.24611e-09/e +5.65054e-09*5.65054e-09/e/e) ; sigma=TMath::Min(20.e-9,sigma) ; //for the soft (<400 MeV) part Double_t mean=1.1e-9 ; if(TMath::Abs(t-mean)<2.*sigma) return kTRUE ; else if(TMath::Abs(t-mean+100.e-9)<2.*sigma) return kTRUE ; return kFALSE ; } //____________________________________________________________________________ Double_t AliPHOSHijingEfficiency::CoreEnergy(AliESDCaloCluster * clu){ //calculate energy of the cluster in the circle with radius distanceCut around the maximum //Can not use already calculated coordinates? //They have incidence correction... const Double_t distanceCut =3.5 ; const Double_t logWeight=4.5 ; Double32_t * elist = clu->GetCellsAmplitudeFraction() ; // Calculates the center of gravity in the local PHOS-module coordinates Float_t wtot = 0; Double_t xc[100]={0} ; Double_t zc[100]={0} ; Double_t x = 0 ; Double_t z = 0 ; Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ; for(Int_t iDigit=0; iDigitAbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ; fPHOSGeo->RelPosInModule(relid, xi, zi); xc[iDigit]=xi ; zc[iDigit]=zi ; if (clu->E()>0 && elist[iDigit]>0) { Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ; x += xc[iDigit] * w ; z += zc[iDigit] * w ; wtot += w ; } } if (wtot>0) { x /= wtot ; z /= wtot ; } Double_t coreE=0. ; for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) { Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ; if(distance < distanceCut) coreE += elist[iDigit] ; } //Apply non-linearity correction return rnlin(coreE) ; } //____________________________________________________________________________ Bool_t AliPHOSHijingEfficiency::AreNeibors(Int_t id1,Int_t id2){ Int_t relid1[4] ; fPHOSGeo->AbsToRelNumbering(id1, relid1) ; Int_t relid2[4] ; fPHOSGeo->AbsToRelNumbering(id2, relid2) ; if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ; Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ; if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side return 1 ; } else { if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) return 0; // Difference in row numbers is too large to look further } return 0 ; } return 0 ; } //____________________________________________________________________________ void AliPHOSHijingEfficiency::Reclusterize(AliESDCaloCluster * clu){ //Re-clusterize to make continues cluster const Int_t oldMulDigit=clu->GetNCells() ; Double32_t * elist = clu->GetCellsAmplitudeFraction() ; UShort_t * dlist = clu->GetCellsAbsId(); Int_t index[oldMulDigit] ; Bool_t used[oldMulDigit] ; for(Int_t i=0; iSetNCells(inClu); //copy UShort_t tmpD[oldMulDigit] ; Double_t tmpE[oldMulDigit] ; for(Int_t i=0; iGetCellsAmplitudeFraction() ; // Calculates the center of gravity in the local PHOS-module coordinates Float_t wtot = 0; Double_t xc[100]={0} ; Double_t zc[100]={0} ; Double_t x = 0 ; Double_t z = 0 ; Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ; const Double_t logWeight=4.5 ; for(Int_t iDigit=0; iDigitAbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ; fPHOSGeo->RelPosInModule(relid, xi, zi); xc[iDigit]=xi ; zc[iDigit]=zi ; if (clu->E()>0 && elist[iDigit]>0) { Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ; x += xc[iDigit] * w ; z += zc[iDigit] * w ; wtot += w ; } } if (wtot>0) { x /= wtot ; z /= wtot ; } wtot = 0. ; Double_t dxx = 0.; Double_t dzz = 0.; Double_t dxz = 0.; Double_t xCut = 0. ; Double_t zCut = 0. ; for(Int_t iDigit=0; iDigitE()>0 && elist[iDigit]>0.) { Double_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ; Double_t xi= xc[iDigit] ; Double_t zi= zc[iDigit] ; if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){ xCut += w * xi ; zCut += w * zi ; dxx += w * xi * xi ; dzz += w * zi * zi ; dxz += w * xi * zi ; wtot += w ; } } } if (wtot>0) { xCut/= wtot ; zCut/= wtot ; dxx /= wtot ; dzz /= wtot ; dxz /= wtot ; dxx -= xCut * xCut ; dzz -= zCut * zCut ; dxz -= xCut * zCut ; m02 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; m20 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; } else { m20=m02=0.; } } //___________________________________________________________________________ void AliPHOSHijingEfficiency::ProcessMC(){ //fill histograms for efficiensy etc. calculation const Double_t rcut = 1. ; //cut for primary particles //---------First pi0/eta----------------------------- char partName[10] ; char hkey[55] ; if(!fStack) return ; for(Int_t i=0;iGetNtrack();i++){ TParticle* particle = fStack->Particle(i); if(particle->GetPdgCode() == 111) snprintf(partName,10,"pi0") ; else if(particle->GetPdgCode() == 221) snprintf(partName,10,"eta") ; else if(particle->GetPdgCode() == 22) snprintf(partName,10,"gamma") ; else continue ; //Primary particle Double_t r=particle->R() ; Double_t pt = particle->Pt() ; //Distribution over vertex FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ; if(r >rcut) continue ; //Total number of pi0 with creation radius <1 cm Double_t w = PrimaryParticleWeight(particle) ; snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCenBin) ; FillHistogram(hkey,pt,w) ; if(TMath::Abs(particle->Y())<0.12){ snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCenBin) ; FillHistogram(hkey,pt,w) ; } snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCenBin) ; FillHistogram(hkey,particle->Y(),w) ; Double_t phi=particle->Phi() ; while(phi<0.)phi+=TMath::TwoPi() ; while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ; snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCenBin) ; FillHistogram(hkey,phi,w) ; } } //___________________________________________________________________________ Int_t AliPHOSHijingEfficiency::FindPrimary(AliESDCaloCluster*clu, Bool_t&sure){ //Finds primary and estimates if it unique one? //First check can it be photon/electron const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle Int_t n=clu->GetNLabels() ; for(Int_t i=0; iParticle(clu->GetLabelAt(i)) ; Int_t pdg = p->GetPdgCode() ; if(pdg==22 || pdg==11 || pdg == -11){ if(p->Energy()>emFraction*clu->E()){ sure=kTRUE ; return clu->GetLabelAt(i); } } } Double_t* Ekin= new Double_t[n] ; for(Int_t i=0; iParticle(clu->GetLabelAt(i)) ; Ekin[i]=p->P() ; // estimate of kinetic energy if(p->GetPdgCode()==-2212 || p->GetPdgCode()==-2112){ Ekin[i]+=1.8 ; //due to annihilation } } Int_t iMax=0; Double_t eMax=0.,eSubMax=0. ; for(Int_t i=0; ieMax){ eSubMax=eMax; eMax=Ekin[i]; iMax=i; } } if(eSubMax>0.8*eMax)//not obvious primary sure=kFALSE; else sure=kTRUE; delete[] Ekin; return clu->GetLabelAt(iMax) ; } //___________________________________________________________________________ Double_t AliPHOSHijingEfficiency::PrimaryWeight(Int_t primary){ //Check who is the primary and introduce weight to correct primary spectrum if(primary<0 || primary>=fStack->GetNtrack()) return 1 ; //trace primaries up to IP TParticle* particle = fStack->Particle(primary); Double_t r=particle->R() ; Int_t mother = particle->GetFirstMother() ; while(mother>-1){ if(r<1. && particle->GetPdgCode()==111) break ; particle = fStack->Particle(mother); mother = particle->GetFirstMother() ; r=particle->R() ; } return TMath::Max(0.,PrimaryParticleWeight(particle)) ; } //________________________________________________________________________ Double_t AliPHOSHijingEfficiency::PrimaryParticleWeight(TParticle * particle){ Int_t pdg = particle->GetPdgCode() ; Int_t type=0 ; if(pdg == 111 || TMath::Abs(pdg)==211){ type =1 ; } else{ if(TMath::Abs(pdg)<1000){ //Kaon-like type =2 ; } else type = 3; //baryons } Double_t x = particle->Pt() ; if(type==1){ if(fCenBin==0) //0-5 return (1.662990+1.140890*x-0.192088*x*x)/(1.-0.806630*x+0.304771*x*x)+0.141690*x ; if(fCenBin==1) //5-10 return (1.474351+0.791492*x-0.066369*x*x)/(1.-0.839338*x+0.317312*x*x)+0.093289*x ; if(fCenBin==2) //10-20 return (1.174728+0.959681*x-0.137695*x*x)/(1.-0.788873*x+0.299538*x*x)+0.128759*x ; if(fCenBin==3) //20-40 return (0.927335+0.475349*x+0.004364*x*x)/(1.-0.817966*x+0.309787*x*x)+0.086899*x ; if(fCenBin==4) //40-60 return (0.676878+0.190680*x+0.077031*x*x)/(1.-0.790623*x+0.305183*x*x)+0.064510*x ; if(fCenBin==5) //60-80 return (0.684726-0.606262*x+0.409963*x*x)/(1.-1.080061*x+0.456933*x*x)+0.005151*x ; } if(type==2){ if(fCenBin==0) //0-5 return (-0.417131+2.253936*x-0.337731*x*x)/(1.-0.909892*x+0.316820*x*x)+0.157312*x ; if(fCenBin==1) //5-10 return (-0.352275+1.844466*x-0.248598*x*x)/(1.-0.897048*x+0.316462*x*x)+0.132461*x ; if(fCenBin==2) //10-20 return (-0.475481+1.975108*x-0.336013*x*x)/(1.-0.801028*x+0.276705*x*x)+0.188164*x ; if(fCenBin==3) //20-40 return (-0.198954+1.068789*x-0.103540*x*x)/(1.-0.848354*x+0.299209*x*x)+0.112939*x ; if(fCenBin==4) //40-60 return (-0.111052+0.664041*x-0.019717*x*x)/(1.-0.804916*x+0.300779*x*x)+0.095784*x ; if(fCenBin==5) //0-5 return (0.202788-0.439832*x+0.564585*x*x)/(1.-1.254029*x+0.679444*x*x)+0.016235*x ; } if(type==3){ if(fCenBin==0) //0-5 return (-1.312732+2.743568*x-0.375775*x*x)/(1.-0.717533*x+0.164694*x*x)+0.164445*x ; if(fCenBin==1) //5-10 return (-1.229425+2.585889*x-0.330164*x*x)/(1.-0.715892*x+0.167386*x*x)+0.133085*x ; if(fCenBin==2) //10-20 return (-1.135677+2.397489*x-0.320355*x*x)/(1.-0.709312*x+0.164350*x*x)+0.146095*x ; if(fCenBin==3) //20-40 return (-0.889993+1.928263*x-0.220785*x*x)/(1.-0.715991*x+0.174729*x*x)+0.095098*x ; if(fCenBin==4) //40-60 return (-0.539237+1.329118*x-0.115439*x*x)/(1.-0.722906*x+0.186832*x*x)+0.059267*x ; if(fCenBin==5) //60-80 return (-0.518126+1.327628*x-0.130881*x*x)/(1.-0.665649*x+0.184300*x*x)+0.081701*x ; } return 1. ; } //________________________________________________________________________ Int_t AliPHOSHijingEfficiency::FindCommonParent(Int_t iPart, Int_t jPart){ //check if there is a common parent for particles i and j // -1: no common parent or wrong iPart/jPart if(iPart==-1 || iPart>=fStack->GetNtrack() || jPart==-1 || jPart>=fStack->GetNtrack()) return -1; Int_t iprim1=iPart; while(iprim1>-1){ Int_t iprim2=jPart; while(iprim2>-1){ if(iprim1==iprim2) return iprim1 ; iprim2=((TParticle *)fStack->Particle(iprim2))->GetFirstMother(); } iprim1=((TParticle *)fStack->Particle(iprim1))->GetFirstMother(); } return -1; } //________________________________________________________________________ Bool_t AliPHOSHijingEfficiency::HaveParent(Int_t iPart, Int_t pdgParent){ //check if there is a common parent for particles i and j // -1: no common parent or wrong iPart/jPart if(iPart==-1 || iPart>=fStack->GetNtrack()) return -1; Int_t iprim1=iPart; while(iprim1>-1){ TParticle * tmp = fStack->Particle(iprim1) ; if(tmp->GetPdgCode()==pdgParent) return kTRUE ; iprim1=tmp->GetFirstMother(); } return kFALSE; } //________________________________________________________________________ Bool_t AliPHOSHijingEfficiency::InPi0mass(Double_t m, Double_t /*pt*/){ return TMath::Abs(m-0.135)<0.007*2.5 ; }