PWG4/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 //_________________________________________________________________________
  23
  24
  25 // --- ROOT system ---
  26 #include "TClonesArray.h"
  27 #include "TRefArray.h"
  28 #include "TH3.h"
  29
  30 //---- AliRoot system ----
  31 #include "AliAnaPi0.h"
  32 #include "AliLog.h"
  33 #include "AliCaloPhoton.h"
  34 #include "AliESDEvent.h"
  35 #include "AliAODEvent.h"
  36 #include "AliESDVertex.h"
  37 #include "AliAODVertex.h"
  38 #include "AliESDCaloCluster.h"
  39 #include "AliAODCaloCluster.h"
  40
  41 ClassImp(AliAnaPi0)
  42
  43 //____________________________________________________________________________
  44   AliAnaPi0::AliAnaPi0() : AliAnalysisTaskSE(),
  45   fNCentrBin(1),fNZvertBin(1),fNrpBin(1),fNPID(2),fNmaxMixEv(10),
  46   fCurCentrBin(0),fCurZvertBin(0),fCurRPBin(0),fPtMin(0.),fZvtxCut(40.),
  47   fMinDist(2.),fMinDist2(4.),fMinDist3(5.),fDispCut(1.5),fTOFCut(5.e-9),fPhotPID(0.6),
  48   fEventsList(0x0),fCurrentEvent(0x0),fOutputList(0x0),fhEtalon(0x0),
  49   fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0)
  50 {
  51   //Default Ctor
  52
  53 }
  54 //____________________________________________________________________________
  55   AliAnaPi0::AliAnaPi0(const char *name): AliAnalysisTaskSE(name),
  56   fNCentrBin(1),fNZvertBin(1),fNrpBin(1),fNPID(9),fNmaxMixEv(10),
  57   fCurCentrBin(0),fCurZvertBin(0),fCurRPBin(0),fPtMin(0.),fZvtxCut(40.),
  58   fMinDist(2.),fMinDist2(4.),fMinDist3(5.),fDispCut(1.5),fTOFCut(5.e-9),fPhotPID(0.6),
  59   fEventsList(0x0),fCurrentEvent(0x0),fOutputList(0x0),fhEtalon(0x0),
  60   fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0)
  61 {
  62   //Ctor
  63   DefineOutput(1,TList::Class());
  64 }
  65
  66 //____________________________________________________________________________
  67 AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnalysisTaskSE(ex),
  68   fNCentrBin(1),fNZvertBin(1),fNrpBin(1),fNPID(9),fNmaxMixEv(10),
  69   fCurCentrBin(0),fCurZvertBin(0),fCurRPBin(0),fPtMin(0.),fZvtxCut(40.),
  70   fMinDist(2.),fMinDist2(4.),fMinDist3(5.),fDispCut(1.5),fTOFCut(5.e-9),fPhotPID(0.6),
  71   fEventsList(0x0),fCurrentEvent(0x0),fOutputList(0x0),fhEtalon(0x0),
  72   fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0)
  73 {
  74   // cpy ctor
  75   //Do not need it
  76 }
  77 //_________________________________________________________________________
  78 AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
  79 {
  80   // assignment operator
  81
  82   if(this == &ex)return *this;
  83   ((AliAnalysisTaskSE *)this)->operator=(ex);
  84
  85   return *this;
  86
  87 }
  88 //____________________________________________________________________________
  89 AliAnaPi0::~AliAnaPi0() {
  90   // Remove event containers
  91   if(fEventsList){
  92     for(Int_t ic=0; ic<fNCentrBin; ic++){
  93       for(Int_t iz=0; iz<fNZvertBin; iz++){
  94         for(Int_t irp=0; irp<fNrpBin; irp++){
  95           fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
  96           delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
  97         }
  98       }
  99     }
 100     delete[] fEventsList;
 101     fEventsList=0 ;
 102   }
 103
 104   if(fhEtalon){
 105     delete fhEtalon ;
 106     fhEtalon = 0;
 107   }
 108   if(fhRe1) delete[] fhRe1 ; //Do not delete histograms!
 109   if(fhRe2) delete[] fhRe2 ; //Do not delete histograms!
 110   if(fhRe3) delete[] fhRe3 ; //Do not delete histograms!
 111   if(fhMi1) delete[] fhMi1 ; //Do not delete histograms!
 112   if(fhMi2) delete[] fhMi2 ; //Do not delete histograms!
 113   if(fhMi3) delete[] fhMi3 ; //Do not delete histograms!
 114
 115 }
 116 //________________________________________________________________________
 117 void  AliAnaPi0::UserCreateOutputObjects()
 118 {
 119   // Create histograms to be saved in output file and
 120   // store them in fOutputContainer
 121
 122   AliDebug(1,"Init inv. mass histograms");
 123   OpenFile(1) ;
 124
 125   //create event containers
 126   fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
 127   for(Int_t ic=0; ic<fNCentrBin; ic++){
 128     for(Int_t iz=0; iz<fNZvertBin; iz++){
 129       for(Int_t irp=0; irp<fNrpBin; irp++){
 130         fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
 131       }
 132     }
 133   }
 134
 135   fOutputList = new TList() ;
 136   fOutputList->SetName(GetName()) ;
 137
 138   fhRe1=new TH3D*[fNCentrBin*fNPID] ;
 139   fhRe2=new TH3D*[fNCentrBin*fNPID] ;
 140   fhRe3=new TH3D*[fNCentrBin*fNPID] ;
 141   fhMi1=new TH3D*[fNCentrBin*fNPID] ;
 142   fhMi2=new TH3D*[fNCentrBin*fNPID] ;
 143   fhMi3=new TH3D*[fNCentrBin*fNPID] ;
 144
 145   char key[255] ;
 146   char title[255] ;
 147   for(Int_t ic=0; ic<fNCentrBin; ic++){
 148     for(Int_t ipid=0; ipid<fNPID; ipid++){
 149       //Distance to bad module 1
 150       sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
 151       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 152       fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 153       fhRe1[ic*fNPID+ipid]->SetName(key) ;
 154       fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
 155       fOutputList->Add(fhRe1[ic*fNPID+ipid]) ;
 156
 157       sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
 158       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 159       fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 160       fhMi1[ic*fNPID+ipid]->SetName(key) ;
 161       fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
 162       fOutputList->Add(fhMi1[ic*fNPID+ipid]) ;
 163
 164       //Distance to bad module 2
 165       sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
 166       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 167       fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 168       fhRe2[ic*fNPID+ipid]->SetName(key) ;
 169       fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
 170       fOutputList->Add(fhRe2[ic*fNPID+ipid]) ;
 171
 172       sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
 173       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 174       fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 175       fhMi2[ic*fNPID+ipid]->SetName(key) ;
 176       fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
 177       fOutputList->Add(fhMi2[ic*fNPID+ipid]) ;
 178
 179       //Distance to bad module 3
 180       sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
 181       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 182       fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 183       fhRe3[ic*fNPID+ipid]->SetName(key) ;
 184       fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
 185       fOutputList->Add(fhRe3[ic*fNPID+ipid]) ;
 186
 187       sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
 188       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 189       fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 190       fhMi3[ic*fNPID+ipid]->SetName(key) ;
 191       fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
 192       fOutputList->Add(fhMi3[ic*fNPID+ipid]) ;
 193     }
 194   }
 195
 196   fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
 197                     fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
 198   fOutputList->Add(fhEvents) ;
 199
 200   //Save parameters used for analysis
 201   TString parList ; //this will be list of parameters used for this analysis.
 202   char onePar[255] ;
 203   sprintf(onePar,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
 204   parList+=onePar ;
 205   sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
 206   parList+=onePar ;
 207   sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
 208   parList+=onePar ;
 209   sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
 210   parList+=onePar ;
 211   sprintf(onePar,"Number of different PID used:  %d \n",fNPID) ;
 212   parList+=onePar ;
 213   sprintf(onePar,"Cuts: \n") ;
 214   parList+=onePar ;
 215   sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
 216   parList+=onePar ;
 217   sprintf(onePar,"Minimal P_t: %f \n", fPtMin) ;
 218   parList+=onePar ;
 219   sprintf(onePar,"fMinDist =%f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 220   parList+=onePar ;
 221   sprintf(onePar,"fMinDist2=%f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 222   parList+=onePar ;
 223   sprintf(onePar,"fMinDist3=%f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 224   parList+=onePar ;
 225   sprintf(onePar,"fDispCut =%f (Cut on dispersion, used in PID evaluation) \n",fDispCut) ;
 226   parList+=onePar ;
 227   sprintf(onePar,"fTOFCut  =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
 228   parList+=onePar ;
 229   sprintf(onePar,"fPhotPID =%f (limit for Baesian PID for photon) \n",fPhotPID) ;
 230   parList+=onePar ;
 231
 232   TObjString *oString= new TObjString(parList) ;
 233   fOutputList->Add(oString);
 234
 235 }
 236 //__________________________________________________
 237 void AliAnaPi0::InitParameters()
 238 {
 239   //Initialize the parameters of the analysis.
 240
 241 }
 242  //__________________________________________________
 243 void AliAnaPi0::Init()
 244 {
 245   //Make here all memory allocations
 246   //create etalon histo for all later histograms
 247   if(!fhEtalon){                                                   //  p_T      alpha   d m_gg
 248     fhEtalon = new TH3D("hEtalon","Histo with binning parameters",20,0.,10.,10,0.,1.,200,0.,1.) ;
 249     fhEtalon->SetXTitle("P_{T} (GeV)") ;
 250     fhEtalon->SetYTitle("#alpha") ;
 251     fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
 252   }
 253
 254 }
 255 //__________________________________________________________________
 256 void AliAnaPi0::Print(const Option_t * /*opt*/) const
 257 {
 258   //Print some relevant parameters set for the analysis
 259   printf("Class AliAnaPi0 for gamma-gamma inv.mass construction \n") ;
 260   printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
 261   printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
 262   printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
 263   printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
 264   printf("Number of different PID used:  %d \n",fNPID) ;
 265   printf("Cuts: \n") ;
 266   printf("Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
 267   printf("Minimal P_t: %f \n", fPtMin) ;
 268   printf("fMinDist =%f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 269   printf("fMinDist2=%f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 270   printf("fMinDist3=%f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 271   printf("fDispCut =%f (Cut on dispersion, used in PID evaluation) \n",fDispCut) ;
 272   printf("fTOFCut  =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
 273   printf("fPhotPID =%f (limit for Baesian PID for photon) \n",fPhotPID) ;
 274   printf("------------------------------------------------------\n") ;
 275
 276 }
 277 //__________________________________________________________________
 278 void AliAnaPi0::UserExec(Option_t *)
 279 {
 280   //Process one event and extract photons
 281   //impose cuts on bad modules and bad runs
 282   //fill internal storage for subsequent histo filling
 283
 284   AliVEvent* event = InputEvent();
 285
 286   //Apply some cuts on event: vertex position and centrality range
 287   Int_t iRun=event->GetRunNumber() ;
 288   if(IsBadRun(iRun))
 289     return ;
 290
 291   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
 292   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event) ;
 293   if(esd){
 294      if(!FillFromESD(esd))
 295       return  ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
 296    }
 297    else if(aod){
 298      if(!FillFromAOD(aod))
 299        return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
 300    }
 301    else{
 302      printf("Input neither ESD nor AOD, do nothing \n") ;
 303      return ;
 304    }
 305
 306
 307   fhEvents->Fill(fCurCentrBin+0.5,fCurZvertBin+0.5,fCurRPBin+0.5) ;
 308
 309   if(fCurrentEvent->GetEntriesFast()>0){
 310     //Reduce size for storing
 311     fCurrentEvent->Expand(fCurrentEvent->GetEntriesFast()) ;
 312     FillHistograms() ;
 313   }
 314
 315   PostData(1, fOutputList);
 316
 317 }
 318 //__________________________________________________________________
 319 Bool_t AliAnaPi0::FillFromESD(AliESDEvent * esd){
 320   //Fill photon list from ESD applying
 321   //some cut should be applyed
 322
 323   //Impose cut on vertex and calculate Zvertex bin
 324   const AliESDVertex *esdV = esd->GetVertex() ;
 325   esdV->GetXYZ(fVert);
 326   if(fVert[2]<-fZvtxCut || fVert[2]> fZvtxCut)
 327     return kFALSE ;
 328   fCurZvertBin=(Int_t)(0.5*fNZvertBin*(fVert[2]+fZvtxCut)/fZvtxCut) ;
 329
 330   //Get Centrality and calculate centrality bin
 331   //Does not exist in ESD yet???????
 332   fCurCentrBin=0 ;
 333
 334   //Get Reaction Plain position and calculate RP bin
 335   //does not exist in ESD yet????
 336   fCurRPBin=0 ;
 337
 338   //************************  PHOS *************************************
 339   TRefArray * caloClustersArr  = new TRefArray();
 340   esd->GetPHOSClusters(caloClustersArr);
 341
 342   const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
 343   //if fCurrentEvent!=0 it was not filled in previous event, still empty and can be used
 344   if(!fCurrentEvent)
 345     fCurrentEvent = new TClonesArray("AliCaloPhoton",kNumberOfPhosClusters) ;
 346   Int_t inList=0 ;
 347
 348   // loop over the PHOS Cluster
 349   for(Int_t i = 0 ; i < kNumberOfPhosClusters ; i++) {
 350     AliESDCaloCluster * calo = (AliESDCaloCluster *) caloClustersArr->At(i) ;
 351
 352     //Make some tests
 353     if(calo->E()<fPtMin)
 354       continue ;
 355
 356     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad in cm
 357     if(distBad<0.)distBad=9999. ; //workout strange convension dist = -1. ;
 358     if(distBad<fMinDist) //In bad channel (cristall size 2.2x2.2 cm)
 359       continue ;
 360
 361     new((*fCurrentEvent)[inList])AliCaloPhoton() ;
 362     AliCaloPhoton * ph = static_cast<AliCaloPhoton*>(fCurrentEvent->At(inList)) ;
 363     inList++ ;
 364
 365     //Set Momentum
 366     calo->GetMomentum(*ph,fVert);
 367
 368     //Dispersion/lambdas
 369     Double_t disp=calo->GetClusterDisp()  ;
 370 //    Double_t m20=calo->GetM20() ;
 371 //    Double_t m02=calo->GetM02() ;
 372     ph->SetDispBit(disp<fDispCut) ;
 373
 374     //TOF
 375     Double_t tof=calo->GetTOF()  ;
 376     ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
 377
 378     //Charged veto
 379 //    Double_t cpvR=calo->GetEmcCpvDistance() ;
 380     Int_t itr=calo->GetNTracksMatched();  //number of track matched
 381     ph->SetCPVBit(itr>0) ;  //Temporary cut, should we evaluate distance?
 382
 383     //Overall PID
 384     Double_t *pid=calo->GetPid();
 385     ph->SetPCAPID(pid[AliPID::kPhoton]>fPhotPID) ;
 386
 387     //Set Distance to Bad channel
 388     if(distBad>fMinDist3)
 389       ph->SetDistToBad(2) ;
 390     else if(distBad>fMinDist2)
 391         ph->SetDistToBad(1) ;
 392     else
 393        ph->SetDistToBad(0) ;
 394
 395   }//loop
 396   return kTRUE ;
 397 }
 398 //__________________________________________________________________
 399 Bool_t AliAnaPi0::FillFromAOD(AliAODEvent * aod){
 400   //Fill photon list from AOD applying
 401   //some cuts
 402
 403   //Impose cut on vertex and calculate Zvertex bin
 404   const AliAODVertex *aodV = aod->GetPrimaryVertex() ;
 405   fVert[0]=aodV->GetX();
 406   fVert[1]=aodV->GetY();
 407   fVert[2]=aodV->GetZ();
 408   if(fVert[2]<-fZvtxCut || fVert[2]> fZvtxCut)
 409     return kFALSE ;
 410   fCurZvertBin=(Int_t)(0.5*fNZvertBin*(fVert[2]+fZvtxCut)/fZvtxCut) ;
 411
 412   //Get Centrality and calculate centrality bin
 413   //Does not exist in ESD yet???????
 414   fCurCentrBin=0 ;
 415
 416   //Get Reaction Plain position and calculate RP bin
 417   //does not exist in ESD yet????
 418   fCurRPBin=0 ;
 419
 420   //************************  PHOS *************************************
 421   const Int_t kNumberOfPhosClusters   = aod->GetNCaloClusters() ;
 422   if(!fCurrentEvent)
 423     fCurrentEvent = new TClonesArray("AliCaloPhoton",kNumberOfPhosClusters) ;
 424   Int_t inList=0 ;
 425
 426   // loop over the PHOS Cluster
 427   for(Int_t i = 0 ; i < kNumberOfPhosClusters ; i++) {
 428     AliAODCaloCluster * calo = aod->GetCaloCluster(i);
 429
 430     //Make some tests
 431     if(calo->E()<fPtMin)
 432       continue ;
 433
 434     Double_t distBad=calo->GetDistToBadChannel() ;
 435     if(distBad<fMinDist) //In bad channel
 436       continue ;
 437
 438     new((*fCurrentEvent)[inList])AliCaloPhoton() ;
 439     AliCaloPhoton * ph = static_cast<AliCaloPhoton*>(fCurrentEvent->At(inList)) ;
 440     inList++ ;
 441
 442     //Set Momentum
 443     calo->GetMomentum(*ph,fVert);
 444
 445     //Dispersion/lambdas
 446     Double_t disp=calo->GetDispersion()  ;
 447 //    Double_t m20=calo->GetM20() ;
 448 //    Double_t m02=calo->GetM02() ;
 449     ph->SetDispBit(disp<fDispCut) ;
 450
 451     //TOF
 452     Double_t tof=calo->GetTOF()  ;
 453     ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
 454
 455     //Charged veto
 456 //    Double_t cpvR=calo->GetEmcCpvDistance() ;
 457     Int_t itr=calo->GetNTracksMatched();  //number of track matched
 458     ph->SetCPVBit(itr>0) ;  //Temporary cut !!!!
 459
 460     //Overall PID
 461     Double_t pid[13];
 462     calo->GetPID(pid);
 463     ph->SetPCAPID(pid[AliAODCluster::kPhoton]>fPhotPID) ;
 464
 465     //Distance to Bad
 466     if(distBad>fMinDist3)
 467       ph->SetDistToBad(2) ;
 468     else if(distBad>fMinDist2)
 469         ph->SetDistToBad(1) ;
 470     else
 471        ph->SetDistToBad(0) ;
 472
 473
 474   }//loop
 475   return kTRUE ;
 476 }
 477 //__________________________________________________________________
 478 void  AliAnaPi0::FillHistograms()
 479 {
 480   //Fill Real and Mixed invariant mass histograms
 481   //then add current event to the buffer and
 482   //remove redundant events from buffer if necessary
 483
 484   //Fill histograms only if list of particles for current event was created
 485   if(fCurrentEvent==0)
 486     return ;
 487
 488   //Fill Real distribution
 489   Int_t nPhot = fCurrentEvent->GetEntriesFast() ;
 490   for(Int_t i1=0; i1<nPhot-1; i1++){
 491     AliCaloPhoton * p1 = static_cast<AliCaloPhoton*>(fCurrentEvent->At(i1)) ;
 492     for(Int_t i2=i1+1; i2<nPhot; i2++){
 493       AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(fCurrentEvent->At(i2)) ;
 494       Double_t m =(*p1 + *p2).M() ;
 495       Double_t pt=(*p1 + *p2).Pt();
 496       Double_t a = TMath::Abs(p1->Energy()-p2->Energy())/(p1->Energy()+p2->Energy()) ;
 497       for(Int_t ipid=0; ipid<fNPID; ipid++){
 498         if(p1->IsPIDOK(ipid)&&p2->IsPIDOK(ipid)){
 499           fhRe1[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 500           if(p1->DistToBad()>0 && p2->DistToBad()>0){
 501             fhRe2[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 502             if(p1->DistToBad()>1 && p2->DistToBad()>1){
 503               fhRe3[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 504             }
 505           }
 506         }
 507       }
 508     }
 509   }
 510
 511   //Fill mixed
 512   TList * evMixList=fEventsList[fCurCentrBin*fNZvertBin*fNrpBin+fCurZvertBin*fNrpBin+fCurRPBin] ;
 513   Int_t nMixed = evMixList->GetSize() ;
 514   for(Int_t ii=0; ii<nMixed; ii++){
 515     TClonesArray* ev2=dynamic_cast<TClonesArray*>(evMixList->At(ii));
 516     Int_t nPhot2=ev2->GetEntriesFast() ;
 517     for(Int_t i1=0; i1<nPhot; i1++){
 518       AliCaloPhoton * p1 = static_cast<AliCaloPhoton*>(fCurrentEvent->At(i1)) ;
 519       for(Int_t i2=0; i2<nPhot2; i2++){
 520         AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(ev2->At(i2)) ;
 521         Double_t m =(*p1 + *p2).M() ;
 522         Double_t pt=(*p1 + *p2).Pt();
 523         Double_t a = TMath::Abs(p1->Energy()-p2->Energy())/(p1->Energy()+p2->Energy()) ;
 524         for(Int_t ipid=0; ipid<fNPID; ipid++){
 525           if(p1->IsPIDOK(ipid)&&p2->IsPIDOK(ipid)){
 526             fhMi1[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 527             if(p1->DistToBad()>0 && p2->DistToBad()>0){
 528               fhMi2[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 529               if(p1->DistToBad()>1 && p2->DistToBad()>1){
 530                 fhMi3[fCurCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 531               }
 532             }
 533           }
 534         }
 535       }
 536     }
 537   }
 538
 539   //Add current event to buffer and Remove redandant events
 540   if(fCurrentEvent->GetEntriesFast()>0){
 541     evMixList->AddFirst(fCurrentEvent) ;
 542     fCurrentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
 543     if(evMixList->GetSize()>=fNmaxMixEv){
 544       TClonesArray * tmp = dynamic_cast<TClonesArray*>(evMixList->Last()) ;
 545       evMixList->RemoveLast() ;
 546       delete tmp ;
 547     }
 548   }
 549   else{ //empty event
 550    delete fCurrentEvent ;
 551    fCurrentEvent=0 ;
 552   }
 553 }
 554 //______________________________________________________________________________
 555 void AliAnaPi0::Terminate(Option_t *)
 556 {
 557 }