PWG4/PartCorrDep/AliAnaPi0.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 // Class to collect two-photon invariant mass distributions for
  19 // extractin raw pi0 yield.
  20 //
  21 //-- Author: Dmitri Peressounko (RRC "KI")
  22 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
  23 //-- and Gustavo Conesa (INFN-Frascati)
  24 //_________________________________________________________________________
  25
  26
  27 // --- ROOT system ---
  28 #include "TH3.h"
  29 //#include "Riostream.h"
  30 #include "TCanvas.h"
  31 #include "TPad.h"
  32 #include "TROOT.h"
  33
  34 //---- AliRoot system ----
  35 #include "AliAnaPi0.h"
  36 #include "AliCaloTrackReader.h"
  37 #include "AliCaloPID.h"
  38 #include "AliStack.h"
  39
  40 #ifdef __PHOSGEO__
  41         #include "AliPHOSGeoUtils.h"
  42 #endif
  43
  44 ClassImp(AliAnaPi0)
  45
  46 //________________________________________________________________________________________________________________________________________________
  47 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
  48 fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
  49 fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
  50 fEventsList(0x0), fhEtalon(0x0),
  51 fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
  52 fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0)
  53 {
  54 //Default Ctor
  55  InitParameters();
  56
  57 }
  58
  59 //________________________________________________________________________________________________________________________________________________
  60 AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),
  61 fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
  62 fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter),
  63 fEventsList(ex.fEventsList), fhEtalon(ex.fhEtalon),
  64 fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
  65 fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY),
  66 fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi)
  67 {
  68         // cpy ctor
  69         //Do not need it
  70 }
  71
  72 //________________________________________________________________________________________________________________________________________________
  73 AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
  74 {
  75         // assignment operator
  76
  77         if(this == &ex)return *this;
  78         ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
  79
  80         fNCentrBin = ex.fNCentrBin  ; fNZvertBin = ex.fNZvertBin  ; fNrpBin = ex.fNrpBin  ;
  81         fNPID = ex.fNPID  ; fNmaxMixEv = ex.fNmaxMixEv  ; fZvtxCut = ex.fZvtxCut  ;  fCalorimeter = ex.fCalorimeter  ;
  82         fEventsList = ex.fEventsList  ;  fhEtalon = ex.fhEtalon  ;
  83         fhRe1 = ex.fhRe1  ; fhMi1 = ex.fhMi1  ; fhRe2 = ex.fhRe2  ; fhMi2 = ex.fhMi2  ;
  84         fhRe3 = ex.fhRe3  ; fhMi3 = ex.fhMi3  ; fhEvents = ex.fhEvents  ;
  85         fhPrimPt = ex.fhPrimPt  ;  fhPrimAccPt = ex.fhPrimAccPt  ;  fhPrimY = ex.fhPrimY  ;
  86         fhPrimAccY = ex.fhPrimAccY  ;  fhPrimPhi = ex.fhPrimPhi  ;  fhPrimAccPhi = ex.fhPrimAccPhi ;
  87
  88         return *this;
  89
  90 }
  91
  92 //________________________________________________________________________________________________________________________________________________
  93 AliAnaPi0::~AliAnaPi0() {
  94         // Remove event containers
  95         if(fEventsList){
  96                 for(Int_t ic=0; ic<fNCentrBin; ic++){
  97                         for(Int_t iz=0; iz<fNZvertBin; iz++){
  98                                 for(Int_t irp=0; irp<fNrpBin; irp++){
  99                                         fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
 100                                         delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
 101                                 }
 102                         }
 103                 }
 104                 delete[] fEventsList;
 105                 fEventsList=0 ;
 106         }
 107
 108 #ifdef __PHOSGEO__
 109     if(fPHOSGeo) delete fPHOSGeo ;
 110 #endif
 111 }
 112
 113 //________________________________________________________________________________________________________________________________________________
 114 void AliAnaPi0::InitParameters()
 115 {
 116 //Init parameters when first called the analysis
 117 //Set default parameters
 118         SetInputAODName("photons");
 119         fNCentrBin = 1;
 120         fNZvertBin = 1;
 121         fNrpBin    = 1;
 122     fNPID      = 9;
 123         fNmaxMixEv = 10;
 124         fZvtxCut   = 40;
 125         fCalorimeter  = "PHOS";
 126 }
 127 //________________________________________________________________________________________________________________________________________________
 128 void AliAnaPi0::Init()
 129 {
 130 //Init some data members needed in analysis
 131
 132 //Histograms binning and range
 133         if(!fhEtalon){                                                   //  p_T      alpha   d m_gg
 134                 fhEtalon = new TH3D("hEtalon","Histo with binning parameters",50,0.,25.,10,0.,1.,200,0.,1.) ;
 135                 fhEtalon->SetXTitle("P_{T} (GeV)") ;
 136                 fhEtalon->SetYTitle("#alpha") ;
 137                 fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
 138         }
 139
 140
 141 //create event containers
 142         fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
 143
 144         for(Int_t ic=0; ic<fNCentrBin; ic++){
 145                 for(Int_t iz=0; iz<fNZvertBin; iz++){
 146                         for(Int_t irp=0; irp<fNrpBin; irp++){
 147                                 fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
 148                         }
 149                 }
 150         }
 151
 152                 //If Geometry library loaded, do geometry selection during analysis.
 153 #ifdef __PHOSGEO__
 154         printf("PHOS geometry initialized!\n");
 155         fPHOSGeo = new AliPHOSGeoUtils("PHOSgeo") ;
 156 #endif
 157
 158 }
 159
 160 //________________________________________________________________________________________________________________________________________________
 161 TList * AliAnaPi0::GetCreateOutputObjects()
 162 {
 163         // Create histograms to be saved in output file and
 164         // store them in fOutputContainer
 165
 166
 167         TList * outputContainer = new TList() ;
 168         outputContainer->SetName(GetName());
 169
 170         fhRe1=new TH3D*[fNCentrBin*fNPID] ;
 171         fhRe2=new TH3D*[fNCentrBin*fNPID] ;
 172         fhRe3=new TH3D*[fNCentrBin*fNPID] ;
 173         fhMi1=new TH3D*[fNCentrBin*fNPID] ;
 174         fhMi2=new TH3D*[fNCentrBin*fNPID] ;
 175         fhMi3=new TH3D*[fNCentrBin*fNPID] ;
 176
 177         char key[255] ;
 178         char title[255] ;
 179
 180         for(Int_t ic=0; ic<fNCentrBin; ic++){
 181                 for(Int_t ipid=0; ipid<fNPID; ipid++){
 182
 183                         //Distance to bad module 1
 184                         sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
 185                         sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 186
 187                         fhEtalon->Clone(key);
 188                         fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 189                         fhRe1[ic*fNPID+ipid]->SetName(key) ;
 190                         fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
 191                         outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
 192
 193                         sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
 194                         sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 195                         fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 196                         fhMi1[ic*fNPID+ipid]->SetName(key) ;
 197                         fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
 198                         outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
 199
 200                         //Distance to bad module 2
 201                         sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
 202                         sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 203                         fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 204                         fhRe2[ic*fNPID+ipid]->SetName(key) ;
 205                         fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
 206                         outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
 207
 208                         sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
 209                         sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 210                         fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 211                         fhMi2[ic*fNPID+ipid]->SetName(key) ;
 212                         fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
 213                         outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
 214
 215                         //Distance to bad module 3
 216                         sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
 217                         sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 218                         fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 219                         fhRe3[ic*fNPID+ipid]->SetName(key) ;
 220                         fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
 221                         outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
 222
 223                         sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
 224                         sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 225                         fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 226                         fhMi3[ic*fNPID+ipid]->SetName(key) ;
 227                         fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
 228                         outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
 229                 }
 230     }
 231
 232
 233         fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
 234                                           fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
 235         outputContainer->Add(fhEvents) ;
 236
 237         //Histograms filled only if MC data is requested
 238         if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
 239                 if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
 240                         fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
 241                         fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",fhEtalon->GetXaxis()->GetNbins(),
 242                                                                    fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
 243                 }
 244                 else{
 245                         fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
 246                         fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",
 247                                                                    fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
 248                 }
 249                 outputContainer->Add(fhPrimPt) ;
 250                 outputContainer->Add(fhPrimAccPt) ;
 251
 252                 fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",100,-5.,5.) ;
 253                 outputContainer->Add(fhPrimY) ;
 254
 255                 fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",100,-5.,5.) ;
 256                 outputContainer->Add(fhPrimAccY) ;
 257
 258                 fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",180,0.,360.) ;
 259                 outputContainer->Add(fhPrimPhi) ;
 260
 261                 fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",180,-0.,360.) ;
 262                 outputContainer->Add(fhPrimAccPhi) ;
 263         }
 264
 265         //Save parameters used for analysis
 266         TString parList ; //this will be list of parameters used for this analysis.
 267         char onePar[255] ;
 268         sprintf(onePar,"--- AliAnaPi0 ---\n") ;
 269         parList+=onePar ;
 270         sprintf(onePar,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
 271         parList+=onePar ;
 272         sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
 273         parList+=onePar ;
 274         sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
 275         parList+=onePar ;
 276         sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
 277         parList+=onePar ;
 278         sprintf(onePar,"Number of different PID used:  %d \n",fNPID) ;
 279         parList+=onePar ;
 280         sprintf(onePar,"Cuts: \n") ;
 281         parList+=onePar ;
 282         sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
 283         parList+=onePar ;
 284         sprintf(onePar,"Calorimeter: %s \n",fCalorimeter.Data()) ;
 285         parList+=onePar ;
 286
 287         TObjString *oString= new TObjString(parList) ;
 288         outputContainer->Add(oString);
 289
 290         return outputContainer;
 291 }
 292
 293 //_________________________________________________________________________________________________________________________________________________
 294 void AliAnaPi0::Print(const Option_t * /*opt*/) const
 295 {
 296         //Print some relevant parameters set for the analysis
 297         AliAnaPartCorrBaseClass::Print(" ");
 298         printf("Class AliAnaPi0 for gamma-gamma inv.mass construction \n") ;
 299         printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
 300         printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
 301         printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
 302         printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
 303         printf("Number of different PID used:  %d \n",fNPID) ;
 304         printf("Cuts: \n") ;
 305         printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
 306         printf("------------------------------------------------------\n") ;
 307 }
 308
 309
 310 //____________________________________________________________________________________________________________________________________________________
 311 void AliAnaPi0::MakeAnalysisFillHistograms()
 312 {
 313         //Process one event and extract photons from AOD branch
 314         // filled with AliAnaPhoton and fill histos with invariant mass
 315
 316         //Apply some cuts on event: vertex position and centrality range
 317         Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
 318         if(IsBadRun(iRun)) return ;
 319
 320         Double_t vert[]={0,0,0} ; //vertex ;
 321         GetReader()->GetVertex(vert);
 322         if(vert[2]<-fZvtxCut || vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
 323
 324         //Get Centrality and calculate centrality bin
 325         //Does not exist in ESD yet???????
 326         Int_t curCentrBin=0 ;
 327
 328         //Get Reaction Plain position and calculate RP bin
 329         //does not exist in ESD yet????
 330         Int_t curRPBin=0 ;
 331
 332         Int_t curZvertBin=(Int_t)(0.5*fNZvertBin*(vert[2]+fZvtxCut)/fZvtxCut) ;
 333
 334         fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
 335
 336         Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
 337         if(GetDebug() > 1) printf("AliAnaPi0::FillHistos: photon entries %d\n", nPhot);
 338
 339         for(Int_t i1=0; i1<nPhot-1; i1++){
 340                 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
 341                 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
 342                 for(Int_t i2=i1+1; i2<nPhot; i2++){
 343                         AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
 344                         TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
 345                         Double_t m  = (photon1 + photon2).M() ;
 346                         Double_t pt = (photon1 + photon2).Pt();
 347                         Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
 348                         if(GetDebug() > 2)
 349                                 printf("AliAnaPi0::FillHistos: Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
 350                                                 p1->Pt(), p2->Pt(), pt,m,a);
 351                         for(Int_t ipid=0; ipid<fNPID; ipid++)
 352                         {
 353                                 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
 354                                         fhRe1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 355                                         if(p1->DistToBad()>0 && p2->DistToBad()>0){
 356                                                 fhRe2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 357                                                 if(p1->DistToBad()>1 && p2->DistToBad()>1){
 358                                                         fhRe3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 359                                                 }
 360                                         }
 361                                 }
 362                         }
 363                 }
 364         }
 365
 366         //Fill mixed
 367
 368         TList * evMixList=fEventsList[curCentrBin*fNZvertBin*fNrpBin+curZvertBin*fNrpBin+curRPBin] ;
 369         Int_t nMixed = evMixList->GetSize() ;
 370         for(Int_t ii=0; ii<nMixed; ii++){
 371                 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
 372                 Int_t nPhot2=ev2->GetEntriesFast() ;
 373                 Double_t m = -999;
 374                 if(GetDebug() > 1) printf("AliAnaPi0::FillHistos: Mixed event %d photon entries %d\n", ii, nPhot);
 375
 376                 for(Int_t i1=0; i1<nPhot; i1++){
 377                         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
 378                         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
 379                         for(Int_t i2=0; i2<nPhot2; i2++){
 380                                 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
 381
 382                                 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
 383                                 m =           (photon1+photon2).M() ;
 384                                 Double_t pt = (photon1 + photon2).Pt();
 385                                 Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
 386                                 if(GetDebug() > 2)
 387                                         printf("AliAnaPi0::FillHistos: Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
 388                                                         p1->Pt(), p2->Pt(), pt,m,a);
 389                                 for(Int_t ipid=0; ipid<fNPID; ipid++){
 390                                 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
 391                                                 fhMi1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 392                                                 if(p1->DistToBad()>0 && p2->DistToBad()>0){
 393                                                         fhMi2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 394                                                         if(p1->DistToBad()>1 && p2->DistToBad()>1){
 395                                                                 fhMi3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 396                                                         }
 397
 398                                                 }
 399                                         }
 400                                 }
 401                         }
 402                 }
 403         }
 404
 405         TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
 406         //Add current event to buffer and Remove redandant events
 407         if(currentEvent->GetEntriesFast()>0){
 408                 evMixList->AddFirst(currentEvent) ;
 409                 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
 410                 if(evMixList->GetSize()>=fNmaxMixEv)
 411                 {
 412                         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
 413                         evMixList->RemoveLast() ;
 414                         delete tmp ;
 415                 }
 416         }
 417         else{ //empty event
 418                 delete currentEvent ;
 419                 currentEvent=0 ;
 420         }
 421
 422         //Acceptance
 423         AliStack * stack = GetMCStack();
 424         if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
 425                 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
 426                         TParticle * prim = stack->Particle(i) ;
 427                         if(prim->GetPdgCode() == 111){
 428                                 Double_t pi0Pt = prim->Pt() ;
 429                                 //printf("pi0, pt %2.2f\n",pi0Pt);
 430                                 Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
 431                                 Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
 432                                 if(TMath::Abs(pi0Y) < 0.5){
 433                                         fhPrimPt->Fill(pi0Pt) ;
 434                                 }
 435                                 fhPrimY  ->Fill(pi0Y) ;
 436                                 fhPrimPhi->Fill(phi) ;
 437
 438                                 //Check if both photons hit Calorimeter
 439                                 Int_t iphot1=prim->GetFirstDaughter() ;
 440                                 Int_t iphot2=prim->GetLastDaughter() ;
 441                                 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
 442                                         TParticle * phot1 = stack->Particle(iphot1) ;
 443                                         TParticle * phot2 = stack->Particle(iphot2) ;
 444                                         if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
 445                                                 //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",
 446                                                 //      phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
 447                                                 Bool_t inacceptance = kFALSE;
 448 #ifdef __PHOSGEO__
 449                                                 Int_t mod ;
 450                                                 Double_t x,z ;
 451                                                 if(fCalorimeter == "PHOS" && fPHOSGeo->ImpactOnEmc(phot1,mod,z,x) && fPHOSGeo->ImpactOnEmc(phot1,mod,z,x))
 452                                                         inacceptance = kTRUE;
 453                                                 //printf("In REAL PHOS acceptance? %d\n",inacceptance);
 454 #else
 455                                                 TLorentzVector lv1, lv2;
 456                                                 phot1->Momentum(lv1);
 457                                                 phot2->Momentum(lv2);
 458                                                 if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter))
 459                                                         inacceptance = kTRUE ;
 460                                                 //printf("In %s fidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 461 #endif
 462                                                 if(inacceptance){
 463                                                         fhPrimAccPt->Fill(pi0Pt) ;
 464                                                         fhPrimAccPhi->Fill(phi) ;
 465                                                         fhPrimAccY->Fill(pi0Y) ;
 466                                                 }//Accepted
 467                                         }// 2 photons
 468                                 }//Check daughters exist
 469                         }// Primary pi0
 470                 }//loop on primaries
 471         }//stack exists and data is MC
 472
 473 }
 474
 475
 476
 477
 478 //____________________________________________________________________________________________________________________________________________________
 479 void AliAnaPi0::Terminate()
 480 {
 481   //Do some calculations and plots from the final histograms.
 482
 483   printf(" *** %s Terminate:", GetName()) ;
 484   printf("        Mgg Real        : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(),   fhRe1[0]->GetRMS() ) ;
 485
 486   TCanvas  * cIM = new TCanvas("cIM", "", 400, 10, 600, 700) ;
 487   cIM->Divide(2, 2);
 488
 489   cIM->cd(1) ;
 490   //gPad->SetLogy();
 491   TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ();
 492   hIMAllPt->SetLineColor(2);
 493   hIMAllPt->SetTitle("No cut on  p_{T, #gamma#gamma} ");
 494   hIMAllPt->Draw();
 495
 496   cIM->cd(2) ;
 497   TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone("IMPt5");
 498   hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
 499   TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D("z");
 500   hIMPt5->SetLineColor(2);
 501   hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
 502   hIMPt5->Draw();
 503
 504   cIM->cd(3) ;
 505   TH3F * hRe1Pt10 =  (TH3F*)fhRe1[0]->Clone("IMPt10");
 506   hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
 507   TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D("z");
 508   hIMPt10->SetLineColor(2);
 509   hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
 510   hIMPt10->Draw();
 511
 512   cIM->cd(4) ;
 513   TH3F * hRe1Pt20 =  (TH3F*)fhRe1[0]->Clone("IMPt20");
 514   hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
 515   TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D("z");
 516   hIMPt20->SetLineColor(2);
 517   hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
 518   hIMPt20->Draw();
 519
 520   cIM->Print("Mgg.eps");
 521
 522
 523  TCanvas  * cPt = new TCanvas("cPt", "", 400, 10, 600, 700) ;
 524   cPt->Divide(2, 2);
 525
 526   cPt->cd(1) ;
 527   //gPad->SetLogy();
 528   TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
 529   hPt->SetLineColor(2);
 530   hPt->SetTitle("No cut on  M_{#gamma#gamma} ");
 531   hPt->Draw();
 532
 533   cPt->cd(2) ;
 534   TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone("PtIM1");
 535   hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
 536   TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
 537   hPtIM1->SetLineColor(2);
 538   hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
 539   hPtIM1->Draw();
 540
 541   cPt->cd(3) ;
 542   TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone("PtIM2");
 543   hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
 544   TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
 545   hPtIM2->SetLineColor(2);
 546   hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
 547   hPtIM2->Draw();
 548
 549   cPt->cd(4) ;
 550   TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone("PtIM3");
 551   hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
 552   TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
 553   hPtIM3->SetLineColor(2);
 554   hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
 555   hPtIM3->Draw();
 556
 557   cPt->Print("Pt.eps");
 558
 559
 560   char line[1024] ;
 561   sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
 562   gROOT->ProcessLine(line);
 563   sprintf(line, ".!rm -fR *.eps");
 564   gROOT->ProcessLine(line);
 565
 566   printf("!! All the eps files are in %s.tar.gz !!!", GetName());
 567
 568 }
 569
 570
 571
 572
 573