AliAnaPi0: Add possibility to cut on opening angle of pairs, to reject background...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 10 Feb 2010 16:59:59 +0000 (16:59 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 10 Feb 2010 16:59:59 +0000 (16:59 +0000)
AliAnaCalorimeterQA:
 - Histograms to display cells energies for different fractions of the super modules
 - Histograms to calculate the measured time  of a cell in a cluster minus the time of cell with highest energy
 -Histogram to accumulate time measured in cells, fill grid of SuperModule. Only for cell energy larger than 0.5 GeV
- Corrected loop of invariant mass calculation, maximum index in second loop was wrongly set to one unit less.

PWG4/PartCorrDep/AliAnaCalorimeterQA.cxx
PWG4/PartCorrDep/AliAnaCalorimeterQA.h
PWG4/PartCorrDep/AliAnaPi0.cxx
PWG4/PartCorrDep/AliAnaPi0.h

index 5dd2b7c..8cdc649 100755 (executable)
@@ -72,11 +72,13 @@ ClassImp(AliAnaCalorimeterQA)
     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),
@@ -122,13 +124,15 @@ AliAnaCalorimeterQA::AliAnaCalorimeterQA(const AliAnaCalorimeterQA & qa) :
   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),
@@ -224,6 +228,9 @@ AliAnaCalorimeterQA & AliAnaCalorimeterQA::operator = (const AliAnaCalorimeterQA
   fh2Phi = qa.fh2Phi;
   fh2Eta = qa.fh2Eta;
        
+  fhCellTimeSpreadRespectToCellMax = qa.fhCellTimeSpreadRespectToCellMax; 
+  fhCellIdCellLargeTimeSpread      = qa.fhCellIdCellLargeTimeSpread;
+       
   fhNCells    = qa.fhNCells;
   fhAmplitude = qa.fhAmplitude;
   fhAmpId     = qa.fhAmpId;
@@ -236,8 +243,9 @@ AliAnaCalorimeterQA & AliAnaCalorimeterQA::operator = (const AliAnaCalorimeterQA
        
   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  ;
@@ -298,12 +306,12 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
     
        //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();
        }
        
@@ -445,7 +453,7 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        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);
@@ -468,17 +476,26 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        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);
@@ -512,15 +529,18 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        }//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++){
                
@@ -554,10 +574,24 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
                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); 
@@ -1300,6 +1334,7 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
        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",
@@ -1320,7 +1355,8 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
                        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();
@@ -1331,6 +1367,43 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
                                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));
@@ -1359,10 +1432,15 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
                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
@@ -1393,6 +1471,8 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
 
                                //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)
@@ -1436,7 +1516,7 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
                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();
                }
                
@@ -1461,9 +1541,17 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
                                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
@@ -1475,14 +1563,15 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
                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));
@@ -1494,8 +1583,14 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
                                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
@@ -1548,7 +1643,7 @@ void AliAnaCalorimeterQA::ClusterHistograms(const TLorentzVector mom, const Int_
        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){
 
@@ -2166,11 +2261,16 @@ void AliAnaCalorimeterQA::ReadHistograms(TList* outputList)
        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++); 
@@ -2184,26 +2284,37 @@ void AliAnaCalorimeterQA::ReadHistograms(TList* outputList)
        }
        
        //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++);        
 
        }
        
index 43cff95..67a5da0 100755 (executable)
@@ -206,6 +206,9 @@ class AliAnaCalorimeterQA : public AliAnaPartCorrBaseClass {
   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
@@ -227,7 +230,9 @@ class AliAnaCalorimeterQA : public AliAnaPartCorrBaseClass {
   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
 
@@ -337,7 +342,7 @@ class AliAnaCalorimeterQA : public AliAnaPartCorrBaseClass {
   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)
 } ;
 
 
index 38ead0c..d13cc34 100755 (executable)
@@ -26,6 +26,7 @@
 
 // --- ROOT system ---
 #include "TH3.h"
+#include "TH2D.h"
 //#include "Riostream.h"
 #include "TCanvas.h"
 #include "TPad.h"
@@ -45,6 +46,7 @@
 #include "AliESDCaloCluster.h"
 #include "AliESDEvent.h"
 #include "AliAODEvent.h"
+#include "AliNeutralMesonSelection.h"
 
 ClassImp(AliAnaPi0)
 
@@ -52,9 +54,11 @@ 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();
@@ -65,11 +69,13 @@ fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0),
 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
@@ -85,12 +91,14 @@ AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & 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; 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;
   
 }
@@ -129,6 +137,8 @@ void AliAnaPi0::InitParameters()
   fNmaxMixEv = 10;
   fZvtxCut   = 40;
   fCalorimeter  = "PHOS";
+  fUseAngleCut = kFALSE;
+       
 }
 //________________________________________________________________________________________________________________________________________________
 //void AliAnaPi0::Init()
@@ -266,6 +276,20 @@ TList * AliAnaPi0::GetCreateOutputObjects()
                    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
@@ -295,21 +319,35 @@ TList * AliAnaPi0::GetCreateOutputObjects()
     
     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] ;
@@ -354,56 +392,57 @@ void AliAnaPi0::Print(const Option_t * /*opt*/) const
   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;
 }
 
 //____________________________________________________________________________________________________________________________________________________
@@ -452,6 +491,15 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
       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)
@@ -492,6 +540,12 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
        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);                   
@@ -530,80 +584,84 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
   
   //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");
   }    
   
 }      
@@ -611,47 +669,47 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
 //________________________________________________________________________
 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++);
+  
 }
 
 
@@ -661,22 +719,22 @@ void AliAnaPi0::Terminate(TList* outputList)
   //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()));
index 0747a52..8176477 100755 (executable)
@@ -15,6 +15,7 @@
 //Root
 class TList;
 class TH3D ;
+class TH2D ;
 
 //Analysis
 class AliAODEvent ;
@@ -67,20 +68,23 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   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
   
@@ -94,7 +98,10 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   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 
@@ -103,8 +110,10 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   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)
 } ;