ClassImp(AliAnaPi0)
-//________________________________________________________________________________________________________________________________________________
+//______________________________________________________
AliAnaPi0::AliAnaPi0() : AliAnaCaloTrackCorrBaseClass(),
fEventsList(0x0),
fCalorimeter(""), fNModules(12),
fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
fFillBadDistHisto(kFALSE), fFillSSCombinations(kFALSE),
-fFillAngleHisto(kFALSE), fFillAsymmetryHisto(kFALSE), fFillOriginHisto(0),
+fFillAngleHisto(kFALSE), fFillAsymmetryHisto(kFALSE), fFillOriginHisto(0), fFillArmenterosThetaStar(0),
//Histograms
fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0),
fhEDensityCluster(0), fhEDensityCell(0), fhEDensityCellvsCluster(0),
fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(),
fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
-fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0)
+fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0),
+fhCosThStarPrimPi0(0), fhCosThStarPrimEta(0)//,
{
-//Default Ctor
- InitParameters();
+ //Default Ctor
+ InitParameters();
+
+ for(Int_t i = 0; i < 4; i++)
+ {
+ fhArmPrimEta[i] = 0;
+ fhArmPrimPi0[i] = 0;
+ }
}
-//________________________________________________________________________________________________________________________________________________
-AliAnaPi0::~AliAnaPi0() {
+//_____________________
+AliAnaPi0::~AliAnaPi0()
+{
// Remove event containers
- if(DoOwnMix() && fEventsList){
+ if(DoOwnMix() && fEventsList)
+ {
for(Int_t ic=0; ic<GetNCentrBin(); ic++)
{
for(Int_t iz=0; iz<GetNZvertBin(); iz++)
}
-//________________________________________________________________________________________________________________________________________________
+//______________________________
void AliAnaPi0::InitParameters()
{
//Init parameters when first called the analysis
}
-//________________________________________________________________________________________________________________________________________________
+//_______________________________________
TObjString * AliAnaPi0::GetAnalysisCuts()
{
//Save parameters used for analysis
return new TObjString(parList) ;
}
-//________________________________________________________________________________________________________________________________________________
+//_________________________________________
TList * AliAnaPi0::GetCreateOutputObjects()
{
// Create histograms to be saved in output file and
//create event containers
fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
-
+
for(Int_t ic=0; ic<GetNCentrBin(); ic++)
{
for(Int_t iz=0; iz<GetNZvertBin(); iz++)
outputContainer->Add(fhReSS[2]) ;
}
- fhEventBin=new TH1I("hEventBin","Number of real pairs per bin(cen,vz,rp)",
- GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
- GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
- fhEventBin->SetXTitle("bin");
- outputContainer->Add(fhEventBin) ;
-
- fhEventMixBin=new TH1I("hEventMixBin","Number of mixed pairs per bin(cen,vz,rp)",
- GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
- GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
- fhEventMixBin->SetXTitle("bin");
- outputContainer->Add(fhEventMixBin) ;
-
+ if(DoOwnMix())
+ {
+ fhEventBin=new TH1I("hEventBin","Number of real pairs per bin(cen,vz,rp)",
+ GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
+ GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
+ fhEventBin->SetXTitle("bin");
+ outputContainer->Add(fhEventBin) ;
+
+
+ fhEventMixBin=new TH1I("hEventMixBin","Number of mixed pairs per bin(cen,vz,rp)",
+ GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
+ GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
+ fhEventMixBin->SetXTitle("bin");
+ outputContainer->Add(fhEventMixBin) ;
+ }
+
if(GetNCentrBin()>1)
{
fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
}//loop combinations
} // SM combinations
+ if(fFillArmenterosThetaStar && IsDataMC())
+ {
+ TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
+ Int_t narmbins = 400;
+ Float_t armmin = 0;
+ Float_t armmax = 0.4;
+
+ for(Int_t i = 0; i < 4; i++)
+ {
+ fhArmPrimPi0[i] = new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
+ Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
+ 200, -1, 1, narmbins,armmin,armmax);
+ fhArmPrimPi0[i]->SetYTitle("p_{T}^{Arm}");
+ fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
+ outputContainer->Add(fhArmPrimPi0[i]) ;
+
+ fhArmPrimEta[i] = new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
+ Form("Armenteros of primary #eta, %s",ebin[i].Data()),
+ 200, -1, 1, narmbins,armmin,armmax);
+ fhArmPrimEta[i]->SetYTitle("p_{T}^{Arm}");
+ fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
+ outputContainer->Add(fhArmPrimEta[i]) ;
+
+ }
+
+ // Same as asymmetry ...
+ fhCosThStarPrimPi0 = new TH2F
+ ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
+ fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
+ fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
+ outputContainer->Add(fhCosThStarPrimPi0) ;
+
+ fhCosThStarPrimEta = new TH2F
+ ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
+ fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
+ fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
+ outputContainer->Add(fhCosThStarPrimEta) ;
+
+ }
+
// for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
//
// printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
//printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
// prim->GetName(), prim->GetPdgCode());
- if( pdg == 111 || pdg == 221){
+ if( pdg == 111 || pdg == 221)
+ {
Double_t pi0Pt = prim->Pt() ;
Double_t pi0E = prim->Energy() ;
if(pi0E == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
//printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
// phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
- TLorentzVector lv1, lv2;
+ TLorentzVector lv1, lv2,lvmeson;
phot1->Momentum(lv1);
phot2->Momentum(lv2);
+ prim ->Momentum(lvmeson);
+
+ if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
Bool_t inacceptance = kFALSE;
if(fCalorimeter == "PHOS")
AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
{
- TLorentzVector lv1, lv2;
+ TLorentzVector lv1, lv2,lvmeson;
lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
+ lvmeson.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
+ if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
+
Bool_t inacceptance = kFALSE;
if(fCalorimeter == "PHOS")
{
} // read AOD MC
}
-//_____________________________________________________________
-void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1, const Int_t index2,
- const Float_t pt1, const Float_t pt2,
- const Int_t ncell1, const Int_t ncell2,
- const Double_t mass, const Double_t pt, const Double_t asym,
- const Double_t deta, const Double_t dphi){
+//__________________________________________________________________________________
+void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg, TLorentzVector meson,
+ TLorentzVector daugh1, TLorentzVector daugh2)
+{
+ // Fill armenteros plots
+
+ // Get pTArm and AlphaArm
+ Float_t momentumSquaredMother = meson.P()*meson.P();
+ Float_t momentumDaughter1AlongMother = 0.;
+ Float_t momentumDaughter2AlongMother = 0.;
+
+ if (momentumSquaredMother > 0.)
+ {
+ momentumDaughter1AlongMother = (daugh1.Px()*meson.Px() + daugh1.Py()*meson.Py()+ daugh1.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
+ momentumDaughter2AlongMother = (daugh2.Px()*meson.Px() + daugh2.Py()*meson.Py()+ daugh2.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
+ }
+
+ Float_t momentumSquaredDaughter1 = daugh1.P()*daugh1.P();
+ Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
+
+ Float_t pTArm = 0.;
+ if (ptArmSquared > 0.)
+ pTArm = sqrt(ptArmSquared);
+
+ Float_t alphaArm = 0.;
+ if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
+ alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
+
+ TLorentzVector daugh1Boost = daugh1;
+ daugh1Boost.Boost(-meson.BoostVector());
+ Float_t cosThStar=TMath::Cos(daugh1Boost.Vect().Angle(meson.Vect()));
+
+ Float_t en = meson.Energy();
+ Int_t ebin = -1;
+ if(en > 8 && en <= 12) ebin = 0;
+ if(en > 12 && en <= 16) ebin = 1;
+ if(en > 16 && en <= 20) ebin = 2;
+ if(en > 20) ebin = 3;
+ if(ebin < 0 || ebin > 3) return ;
+
+
+ if(pdg==111)
+ {
+ fhCosThStarPrimPi0->Fill(en,cosThStar);
+ fhArmPrimPi0[ebin]->Fill(alphaArm,pTArm);
+ }
+ else
+ {
+ fhCosThStarPrimEta->Fill(en,cosThStar);
+ fhArmPrimEta[ebin]->Fill(alphaArm,pTArm);
+ }
+
+ if(GetDebug() > 2 )
+ {
+ Float_t asym = 0;
+ if(daugh1.E()+daugh2.E() > 0) asym = TMath::Abs(daugh1.E()-daugh2.E())/(daugh1.E()+daugh2.E());
+
+ printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
+ en,alphaArm,pTArm,cosThStar,asym);
+ }
+}
+
+//_______________________________________________________________________________________
+void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1, Int_t index2,
+ Float_t pt1, Float_t pt2,
+ Int_t ncell1, Int_t ncell2,
+ Double_t mass, Double_t pt, Double_t asym,
+ Double_t deta, Double_t dphi)
+{
//Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
//Adjusted for Pythia, need to see what to do for other generators.
//Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
}
}
-//____________________________________________________________________________________________________________________________________________________
+//__________________________________________
void AliAnaPi0::MakeAnalysisFillHistograms()
{
//Process one event and extract photons from AOD branch
//Int_t curVzBin = GetEventVzBin();
//Int_t curRPBin = GetEventRPBin();
Int_t eventbin = GetEventMixBin();
-
+
+ if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
+ {
+ printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mix Bin not exepcted: cen bin %d, z bin %d, rp bin %d, total bin %d, Event Centrality %d, z vertex %2.3f, Reaction Plane %2.3f\n",GetEventCentralityBin(),GetEventVzBin(), GetEventRPBin(),eventbin,GetEventCentrality(),vert[2],GetEventPlaneAngle());
+ return;
+ }
+
//Get shower shape information of clusters
TObjArray *clusters = 0;
if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
if (evtIndex1 != currentEvtIndex)
{
//Fill event bin info
- fhEventBin->Fill(eventbin) ;
+ if(DoOwnMix()) fhEventBin->Fill(eventbin) ;
if(GetNCentrBin() > 1)
{
fhCentrality->Fill(curCentrBin);
{
if(a < fAsymCuts[iasym])
{
- Int_t index = ((GetEventCentralityBin()*fNPIDBits)+ipid)*fNAsymCuts + iasym;
+ Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
- if(index < 0 || index >= curCentrBin*fNPIDBits*fNAsymCuts) continue ;
+ if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
fhMi1 [index]->Fill(pt,m) ;
if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
}
-//____________________________________________________________________________________________________________________________________________________
+//________________________________________________________________________
Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
{
// retieves the event index and checks the vertex
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff