fhEChargedNoOut(0),fhPtChargedNoOut(0),fhPhiChargedNoOut(0),fhEtaChargedNoOut(0), fhEtaPhiChargedNoOut(0),
fhDeltaE(0), fhDeltaPt(0),fhDeltaPhi(0),fhDeltaEta(0), fhRatioE(0), fhRatioPt(0),fhRatioPhi(0),fhRatioEta(0),
fh2E(0),fh2Pt(0),fh2Phi(0),fh2Eta(0), fhIM(0), fhIMCellCut(0),fhAsym(0),
- fhNCellsPerCluster(0),fhNCellsPerClusterMIP(0), fhNClusters(0), fhNCells(0),
- fhAmplitude(0), fhAmpId(0), fhTime(0), fhTimeId(0), fhTimeAmp(0),
+ fhNCellsPerCluster(0),fhNCellsPerClusterMIP(0), fhNClusters(0),
+ fhCellTimeSpreadRespectToCellMax(0),fhCellIdCellLargeTimeSpread(0),
+ fhNCells(0), fhAmplitude(0), fhAmpId(0), fhTime(0), fhTimeId(0), fhTimeAmp(0),
fhCaloCorrNClusters(0), fhCaloCorrEClusters(0), fhCaloCorrNCells(0), fhCaloCorrECells(0),
- fhEMod(0),fhNClustersMod(0), fhNCellsPerClusterMod(0), fhNCellsMod(0),
- fhGridCellsMod(0), fhGridCellsEMod(0), fhAmplitudeMod(0), fhIMMod(0), fhIMCellCutMod(0),
+ fhEMod(0), fhNClustersMod(0), fhNCellsPerClusterMod(0), fhNCellsMod(0),
+ fhGridCellsMod(0), fhGridCellsEMod(0), fhGridCellsTimeMod(0),
+ fhAmplitudeMod(0), fhAmplitudeModFraction(0),fhIMMod(0), fhIMCellCutMod(0),
fhGenGamPt(0),fhGenGamEta(0),fhGenGamPhi(0),fhGenPi0Pt(0),fhGenPi0Eta(0),fhGenPi0Phi(0),
fhGenEtaPt(0),fhGenEtaEta(0),fhGenEtaPhi(0),fhGenOmegaPt(0),fhGenOmegaEta(0),fhGenOmegaPhi(0),
fhGenElePt(0),fhGenEleEta(0),fhGenElePhi(0), fhEMVxyz(0), fhEMR(0), fhHaVxyz(0), fhHaR(0),
fhRatioE(qa.fhRatioE), fhRatioPt(qa.fhRatioPt), fhRatioPhi(qa.fhRatioPhi), fhRatioEta(qa.fhRatioEta),
fh2E(qa.fh2E), fh2Pt(qa.fh2Pt), fh2Phi(qa.fh2Phi),fh2Eta(qa.fh2Eta),
fhIM(qa.fhIM), fhIMCellCut(qa.fhIMCellCut), fhAsym(qa.fhAsym),
- fhNCellsPerCluster(qa.fhNCellsPerCluster), fhNCellsPerClusterMIP(qa.fhNCellsPerClusterMIP),
- fhNClusters(qa.fhNClusters), fhNCells(qa.fhNCells), fhAmplitude(qa.fhAmplitude), fhAmpId(fhAmpId),
+ fhNCellsPerCluster(qa.fhNCellsPerCluster), fhNCellsPerClusterMIP(qa.fhNCellsPerClusterMIP),fhNClusters(qa.fhNClusters),
+ fhCellTimeSpreadRespectToCellMax(qa.fhCellTimeSpreadRespectToCellMax),fhCellIdCellLargeTimeSpread(qa.fhCellIdCellLargeTimeSpread),
+ fhNCells(qa.fhNCells), fhAmplitude(qa.fhAmplitude), fhAmpId(fhAmpId),
fhTime(qa.fhTime), fhTimeId(qa.fhTimeId),fhTimeAmp(qa.fhTimeAmp),
fhCaloCorrNClusters(qa.fhCaloCorrNClusters), fhCaloCorrEClusters(qa.fhCaloCorrEClusters),
fhCaloCorrNCells(qa.fhCaloCorrNCells), fhCaloCorrECells(qa.fhCaloCorrECells),
- fhEMod(qa.fhEMod),fhNClustersMod(qa.fhNClustersMod), fhNCellsPerClusterMod(qa.fhNCellsPerClusterMod), fhNCellsMod(qa.fhNCellsMod),
- fhGridCellsMod(qa.fhGridCellsMod), fhGridCellsEMod(qa.fhGridCellsEMod), fhAmplitudeMod(qa.fhAmplitudeMod),
+ fhEMod(qa.fhEMod), fhNClustersMod(qa.fhNClustersMod), fhNCellsPerClusterMod(qa.fhNCellsPerClusterMod), fhNCellsMod(qa.fhNCellsMod),
+ fhGridCellsMod(qa.fhGridCellsMod), fhGridCellsEMod(qa.fhGridCellsEMod),
+ fhGridCellsTimeMod(qa.fhGridCellsTimeMod), fhAmplitudeMod(qa.fhAmplitudeMod), fhAmplitudeModFraction(qa.fhAmplitudeModFraction),
fhIMMod(qa.fhIMMod),fhIMCellCutMod(qa.fhIMCellCutMod),
fhGenGamPt(qa.fhGenGamPt), fhGenGamEta(qa.fhGenGamEta), fhGenGamPhi(qa.fhGenGamPhi),
fhGenPi0Pt(qa.fhGenPi0Pt), fhGenPi0Eta(qa.fhGenPi0Eta), fhGenPi0Phi(qa.fhGenPi0Phi),
fh2Phi = qa.fh2Phi;
fh2Eta = qa.fh2Eta;
+ fhCellTimeSpreadRespectToCellMax = qa.fhCellTimeSpreadRespectToCellMax;
+ fhCellIdCellLargeTimeSpread = qa.fhCellIdCellLargeTimeSpread;
+
fhNCells = qa.fhNCells;
fhAmplitude = qa.fhAmplitude;
fhAmpId = qa.fhAmpId;
fhEMod = qa.fhEMod; fhNClustersMod = qa.fhNClustersMod;
fhNCellsPerClusterMod = qa.fhNCellsPerClusterMod; fhNCellsMod = qa.fhNCellsMod;
- fhGridCellsMod = qa.fhGridCellsMod; fhGridCellsEMod = qa.fhGridCellsEMod;
- fhAmplitudeMod = qa.fhAmplitudeMod; fhIMMod=qa.fhIMMod; fhIMCellCutMod=qa.fhIMCellCutMod;
+ fhGridCellsMod = qa.fhGridCellsMod; fhGridCellsEMod = qa.fhGridCellsEMod; fhGridCellsTimeMod = qa.fhGridCellsTimeMod;
+ fhAmplitudeMod = qa.fhAmplitudeMod; fhAmplitudeModFraction = qa.fhAmplitudeModFraction;
+ fhIMMod=qa.fhIMMod; fhIMCellCutMod=qa.fhIMCellCutMod;
fhGenGamPt = qa.fhGenGamPt ; fhGenGamEta = qa.fhGenGamEta ; fhGenGamPhi = qa.fhGenGamPhi ;
fhGenPi0Pt = qa.fhGenPi0Pt ; fhGenPi0Eta = qa.fhGenPi0Eta ; fhGenPi0Phi = qa.fhGenPi0Phi ;
//If Geometry library loaded, do geometry selection during analysis.
if(fCalorimeter=="PHOS"){
- if(!GetReader()->GetPHOSGeometry()) printf("AliAnaPi0::GetCreateOutputObjects() - Initialize PHOS geometry!\n");
+ if(!GetReader()->GetPHOSGeometry()) printf("AliAnaCalorimeterQA::GetCreateOutputObjects() - Initialize PHOS geometry!\n");
GetReader()->InitPHOSGeometry();
}
else if(fCalorimeter=="EMCAL"){
- if(!GetReader()->GetEMCALGeometry()) printf("AliAnaPi0::GetCreateOutputObjects() - Initialize EMCAL geometry!\n");
+ if(!GetReader()->GetEMCALGeometry()) printf("AliAnaCalorimeterQA::GetCreateOutputObjects() - Initialize EMCAL geometry!\n");
GetReader()->InitEMCALGeometry();
}
fhAsym->SetYTitle("Asymmetry");
outputContainer->Add(fhAsym);
- fhNCellsPerCluster = new TH2F ("hNCellsPerCluster","# cells per cluster vs cluster energy", nptbins,ptmin,ptmax, nbins,nmin,nmax);
+ fhNCellsPerCluster = new TH2F ("hNCellsPerCluster","# cells per cluster vs cluster energy",nptbins,ptmin,ptmax, nbins,nmin,nmax);
fhNCellsPerCluster->SetXTitle("E (GeV)");
fhNCellsPerCluster->SetYTitle("n cells");
outputContainer->Add(fhNCellsPerCluster);
fhAmplitude->SetXTitle("Cell Energy (GeV)");
outputContainer->Add(fhAmplitude);
- fhAmpId = new TH2F ("hAmpId","Cell Energy", nptbins*2,ptmin,ptmax,rowmax*colmax,0,rowmax*colmax);
+ fhAmpId = new TH2F ("hAmpId","Cell Energy", nptbins*4,ptmin,ptmax*2,rowmax*colmax*fNModules,0,rowmax*colmax*fNModules);
fhAmpId->SetXTitle("Cell Energy (GeV)");
outputContainer->Add(fhAmpId);
//Cell Time histograms, time only available in ESDs
if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
+
+ fhCellTimeSpreadRespectToCellMax = new TH1F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", 100,-200,200);
+ fhCellTimeSpreadRespectToCellMax->SetXTitle("#Delta t (ns)");
+ outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
+
+ fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","", colmax*rowmax*fNModules,0,colmax*rowmax*fNModules);
+ fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
+ outputContainer->Add(fhCellIdCellLargeTimeSpread);
+
fhTime = new TH1F ("hTime","Cell Time",ntimebins,timemin,timemax);
fhTime->SetXTitle("Cell Time (ns)");
outputContainer->Add(fhTime);
- fhTimeId = new TH2F ("hTimeId","Cell Time vs Absolute Id",ntimebins,timemin,timemax,rowmax*colmax,0,rowmax*colmax);
+ fhTimeId = new TH2F ("hTimeId","Cell Time vs Absolute Id",ntimebins,timemin,timemax,rowmax*colmax*fNModules,0,rowmax*colmax*fNModules);
fhTimeId->SetXTitle("Cell Time (ns)");
fhTimeId->SetYTitle("Cell Absolute Id");
outputContainer->Add(fhTimeId);
}//correlate calorimeters
//Module histograms
- fhEMod = new TH1F*[fNModules];
- fhNClustersMod = new TH1F*[fNModules];
- fhNCellsPerClusterMod = new TH2F*[fNModules];
- fhNCellsMod = new TH1F*[fNModules];
- fhGridCellsMod = new TH2F*[fNModules];
- fhGridCellsEMod = new TH2F*[fNModules];
- fhAmplitudeMod = new TH1F*[fNModules];
- fhIMMod = new TH2F*[fNModules];
- fhIMCellCutMod = new TH2F*[fNModules];
+ fhEMod = new TH1F*[fNModules];
+ fhNClustersMod = new TH1F*[fNModules];
+ fhNCellsPerClusterMod = new TH2F*[fNModules];
+ fhNCellsMod = new TH1F*[fNModules];
+ fhGridCellsMod = new TH2F*[fNModules];
+ fhGridCellsEMod = new TH2F*[fNModules];
+ fhGridCellsTimeMod = new TH2F*[fNModules];
+ fhAmplitudeMod = new TH1F*[fNModules];
+ if(fCalorimeter=="EMCAL")
+ fhAmplitudeModFraction = new TH1F*[fNModules*3];
+ fhIMMod = new TH2F*[fNModules];
+ fhIMCellCutMod = new TH2F*[fNModules];
for(Int_t imod = 0; imod < fNModules; imod++){
fhGridCellsEMod[imod]->SetXTitle("column (eta direction)");
outputContainer->Add(fhGridCellsEMod[imod]);
+ fhGridCellsTimeMod[imod] = new TH2F (Form("hGridCellsTime_Mod%d",imod),Form("Accumulated time in grid of cells in Module %d, with E > 0.5 GeV",imod),
+ colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
+ fhGridCellsTimeMod[imod]->SetYTitle("row (phi direction)");
+ fhGridCellsTimeMod[imod]->SetXTitle("column (eta direction)");
+ outputContainer->Add(fhGridCellsTimeMod[imod]);
+
fhAmplitudeMod[imod] = new TH1F (Form("hAmplitude_Mod%d",imod),Form("Cell Energy in Module %d",imod), nptbins*2,ptmin,ptmax);
fhAmplitudeMod[imod]->SetXTitle("Cell Energy (GeV)");
outputContainer->Add(fhAmplitudeMod[imod]);
+ if(fCalorimeter == "EMCAL"){
+ for(Int_t ifrac = 0; ifrac < 3; ifrac++){
+ fhAmplitudeModFraction[imod*3+ifrac] = new TH1F (Form("hAmplitude_Mod%d_Frac%d",imod,ifrac),Form("Cell reconstructed Energy in Module %d, Fraction %d ",imod,ifrac), nptbins,ptmin,ptmax);
+ fhAmplitudeModFraction[imod*3+ifrac]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhAmplitudeModFraction[imod*3+ifrac]);
+ }
+ }
+
fhIMMod[imod] = new TH2F (Form("hIM_Mod%d",imod),
Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d",imod),
nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
GetReader()->GetVertex(v);
TObject * track = 0x0;
//Loop over CaloClusters
+ //if(nCaloClusters > 0)printf("QA : Vertex Cut passed %f, cut %f, entries %d, %s\n",v[2], 40., nCaloClusters, fCalorimeter.Data());
for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
if(GetDebug() > 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
nLabel = clus->GetNLabels();
if(clus->GetLabels()) labels = (clus->GetLabels())->GetArray();
//Cells per cluster
- nCaloCellsPerCluster = clus->GetNCells();
+ nCaloCellsPerCluster = clus->GetNCells();
+ //if(mom.E() > 10 && nCaloCellsPerCluster == 1 ) printf("%s:************** E = %f ********** ncells = %d\n",fCalorimeter.Data(), mom.E(),nCaloCellsPerCluster);
//matched cluster with tracks
nTracksMatched = clus->GetNTracksMatched();
trackIndex = clus->GetTrackMatched();
if(nTracksMatched == 1) nTracksMatched = 0;
track = 0;
}
+
+
+ //Cells in cluster, check time of cells respect to max energy cell
+ if(nCaloCellsPerCluster > 1){
+ AliESDCaloCells * cell = 0x0;
+ if(fCalorimeter == "PHOS") cell = ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
+ else cell = ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
+ //Get list of contributors
+ UShort_t * indexList = clus->GetCellsAbsId() ;
+ //Get maximum energy cell
+ Float_t emax = -1;
+ Double_t tmax = -1;
+ Int_t imax = -1;
+ Int_t absId = -1 ;
+ //printf("nCaloCellsPerCluster %d\n",nCaloCellsPerCluster);
+ for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
+ // printf("Index %d\n",ipos);
+ absId = indexList[ipos];
+ if(cell->GetCellAmplitude(absId) > emax) {
+ imax = ipos;
+ emax = cell->GetCellAmplitude(absId);
+ tmax = cell->GetCellTime(absId);
+ }
+ }// max search loop
+
+ for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
+ if(imax == ipos) continue;
+ absId = indexList[ipos];
+ Float_t diff = (tmax-cell->GetCellTime(absId))*1e9;
+ fhCellTimeSpreadRespectToCellMax->Fill(diff);
+ if(TMath::Abs(TMath::Abs(diff) > 100)) fhCellIdCellLargeTimeSpread->Fill(absId);
+ }// fill cell-cluster histogram loop
+
+
+ }//check time of cells respect to max energy cell
+
+
}
else{
AliAODCaloCluster* clus = (AliAODCaloCluster*) (caloClusters->At(iclus));
ClusterHistograms(mom, nCaloCellsPerCluster, nModule, nTracksMatched, track, labels, nLabel);
if(GetDebug()>1) printf("Invariant mass \n");
//Invariant mass
+
+ //do not do for bad vertex
+ Float_t fZvtxCut = 40. ;
+ if(v[2]<-fZvtxCut || v[2]> fZvtxCut) continue ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
+
Int_t nModule2 = -1;
Int_t nCaloCellsPerCluster2=0;
if (nCaloClusters > 1 ) {
- for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters-1 ; jclus++) {
+ for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++) {
if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
AliESDCaloCluster* clus2 = (AliESDCaloCluster*) (caloClusters->At(jclus));
//Get cluster kinematics
//Fill invariant mass histograms
//All modules
+
+ //printf("QA : Fill inv mass histo: pt1 %f, pt2 %f, pt12 %f, mass %f, calo %s \n",mom.Pt(),mom2.Pt(),(mom+mom2).Pt(),(mom+mom2).M(), fCalorimeter.Data());
fhIM ->Fill((mom+mom2).Pt(),(mom+mom2).M());
//Single module
if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
else cell = ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
if(!cell) {
- printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No CELLS available for analysis");
+ printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No %s ESD CELLS available for analysis\n",fCalorimeter.Data());
abort();
}
fhTimeId ->Fill(time,id);
fhTimeAmp ->Fill(amp,time);
fhAmplitudeMod[nModule]->Fill(cell->GetAmplitude(iCell));
+ if(fCalorimeter=="EMCAL"){
+ Int_t ifrac = 0;
+ if(icol > 15 && icol < 32) ifrac = 1;
+ else if(icol > 31) ifrac = 2;
+ fhAmplitudeModFraction[nModule*3+ifrac]->Fill(cell->GetAmplitude(iCell));
+ }
nCellsInModule[nModule]++;
- fhGridCellsMod[nModule] ->Fill(icol,irow);
- fhGridCellsEMod[nModule]->Fill(icol,irow,amp);
+ fhGridCellsMod[nModule] ->Fill(icol,irow);
+ fhGridCellsEMod[nModule] ->Fill(icol,irow,amp);
+ if(amp > 0.5)fhGridCellsTimeMod[nModule]->Fill(icol,irow,time);
+
}//nmodules
}//cell loop
}//ESD
else cell = ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
if(!cell) {
- printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No CELLS available for analysis");
- abort();
+ printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No %s AOD CELLS available for analysis\n",fCalorimeter.Data());
+ //abort();
+ return;
}
ncells = cell->GetNumberOfCells() ;
fhNCells->Fill(ncells) ;
if(GetDebug() > 0)
- printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In ESD %s cell entries %d\n", fCalorimeter.Data(), ncells);
+ printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In AOD %s cell entries %d\n", fCalorimeter.Data(), ncells);
for (Int_t iCell = 0; iCell < ncells; iCell++) {
if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell));
fhAmplitude->Fill(amp);
fhAmpId ->Fill(amp,id);
fhAmplitudeMod[nModule]->Fill(cell->GetAmplitude(iCell));
+ if(fCalorimeter=="EMCAL"){
+ Int_t ifrac = 0;
+ if(icol > 15 && icol < 32) ifrac = 1;
+ else if(icol > 31) ifrac = 2;
+ fhAmplitudeModFraction[nModule*3+ifrac]->Fill(cell->GetAmplitude(iCell));
+ }
nCellsInModule[nModule]++;
- fhGridCellsMod[nModule]->Fill(icol,irow);
+ fhGridCellsMod[nModule] ->Fill(icol,irow);
fhGridCellsEMod[nModule]->Fill(icol,irow,amp);
}//nmodules
}//cell loop
fhNCellsPerCluster ->Fill(e, nCaloCellsPerCluster);
fhNCellsPerClusterMIP->Fill(e, nCaloCellsPerCluster);
if(nModule < fNModules) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
-
+
//Fill histograms only possible when simulation
if(IsDataMC() && nLabels > 0 && labels){
fhNCellsPerCluster = (TH2F *) outputList->At(index++);
fhNCellsPerClusterMIP = (TH2F *) outputList->At(index++);
fhNClusters = (TH1F *) outputList->At(index++);
+
fhNCells = (TH1F *) outputList->At(index++);
fhAmplitude = (TH1F *) outputList->At(index++);
fhAmpId = (TH2F *) outputList->At(index++);
if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
+
+ fhCellTimeSpreadRespectToCellMax = (TH1F *) outputList->At(index++);
+ fhCellIdCellLargeTimeSpread = (TH1F *) outputList->At(index++);
+
fhTime = (TH1F *) outputList->At(index++);
fhTimeId = (TH2F *) outputList->At(index++);
fhTimeAmp = (TH2F *) outputList->At(index++);
}
//Module histograms
- fhEMod = new TH1F*[fNModules];
- fhNClustersMod = new TH1F*[fNModules];
- fhNCellsPerClusterMod = new TH2F*[fNModules];
- fhNCellsMod = new TH1F*[fNModules];
- fhGridCellsMod = new TH2F*[fNModules];
- fhGridCellsEMod = new TH2F*[fNModules];
- fhAmplitudeMod = new TH1F*[fNModules];
- fhIMMod = new TH2F*[fNModules];
- fhIMCellCutMod = new TH2F*[fNModules];
+ fhEMod = new TH1F*[fNModules];
+ fhNClustersMod = new TH1F*[fNModules];
+ fhNCellsPerClusterMod = new TH2F*[fNModules];
+ fhNCellsMod = new TH1F*[fNModules];
+ fhGridCellsMod = new TH2F*[fNModules];
+ fhGridCellsEMod = new TH2F*[fNModules];
+ if(GetReader()->GetDataType()==AliCaloTrackReader::kESD)
+ fhGridCellsTimeMod = new TH2F*[fNModules];
+ fhAmplitudeMod = new TH1F*[fNModules];
+ if(fCalorimeter=="EMCAL")
+ fhAmplitudeModFraction = new TH1F*[fNModules*3];
+ fhIMMod = new TH2F*[fNModules];
+ fhIMCellCutMod = new TH2F*[fNModules];
for(Int_t imod = 0 ; imod < fNModules; imod++){
- fhEMod[imod] = (TH1F *) outputList->At(index++);
- fhNClustersMod[imod] = (TH1F *) outputList->At(index++);
- fhNCellsPerClusterMod[imod] = (TH2F *) outputList->At(index++);
- fhNCellsMod[imod] = (TH1F *) outputList->At(index++);
- fhGridCellsMod[imod] = (TH2F *) outputList->At(index++);
- fhGridCellsEMod[imod] = (TH2F *) outputList->At(index++);
- fhAmplitudeMod[imod] = (TH1F *) outputList->At(index++);
- fhIMMod[imod] = (TH2F *) outputList->At(index++);
- fhIMCellCutMod[imod] = (TH2F *) outputList->At(index++);
+ fhEMod[imod] = (TH1F *) outputList->At(index++);
+ fhNClustersMod[imod] = (TH1F *) outputList->At(index++);
+ fhNCellsPerClusterMod[imod] = (TH2F *) outputList->At(index++);
+ fhNCellsMod[imod] = (TH1F *) outputList->At(index++);
+ fhGridCellsMod[imod] = (TH2F *) outputList->At(index++);
+ fhGridCellsEMod[imod] = (TH2F *) outputList->At(index++);
+ if(GetReader()->GetDataType()==AliCaloTrackReader::kESD)
+ fhGridCellsTimeMod[imod] = (TH2F *) outputList->At(index++);
+ fhAmplitudeMod[imod] = (TH1F *) outputList->At(index++);
+ if(fCalorimeter=="EMCAL"){
+ for(Int_t ifrac = 0; ifrac < 3; ifrac++){
+ fhAmplitudeModFraction[imod*3+ifrac] = (TH1F *) outputList->At(index++);
+ }
+ }
+ fhIMMod[imod] = (TH2F *) outputList->At(index++);
+ fhIMCellCutMod[imod] = (TH2F *) outputList->At(index++);
}
TH2F * fhNCellsPerClusterMIP; //! N cells per cluster, finer fixed pT bin for MIP search.
TH1F * fhNClusters; //! Number of clusters
+ TH1F * fhCellTimeSpreadRespectToCellMax; //! Difference of the time of cell with maximum dep energy and the rest of cells
+ TH1F * fhCellIdCellLargeTimeSpread; //! Cells with large time respect to max (diff > 100 ns)
+
//Calo Cells
TH1F * fhNCells; //! Number of towers/crystals with signal
TH1F * fhAmplitude; //! Amplitude measured in towers/crystals
TH1F ** fhNCellsMod ; //! Number of towers/crystals with signal different module, Reco
TH2F ** fhGridCellsMod ; //! Cells ordered in column/row for different module, Reco
TH2F ** fhGridCellsEMod ; //! Cells ordered in column/row for different module, weighted with energy, Reco
+ TH2F ** fhGridCellsTimeMod ; //! Cells ordered in column/row for different module, weighted with time, Reco
TH1F ** fhAmplitudeMod ; //! Amplitude measured in towers/crystals different module, Reco
+ TH1F ** fhAmplitudeModFraction; //! Amplitude measured in towers/crystals different fractions of module, Reco
TH2F ** fhIMMod; //! cluster pairs invariant mass, different module,
TH2F ** fhIMCellCutMod; //! cluster pairs invariant mass, n cells > 1 per cluster, different module
TH2F *fhMCChHad1pOverER02; //! p/E for track-cluster matches, dR > 0.2, MC charged hadrons
TH2F *fhMCNeutral1pOverER02; //! p/E for track-cluster matches, dR > 0.2, MC neutral
- ClassDef(AliAnaCalorimeterQA,5)
+ ClassDef(AliAnaCalorimeterQA,6)
} ;
// --- ROOT system ---
#include "TH3.h"
+#include "TH2D.h"
//#include "Riostream.h"
#include "TCanvas.h"
#include "TPad.h"
#include "AliESDCaloCluster.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
+#include "AliNeutralMesonSelection.h"
ClassImp(AliAnaPi0)
AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
-fNModules(12), fEventsList(0x0), //fhEtalon(0x0),
+fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), //fhEtalon(0x0),
fhReMod(0x0), fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
-fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0)
+fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
+fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
+fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
{
//Default Ctor
InitParameters();
AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),
fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter),
-fNModules(ex.fNModules), fEventsList(ex.fEventsList), //fhEtalon(ex.fhEtalon),
+fNModules(ex.fNModules), fUseAngleCut(ex.fUseAngleCut), fEventsList(ex.fEventsList), //fhEtalon(ex.fhEtalon),
fhReMod(ex.fhReMod), fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),
fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
+fhRealOpeningAngle(ex.fhRealOpeningAngle),fhRealCosOpeningAngle(ex.fhRealCosOpeningAngle),
fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY),
-fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi)
+fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi),
+fhPrimOpeningAngle(ex.fhPrimOpeningAngle),fhPrimCosOpeningAngle(ex.fhPrimCosOpeningAngle)
{
// cpy ctor
//Do not need it
fNCentrBin = ex.fNCentrBin ; fNZvertBin = ex.fNZvertBin ; fNrpBin = ex.fNrpBin ;
fNPID = ex.fNPID ; fNmaxMixEv = ex.fNmaxMixEv ; fZvtxCut = ex.fZvtxCut ; fCalorimeter = ex.fCalorimeter ;
- fNModules = ex.fNModules; fEventsList = ex.fEventsList ; //fhEtalon = ex.fhEtalon ;
+ fNModules = ex.fNModules; fEventsList = ex.fEventsList ; //fhEtalon = ex.fhEtalon ;
fhRe1 = ex.fhRe1 ; fhMi1 = ex.fhMi1 ; fhRe2 = ex.fhRe2 ; fhMi2 = ex.fhMi2 ; fhReMod = ex.fhReMod;
- fhRe3 = ex.fhRe3 ; fhMi3 = ex.fhMi3 ; fhEvents = ex.fhEvents ;
+ fhRe3 = ex.fhRe3 ; fhMi3 = ex.fhMi3 ; fhEvents = ex.fhEvents ; fUseAngleCut = ex.fUseAngleCut;
fhPrimPt = ex.fhPrimPt ; fhPrimAccPt = ex.fhPrimAccPt ; fhPrimY = ex.fhPrimY ;
fhPrimAccY = ex.fhPrimAccY ; fhPrimPhi = ex.fhPrimPhi ; fhPrimAccPhi = ex.fhPrimAccPhi ;
-
+ fhRealOpeningAngle = ex.fhRealOpeningAngle; fhRealCosOpeningAngle = ex.fhRealCosOpeningAngle;
+ fhPrimOpeningAngle = ex.fhPrimOpeningAngle; fhPrimCosOpeningAngle = ex.fhPrimCosOpeningAngle;
+
return *this;
}
fNmaxMixEv = 10;
fZvtxCut = 40;
fCalorimeter = "PHOS";
+ fUseAngleCut = kFALSE;
+
}
//________________________________________________________________________________________________________________________________________________
//void AliAnaPi0::Init()
fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
outputContainer->Add(fhEvents) ;
+
+ fhRealOpeningAngle = new TH2D
+ ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5);
+ fhRealOpeningAngle->SetYTitle("#theta(rad)");
+ fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
+ outputContainer->Add(fhRealOpeningAngle) ;
+
+ fhRealCosOpeningAngle = new TH2D
+ ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1);
+ fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
+ fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
+ outputContainer->Add(fhRealCosOpeningAngle) ;
+
+
//Histograms filled only if MC data is requested
if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
// if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
outputContainer->Add(fhPrimAccPhi) ;
+
+
+ fhPrimOpeningAngle = new TH2D
+ ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5);
+ fhPrimOpeningAngle->SetYTitle("#theta(rad)");
+ fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
+ outputContainer->Add(fhPrimOpeningAngle) ;
+
+ fhPrimCosOpeningAngle = new TH2D
+ ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1);
+ fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
+ fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
+ outputContainer->Add(fhPrimCosOpeningAngle) ;
+
}
-
+
for(Int_t imod=0; imod<fNModules; imod++){
- //Module dependent invariant mass
- sprintf(key,"hReMod_%d",imod) ;
- sprintf(title,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
- //fhEtalon->Clone(key);
- //fhReMod[imod]=(TH3D*)fhEtalon->Clone(key) ;
- //fhReMod[imod]->SetName(key) ;
- //fhReMod[imod]->SetTitle(title) ;
- fhReMod[imod] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
- outputContainer->Add(fhReMod[imod]) ;
+ //Module dependent invariant mass
+ sprintf(key,"hReMod_%d",imod) ;
+ sprintf(title,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
+ //fhEtalon->Clone(key);
+ //fhReMod[imod]=(TH3D*)fhEtalon->Clone(key) ;
+ //fhReMod[imod]->SetName(key) ;
+ //fhReMod[imod]->SetTitle(title) ;
+ fhReMod[imod] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
+ outputContainer->Add(fhReMod[imod]) ;
}
-
-
+
+
//Save parameters used for analysis
TString parList ; //this will be list of parameters used for this analysis.
char onePar[255] ;
printf("Number of different PID used: %d \n",fNPID) ;
printf("Cuts: \n") ;
printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
- printf("Number of modules: %d \n",fNModules) ;
+ printf("Number of modules: %d \n",fNModules) ;
+ printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
printf("------------------------------------------------------\n") ;
}
//____________________________________________________________________________________________________________________________________________________
Int_t AliAnaPi0::GetModuleNumber(AliAODPWG4Particle * particle)
{
- //Get the EMCAL/PHOS module number that corresponds to this particle
-
- Int_t absId = -1;
- if(fCalorimeter=="EMCAL"){
- GetReader()->GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
- if(GetDebug() > 2)
- printf("AliAnaPi0::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
- particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId));
- return GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
- }//EMCAL
- else{//PHOS
- Int_t relId[4];
- if(!strcmp((GetReader()->GetInputEvent())->GetName(),"AliESDEvent")) {
- AliESDCaloCluster *cluster = ((AliESDEvent*)GetReader()->GetInputEvent())->GetCaloCluster(particle->GetCaloLabel(0));
- if ( cluster->GetNCells() > 0) {
- absId = cluster->GetCellAbsId(0);
- if(GetDebug() > 2)
- printf("AliAnaPi0::GetModuleNumber(ESD) - PHOS: cluster eta %f, phi %f, e %f, e cluster %f, absId %d\n",
- particle->Eta(), particle->Phi()*TMath::RadToDeg(), particle->E(), cluster->E(), absId);
- }
- else return -1;
- }//ESDs
- else{
- AliAODCaloCluster *cluster = ((AliAODEvent*)GetReader()->GetInputEvent())->GetCaloCluster(particle->GetCaloLabel(0));
- if ( cluster->GetNCells() > 0) {
- absId = cluster->GetCellAbsId(0);
- if(GetDebug() > 2)
- printf("AliAnaPi0::GetModuleNumber(AOD) - PHOS: cluster eta %f, phi %f, e %f, e cluster %f, absId %d\n",
- particle->Eta(), particle->Phi()*TMath::RadToDeg(), particle->E(), cluster->E(), absId);
- }
- else return -1;
- }//AODs
-
- if ( absId >= 0) {
- GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
- if(GetDebug() > 2)
- printf("PHOS: Module %d\n",relId[0]-1);
- return relId[0]-1;
- }
- else return -1;
- }//PHOS
-
- return -1;
+ //Get the EMCAL/PHOS module number that corresponds to this particle
+
+ Int_t absId = -1;
+ if(fCalorimeter=="EMCAL"){
+ GetReader()->GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
+ if(GetDebug() > 2)
+ printf("AliAnaPi0::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
+ particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId));
+ return GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
+ }//EMCAL
+ else{//PHOS
+ Int_t relId[4];
+ if(!strcmp((GetReader()->GetInputEvent())->GetName(),"AliESDEvent")) {
+ AliESDCaloCluster *cluster = ((AliESDEvent*)GetReader()->GetInputEvent())->GetCaloCluster(particle->GetCaloLabel(0));
+ if ( cluster->GetNCells() > 0) {
+ absId = cluster->GetCellAbsId(0);
+ if(GetDebug() > 2)
+ printf("AliAnaPi0::GetModuleNumber(ESD) - PHOS: cluster eta %f, phi %f, e %f, e cluster %f, absId %d\n",
+ particle->Eta(), particle->Phi()*TMath::RadToDeg(), particle->E(), cluster->E(), absId);
+ }
+ else return -1;
+ }//ESDs
+ else{
+ AliAODCaloCluster *cluster = ((AliAODEvent*)GetReader()->GetInputEvent())->GetCaloCluster(particle->GetCaloLabel(0));
+ if ( cluster->GetNCells() > 0) {
+ absId = cluster->GetCellAbsId(0);
+ if(GetDebug() > 2)
+ printf("AliAnaPi0::GetModuleNumber(AOD) - PHOS: cluster eta %f, phi %f, e %f, e cluster %f, absId %d\n",
+ particle->Eta(), particle->Phi()*TMath::RadToDeg(), particle->E(), cluster->E(), absId);
+ }
+ else return -1;
+ }//AODs
+
+ if ( absId >= 0) {
+ GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
+ if(GetDebug() > 2)
+ printf("PHOS: Module %d\n",relId[0]-1);
+ return relId[0]-1;
+ }
+ else return -1;
+ }//PHOS
+
+ return -1;
}
//____________________________________________________________________________________________________________________________________________________
if(GetDebug() > 2)
printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
p1->Pt(), p2->Pt(), pt,m,a);
+
+ //Check if opening angle is too large or too small compared to what is expected
+ Double_t angle = photon1.Angle(photon2.Vect());
+ //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
+ //printf("angle %f\n",angle);
+ if(fUseAngleCut && angle < 0.1) continue;
+ fhRealOpeningAngle ->Fill(pt,angle);
+ fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
+
//Fill module dependent histograms
//if(module1==module2) printf("mod1 %d\n",module1);
if(module1==module2 && module1 >=0 && module1<fNModules)
m = (photon1+photon2).M() ;
Double_t pt = (photon1 + photon2).Pt();
Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
+
+ //Check if opening angle is too large or too small compared to what is expected
+ Double_t angle = photon1.Angle(photon2.Vect());
+ //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
+ if(fUseAngleCut && angle < 0.1) continue;
+
if(GetDebug() > 2)
printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
p1->Pt(), p2->Pt(), pt,m,a);
//Acceptance
if(IsDataMC() && GetReader()->ReadStack()){
- AliStack * stack = GetMCStack();
- if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
- for(Int_t i=0 ; i<stack->GetNprimary(); i++){
- TParticle * prim = stack->Particle(i) ;
- if(prim->GetPdgCode() == 111){
- Double_t pi0Pt = prim->Pt() ;
- //printf("pi0, pt %2.2f\n",pi0Pt);
- if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
- Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
- Double_t phi = TMath::RadToDeg()*prim->Phi() ;
- if(TMath::Abs(pi0Y) < 0.5){
- fhPrimPt->Fill(pi0Pt) ;
- }
- fhPrimY ->Fill(pi0Y) ;
- fhPrimPhi->Fill(phi) ;
-
- //Check if both photons hit Calorimeter
- Int_t iphot1=prim->GetFirstDaughter() ;
- Int_t iphot2=prim->GetLastDaughter() ;
- if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
- TParticle * phot1 = stack->Particle(iphot1) ;
- TParticle * phot2 = stack->Particle(iphot2) ;
- if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
- //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());
- Bool_t inacceptance = kFALSE;
-
- if(fCalorimeter == "PHOS"){
- if(GetReader()->GetPHOSGeometry()){
- Int_t mod ;
- Double_t x,z ;
- if(GetReader()->GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetReader()->GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
- inacceptance = kTRUE;
- if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
- }
- else{
- TLorentzVector lv1, lv2;
- phot1->Momentum(lv1);
- phot2->Momentum(lv2);
- if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
- inacceptance = kTRUE ;
- if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
- }
-
- }
- else if(fCalorimeter == "EMCAL"){
- if(GetReader()->GetEMCALGeometry()){
- if(GetReader()->GetEMCALGeometry()->Impact(phot1) && GetReader()->GetEMCALGeometry()->Impact(phot2))
- inacceptance = kTRUE;
- if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
- }
- else{
- TLorentzVector lv1, lv2;
- phot1->Momentum(lv1);
- phot2->Momentum(lv2);
- if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
- inacceptance = kTRUE ;
- if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
- }
- }
+ AliStack * stack = GetMCStack();
+ if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
+ for(Int_t i=0 ; i<stack->GetNprimary(); i++){
+ TParticle * prim = stack->Particle(i) ;
+ if(prim->GetPdgCode() == 111){
+ Double_t pi0Pt = prim->Pt() ;
+ //printf("pi0, pt %2.2f\n",pi0Pt);
+ if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
+ Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
+ Double_t phi = TMath::RadToDeg()*prim->Phi() ;
+ if(TMath::Abs(pi0Y) < 0.5){
+ fhPrimPt->Fill(pi0Pt) ;
+ }
+ fhPrimY ->Fill(pi0Y) ;
+ fhPrimPhi->Fill(phi) ;
+
+ //Check if both photons hit Calorimeter
+ Int_t iphot1=prim->GetFirstDaughter() ;
+ Int_t iphot2=prim->GetLastDaughter() ;
+ if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
+ TParticle * phot1 = stack->Particle(iphot1) ;
+ TParticle * phot2 = stack->Particle(iphot2) ;
+ if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
+ //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;
+ phot1->Momentum(lv1);
+ phot2->Momentum(lv2);
+
+ Bool_t inacceptance = kFALSE;
+ if(fCalorimeter == "PHOS"){
+ if(GetReader()->GetPHOSGeometry()){
+ Int_t mod ;
+ Double_t x,z ;
+ if(GetReader()->GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetReader()->GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
+ inacceptance = kTRUE;
+ if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ else{
+
+ if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
+ inacceptance = kTRUE ;
+ if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+
+ }
+ else if(fCalorimeter == "EMCAL"){
+ if(GetReader()->GetEMCALGeometry()){
+ if(GetReader()->GetEMCALGeometry()->Impact(phot1) && GetReader()->GetEMCALGeometry()->Impact(phot2))
+ inacceptance = kTRUE;
+ if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ else{
+ if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
+ inacceptance = kTRUE ;
+ if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ }
+
+ if(inacceptance){
- if(inacceptance){
- fhPrimAccPt->Fill(pi0Pt) ;
- fhPrimAccPhi->Fill(phi) ;
- fhPrimAccY->Fill(pi0Y) ;
- }//Accepted
- }// 2 photons
- }//Check daughters exist
- }// Primary pi0
- }//loop on primaries
- }//stack exists and data is MC
+ fhPrimAccPt->Fill(pi0Pt) ;
+ fhPrimAccPhi->Fill(phi) ;
+ fhPrimAccY->Fill(pi0Y) ;
+ Double_t angle = lv1.Angle(lv2.Vect());
+ fhPrimOpeningAngle ->Fill(pi0Pt,angle);
+ fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
+
+ }//Accepted
+ }// 2 photons
+ }//Check daughters exist
+ }// Primary pi0
+ }//loop on primaries
+ }//stack exists and data is MC
}//read stack
else if(GetReader()->ReadAODMCParticles()){
- if(GetDebug() >= 0) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
+ if(GetDebug() >= 0) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
}
}
//________________________________________________________________________
void AliAnaPi0::ReadHistograms(TList* outputList)
{
- // Needed when Terminate is executed in distributed environment
- // Refill analysis histograms of this class with corresponding histograms in output list.
-
- // Histograms of this analsys are kept in the same list as other analysis, recover the position of
- // the first one and then add the next.
- Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
-
- if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
- if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
- if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
- if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
- if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
- if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
- if(!fhReMod) fhReMod = new TH3D*[fNModules] ;
-
- for(Int_t ic=0; ic<fNCentrBin; ic++){
- for(Int_t ipid=0; ipid<fNPID; ipid++){
- fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
- fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
- fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
- fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
- fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
- fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
- }
- }
-
- fhEvents = (TH3D *) outputList->At(index++);
-
- //Histograms filled only if MC data is requested
- if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
- fhPrimPt = (TH1D*) outputList->At(index++);
- fhPrimAccPt = (TH1D*) outputList->At(index++);
- fhPrimY = (TH1D*) outputList->At(index++);
- fhPrimAccY = (TH1D*) outputList->At(index++);
- fhPrimPhi = (TH1D*) outputList->At(index++);
- fhPrimAccPhi = (TH1D*) outputList->At(index++);
- }
-
- for(Int_t imod=0; imod < fNModules; imod++)
- fhReMod[imod] = (TH3D*) outputList->At(index++);
-
+ // Needed when Terminate is executed in distributed environment
+ // Refill analysis histograms of this class with corresponding histograms in output list.
+
+ // Histograms of this analsys are kept in the same list as other analysis, recover the position of
+ // the first one and then add the next.
+ Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
+
+ if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
+ if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
+ if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
+ if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
+ if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
+ if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
+ if(!fhReMod) fhReMod = new TH3D*[fNModules] ;
+
+ for(Int_t ic=0; ic<fNCentrBin; ic++){
+ for(Int_t ipid=0; ipid<fNPID; ipid++){
+ fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
+ fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
+ fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
+ fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
+ fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
+ fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
+ }
+ }
+
+ fhEvents = (TH3D *) outputList->At(index++);
+
+ //Histograms filled only if MC data is requested
+ if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
+ fhPrimPt = (TH1D*) outputList->At(index++);
+ fhPrimAccPt = (TH1D*) outputList->At(index++);
+ fhPrimY = (TH1D*) outputList->At(index++);
+ fhPrimAccY = (TH1D*) outputList->At(index++);
+ fhPrimPhi = (TH1D*) outputList->At(index++);
+ fhPrimAccPhi = (TH1D*) outputList->At(index++);
+ }
+
+ for(Int_t imod=0; imod < fNModules; imod++)
+ fhReMod[imod] = (TH3D*) outputList->At(index++);
+
}
//Do some calculations and plots from the final histograms.
printf(" *** %s Terminate:\n", GetName()) ;
-
+
//Recover histograms from output histograms list, needed for distributed analysis.
ReadHistograms(outputList);
-
+
if (!fhRe1) {
- printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
- return;
+ printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
+ return;
}
-
+
printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
-
+
char nameIM[128];
sprintf(nameIM,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
cIM->Divide(2, 2);
-
+
cIM->cd(1) ;
//gPad->SetLogy();
TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPtAll_%s",fCalorimeter.Data()));
//Root
class TList;
class TH3D ;
+class TH2D ;
//Analysis
class AliAODEvent ;
Int_t GetNPID() const {return fNPID ; }
void SetNPID(Int_t n) {fNPID = n ; }
+ void SwitchOnAngleSelection() {fUseAngleCut = kTRUE ; }
+ void SwitchOffAngleSelection() {fUseAngleCut = kFALSE ; }
+
private:
Bool_t IsBadRun(Int_t /*iRun*/) const {return kFALSE;} //Tests if this run bad according to private list
private:
- Int_t fNCentrBin ; // Number of bins in event container for centrality
- Int_t fNZvertBin ; // Number of bins in event container for vertex position
- Int_t fNrpBin ; // Number of bins in event container for reaction plain
- Int_t fNPID ; // Number of possible PID combinations
- Int_t fNmaxMixEv ; // Maximal number of events stored in buffer for mixing
- Float_t fZvtxCut ; // Cut on vertex position
- TString fCalorimeter ; // Select Calorimeter for IM
- Int_t fNModules ; // Number of EMCAL/PHOS modules, set as many histogras as modules
-
- TList ** fEventsList ; //! containers for photons in stored events
+ Int_t fNCentrBin ; // Number of bins in event container for centrality
+ Int_t fNZvertBin ; // Number of bins in event container for vertex position
+ Int_t fNrpBin ; // Number of bins in event container for reaction plain
+ Int_t fNPID ; // Number of possible PID combinations
+ Int_t fNmaxMixEv ; // Maximal number of events stored in buffer for mixing
+ Float_t fZvtxCut ; // Cut on vertex position
+ TString fCalorimeter ; // Select Calorimeter for IM
+ Int_t fNModules ; // Number of EMCAL/PHOS modules, set as many histogras as modules
+ Bool_t fUseAngleCut ; // Select pairs depending on their opening angle
+ TList ** fEventsList ; //! Containers for photons in stored events
//Histograms
TH3D ** fhMi2 ; //!MIXED two-photon invariant mass distribution for different centralities and PID
TH3D ** fhRe3 ; //!REAL two-photon invariant mass distribution for different centralities and PID
TH3D ** fhMi3 ; //!MIXED two-photon invariant mass distribution for different centralities and PID
- TH3D * fhEvents; //!Number of events per centrality, RP, zbin
+ TH3D * fhEvents; //!Number of events per centrality, RP, zbin
+
+ TH2D * fhRealOpeningAngle ; //! Opening angle of pair versus pair energy
+ TH2D * fhRealCosOpeningAngle ; //! Cosinus of opening angle of pair version pair energy
//Acceptance
TH1D * fhPrimPt ; //! Spectrum of Primary
TH1D * fhPrimAccY ; //! Rapidity distribution of primary with accepted daughters
TH1D * fhPrimPhi ; //! Azimutal distribution of primary particles
TH1D * fhPrimAccPhi; //! Azimutal distribution of primary with accepted daughters
-
- ClassDef(AliAnaPi0,7)
+ TH2D * fhPrimOpeningAngle ; //! Opening angle of pair versus pair energy, primaries
+ TH2D * fhPrimCosOpeningAngle ; //! Cosinus of opening angle of pair version pair energy, primaries
+
+ ClassDef(AliAnaPi0,8)
} ;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff