PWGGA/CaloTrackCorrelations/AliAnaClusterPileUp.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 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
  16 //_________________________________________________________________________
  17 //
  18 // Class for the study of Pile-up effect on
  19 // Calorimeter clusters.
  20 // Open time cuts in reader.
  21 //
  22 //-- Author: Gustavo Conesa (CNRS-LPSC-Grenoble)
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25 // --- ROOT system ---
  26 #include <TH2F.h>
  27 #include <TClonesArray.h>
  28 #include <TObjString.h>
  29
  30 // --- Analysis system ---
  31 #include "AliAnaClusterPileUp.h"
  32 #include "AliCaloTrackReader.h"
  33 #include "AliFiducialCut.h"
  34 #include "AliVCluster.h"
  35 #include "AliAODEvent.h"
  36 #include "AliESDEvent.h"
  37
  38 ClassImp(AliAnaClusterPileUp)
  39
  40 //___________________________________________
  41 AliAnaClusterPileUp::AliAnaClusterPileUp() :
  42 AliAnaCaloTrackCorrBaseClass(),
  43 fCalorimeter(""),                     fNCellsCut(0),
  44 // Histograms
  45 fhTimePtNoCut(0),                     fhTimePtSPD(0),
  46 fhTimeNPileUpVertSPD(0),              fhTimeNPileUpVertTrack(0),
  47 fhTimeNPileUpVertContributors(0),
  48 fhTimePileUpMainVertexZDistance(0),   fhTimePileUpMainVertexZDiamond(0),
  49 fhClusterMultSPDPileUp(),             fhClusterMultNoPileUp(),
  50 fhEtaPhiBC0(0),  fhEtaPhiBCPlus(0),   fhEtaPhiBCMinus(0),
  51 fhEtaPhiBC0PileUpSPD(0),
  52 fhEtaPhiBCPlusPileUpSPD(0),           fhEtaPhiBCMinusPileUpSPD(0),
  53 fhPtNPileUpSPDVtx(0),                 fhPtNPileUpTrkVtx(0),
  54 fhPtNPileUpSPDVtxTimeCut(0),          fhPtNPileUpTrkVtxTimeCut(0),
  55 fhPtNPileUpSPDVtxTimeCut2(0),         fhPtNPileUpTrkVtxTimeCut2(0)
  56 {
  57   //default ctor
  58
  59   for(Int_t i = 0; i < 7; i++)
  60   {
  61     fhPtPileUp       [i] = 0;
  62     fhPtNeutralPileUp[i] = 0;
  63
  64     fhLambda0PileUp       [i] = 0;
  65     fhLambda0NeutralPileUp[i] = 0;
  66
  67     fhClusterEFracLongTimePileUp  [i] = 0;
  68
  69     fhClusterCellTimePileUp       [i] = 0;
  70     fhClusterTimeDiffPileUp       [i] = 0;
  71     fhClusterTimeDiffNeutralPileUp[i] = 0;
  72
  73   }
  74
  75   for(Int_t i = 0; i < 4; i++)
  76   {
  77     fhClusterMultSPDPileUp[i] = 0;
  78     fhClusterMultNoPileUp [i] = 0;
  79   }
  80
  81   //Initialize parameters
  82   InitParameters();
  83
  84 }
  85
  86 //___________________________________________
  87 TObjString *  AliAnaClusterPileUp::GetAnalysisCuts()
  88 {
  89   //Save parameters used for analysis
  90   TString parList ; //this will be list of parameters used for this analysis.
  91   const Int_t buffersize = 255;
  92   char onePar[buffersize] ;
  93
  94   snprintf(onePar,buffersize,"--- AliAnaClusterPileUp ---\n") ;
  95   parList+=onePar ;
  96   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
  97   parList+=onePar ;
  98
  99   //Get parameters set in base class.
 100   //parList += GetBaseParametersList() ;
 101
 102   return new TObjString(parList) ;
 103 }
 104
 105 //________________________________________________________________________
 106 TList *  AliAnaClusterPileUp::GetCreateOutputObjects()
 107 {
 108   // Create histograms to be saved in output file and
 109   // store them in outputContainer
 110   TList * outputContainer = new TList() ;
 111   outputContainer->SetName("PhotonHistos") ;
 112
 113   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();  Float_t ptimecluster  = GetHistogramRanges()->GetHistoPtMax();  Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 114   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
 115   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 116   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax   = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin   = GetHistogramRanges()->GetHistoShowerShapeMin();
 117   Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
 118
 119
 120   fhTimePtNoCut  = new TH2F ("hTimePt_NoCut","time of cluster vs pT of clusters, no event selection", nptbins,ptmin,ptimecluster, ntimebins,timemin,timemax);
 121   fhTimePtNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 122   fhTimePtNoCut->SetYTitle("#it{time} (ns)");
 123   outputContainer->Add(fhTimePtNoCut);
 124
 125   fhTimePtSPD  = new TH2F ("hTimePt_SPD","time of cluster vs pT of clusters, SPD Pile-up events", nptbins,ptmin,ptimecluster, ntimebins,timemin,timemax);
 126   fhTimePtSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 127   fhTimePtSPD->SetYTitle("#it{time} (ns)");
 128   outputContainer->Add(fhTimePtSPD);
 129
 130   fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,20,0,20);
 131   fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
 132   fhTimeNPileUpVertSPD->SetXTitle("#it{time} (ns)");
 133   outputContainer->Add(fhTimeNPileUpVertSPD);
 134
 135   fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 20,0,20 );
 136   fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
 137   fhTimeNPileUpVertTrack->SetXTitle("#it{time} (ns)");
 138   outputContainer->Add(fhTimeNPileUpVertTrack);
 139
 140   fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
 141   fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
 142   fhTimeNPileUpVertContributors->SetXTitle("#it{time} (ns)");
 143   outputContainer->Add(fhTimeNPileUpVertContributors);
 144
 145   fhTimePileUpMainVertexZDistance  = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50);
 146   fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
 147   fhTimePileUpMainVertexZDistance->SetXTitle("#it{time} (ns)");
 148   outputContainer->Add(fhTimePileUpMainVertexZDistance);
 149
 150   fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
 151   fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
 152   fhTimePileUpMainVertexZDiamond->SetXTitle("#it{time} (ns)");
 153   outputContainer->Add(fhTimePileUpMainVertexZDiamond);
 154
 155   fhEtaPhiBC0  = new TH2F ("hEtaPhiBC0","eta-phi for clusters tof corresponding to BC=0",netabins,etamin,etamax, nphibins,phimin,phimax);
 156   fhEtaPhiBC0->SetXTitle("#eta ");
 157   fhEtaPhiBC0->SetYTitle("#phi (rad)");
 158   outputContainer->Add(fhEtaPhiBC0);
 159
 160   fhEtaPhiBCPlus  = new TH2F ("hEtaPhiBCPlus","eta-phi for clusters tof corresponding to BC>0",netabins,etamin,etamax, nphibins,phimin,phimax);
 161   fhEtaPhiBCPlus->SetXTitle("#eta ");
 162   fhEtaPhiBCPlus->SetYTitle("#phi (rad)");
 163   outputContainer->Add(fhEtaPhiBCPlus);
 164
 165   fhEtaPhiBCMinus  = new TH2F ("hEtaPhiBCMinus","eta-phi for clusters tof corresponding to BC<0",netabins,etamin,etamax, nphibins,phimin,phimax);
 166   fhEtaPhiBCMinus->SetXTitle("#eta ");
 167   fhEtaPhiBCMinus->SetYTitle("#phi (rad)");
 168   outputContainer->Add(fhEtaPhiBCMinus);
 169
 170   fhEtaPhiBC0PileUpSPD  = new TH2F ("hEtaPhiBC0PileUpSPD","eta-phi for clusters tof corresponding to BC=0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
 171   fhEtaPhiBC0PileUpSPD->SetXTitle("#eta ");
 172   fhEtaPhiBC0PileUpSPD->SetYTitle("#phi (rad)");
 173   outputContainer->Add(fhEtaPhiBC0PileUpSPD);
 174
 175   fhEtaPhiBCPlusPileUpSPD  = new TH2F ("hEtaPhiBCPlusPileUpSPD","eta-phi for clusters tof corresponding to BC>0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
 176   fhEtaPhiBCPlusPileUpSPD->SetXTitle("#eta ");
 177   fhEtaPhiBCPlusPileUpSPD->SetYTitle("#phi (rad)");
 178   outputContainer->Add(fhEtaPhiBCPlusPileUpSPD);
 179
 180   fhEtaPhiBCMinusPileUpSPD  = new TH2F ("hEtaPhiBCMinusPileUpSPD","eta-phi for clusters tof corresponding to BC<0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
 181   fhEtaPhiBCMinusPileUpSPD->SetXTitle("#eta ");
 182   fhEtaPhiBCMinusPileUpSPD->SetYTitle("#phi (rad)");
 183   outputContainer->Add(fhEtaPhiBCMinusPileUpSPD);
 184
 185   TString title[] = {"no |t diff| cut","|t diff|<20 ns","|t diff|>20 ns","|t diff|>40 ns"};
 186   TString name [] = {"TDiffNoCut","TDiffSmaller20ns","TDiffLarger20ns","TDiffLarger40ns"};
 187   for(Int_t i = 0; i < 4; i++)
 188   {
 189     fhClusterMultSPDPileUp[i] = new TH2F(Form("fhClusterMultSPDPileUp_%s", name[i].Data()),
 190                                          Form("Number of clusters per pile up event with #it{E} > 0.5 and %s respect cluster max vs cluster max E ",title[i].Data()),
 191                                          nptbins,ptmin,ptimecluster,100,0,100);
 192     fhClusterMultSPDPileUp[i]->SetYTitle("n clusters ");
 193     fhClusterMultSPDPileUp[i]->SetXTitle("#it{E}_{cluster max} (GeV)");
 194     outputContainer->Add(fhClusterMultSPDPileUp[i]) ;
 195
 196     fhClusterMultNoPileUp[i] = new TH2F(Form("fhClusterMultNoPileUp_%s", name[i].Data()),
 197                                         Form("Number of clusters per non pile up event with #it{E} > 0.5 and %s respect cluster max vs cluster max E ",title[i].Data()),
 198                                         nptbins,ptmin,ptimecluster,100,0,100);
 199     fhClusterMultNoPileUp[i]->SetYTitle("n clusters ");
 200     fhClusterMultNoPileUp[i]->SetXTitle("#it{E}_{cluster max} (GeV)");
 201     outputContainer->Add(fhClusterMultNoPileUp[i]) ;
 202   }
 203
 204   fhPtNPileUpSPDVtx  = new TH2F ("hPt_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
 205                                  nptbins,ptmin,ptimecluster,20,0,20);
 206   fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
 207   fhPtNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 208   outputContainer->Add(fhPtNPileUpSPDVtx);
 209
 210   fhPtNPileUpTrkVtx  = new TH2F ("hPt_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
 211                                  nptbins,ptmin,ptimecluster, 20,0,20 );
 212   fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
 213   fhPtNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 214   outputContainer->Add(fhPtNPileUpTrkVtx);
 215
 216   fhPtNPileUpSPDVtxTimeCut  = new TH2F ("hPt_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
 217                                         nptbins,ptmin,ptimecluster,20,0,20);
 218   fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
 219   fhPtNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 220   outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
 221
 222   fhPtNPileUpTrkVtxTimeCut  = new TH2F ("hPt_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
 223                                         nptbins,ptmin,ptimecluster, 20,0,20 );
 224   fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
 225   fhPtNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 226   outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
 227
 228   fhPtNPileUpSPDVtxTimeCut2  = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
 229                                          nptbins,ptmin,ptimecluster,20,0,20);
 230   fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
 231   fhPtNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 232   outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
 233
 234   fhPtNPileUpTrkVtxTimeCut2  = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
 235                                          nptbins,ptmin,ptimecluster, 20,0,20 );
 236   fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
 237   fhPtNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 238   outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
 239
 240
 241   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
 242
 243   for(Int_t i = 0 ; i < 7 ; i++)
 244   {
 245     fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
 246                               Form("Cluster  #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptimecluster);
 247     fhPtPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 248     outputContainer->Add(fhPtPileUp[i]);
 249
 250     fhPtNeutralPileUp[i]  = new TH1F(Form("hPtNeutralPileUp%s",pileUpName[i].Data()),
 251                                      Form("Neutral clusters #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptimecluster);
 252     fhPtNeutralPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 253     outputContainer->Add(fhPtNeutralPileUp[i]);
 254
 255     fhClusterEFracLongTimePileUp[i]  = new TH2F(Form("hClusterEFracLongTimePileUp%s",pileUpName[i].Data()),
 256                                                 Form("Cluster E vs fraction of cluster energy from large T cells, %s Pile-Up event",pileUpName[i].Data()),
 257                                                 nptbins,ptmin,ptimecluster,200,0,1);
 258     fhClusterEFracLongTimePileUp[i]->SetXTitle("#it{E} (GeV)");
 259     fhClusterEFracLongTimePileUp[i]->SetYTitle("E(large time) / E");
 260     outputContainer->Add(fhClusterEFracLongTimePileUp[i]);
 261
 262     fhClusterCellTimePileUp[i]  = new TH2F(Form("hClusterCellTimePileUp%s",pileUpName[i].Data()),
 263                                            Form("Cluster E vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
 264                                            nptbins,ptmin,ptimecluster,ntimebins,timemin,timemax);
 265     fhClusterCellTimePileUp[i]->SetXTitle("#it{E} (GeV)");
 266     fhClusterCellTimePileUp[i]->SetYTitle("t_{cell} (ns)");
 267     outputContainer->Add(fhClusterCellTimePileUp[i]);
 268
 269     fhClusterTimeDiffPileUp[i]  = new TH2F(Form("hClusterTimeDiffPileUp%s",pileUpName[i].Data()),
 270                                            Form("Cluster E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
 271                                            nptbins,ptmin,ptimecluster,400,-200,200);
 272     fhClusterTimeDiffPileUp[i]->SetXTitle("#it{E} (GeV)");
 273     fhClusterTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
 274     outputContainer->Add(fhClusterTimeDiffPileUp[i]);
 275
 276     fhClusterTimeDiffNeutralPileUp[i]  = new TH2F(Form("hClusterTimeDiffNeutralPileUp%s",pileUpName[i].Data()),
 277                                                   Form("Neutral clusters E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
 278                                                   nptbins,ptmin,ptimecluster,400,-200,200);
 279     fhClusterTimeDiffNeutralPileUp[i]->SetXTitle("#it{E} (GeV)");
 280     fhClusterTimeDiffNeutralPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
 281     outputContainer->Add(fhClusterTimeDiffNeutralPileUp[i]);
 282
 283     fhLambda0PileUp[i]  = new TH2F(Form("hLambda0PileUp%s",pileUpName[i].Data()),
 284                                    Form("Cluster E vs #lambda^{2}_{0} in cluster, %s Pile-Up event",pileUpName[i].Data()),
 285                                    nptbins,ptmin,ptimecluster,ssbins,ssmin,ssmax);
 286     fhLambda0PileUp[i]->SetXTitle("#it{E} (GeV)");
 287     fhLambda0PileUp[i]->SetYTitle("#lambda^{2}_{0}");
 288     outputContainer->Add(fhLambda0PileUp[i]);
 289
 290     fhLambda0NeutralPileUp[i]  = new TH2F(Form("hLambda0NeutralPileUp%s",pileUpName[i].Data()),
 291                                           Form("Neutral clusters E vs #lambda^{2}_{0}in cluster, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptimecluster,ssbins,ssmin,ssmax);
 292     fhLambda0NeutralPileUp[i]->SetXTitle("#it{E} (GeV)");
 293     fhLambda0NeutralPileUp[i]->SetYTitle("#lambda^{2}_{0}");
 294     outputContainer->Add(fhLambda0NeutralPileUp[i]);
 295
 296   }
 297
 298   return outputContainer ;
 299
 300 }
 301
 302 //_______________________
 303 void AliAnaClusterPileUp::Init()
 304 {
 305   //Init
 306
 307   //Do some checks
 308   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn())
 309     AliFatal("You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!");
 310
 311   if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn())
 312     AliFatal("You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!");
 313
 314   if(GetReader()->GetDataType() == AliCaloTrackReader::kMC)
 315     AliFatal("You want to use MC data in analysis but this is not possible in pile-up!!");
 316
 317 }
 318
 319 //____________________________________________________________________________
 320 void AliAnaClusterPileUp::InitParameters()
 321 {
 322
 323   //Initialize the parameters of the analysis.
 324   AddToHistogramsName("AnaClusterPileUp_");
 325
 326   fCalorimeter = "EMCAL" ;
 327         fNCellsCut   = 2;
 328 }
 329
 330 //__________________________________________________________________
 331 void  AliAnaClusterPileUp::MakeAnalysisFillHistograms()
 332 {
 333   // Do cluster analysis
 334   // Remember to open time cuts in reader
 335
 336   //Select the calorimeter
 337   TObjArray * pl = 0x0;
 338   AliVCaloCells* cells    = 0;
 339   if      (fCalorimeter == "PHOS" )
 340   {
 341     pl    = GetPHOSClusters();
 342     cells = GetPHOSCells();
 343   }
 344   else if (fCalorimeter == "EMCAL")
 345   {
 346     pl    = GetEMCALClusters();
 347     cells = GetEMCALCells();
 348   }
 349
 350   if(!pl)
 351   {
 352     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 353     return;
 354   }
 355
 356   AliVEvent  * event = GetReader()->GetInputEvent();
 357   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
 358   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
 359
 360   //-------------------
 361   // N pile up vertices
 362   Int_t nVtxSPD = -1;
 363   Int_t nVtxTrk = -1;
 364
 365   if      (esdEv)
 366   {
 367                 nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
 368                 nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
 369   }//ESD
 370   else if (aodEv)
 371   {
 372                 nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
 373                 nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
 374   }//AOD
 375
 376   //-------------------
 377   // Loop on clusters
 378   Int_t nCaloClusters = pl->GetEntriesFast();
 379
 380   if(GetDebug() > 0) printf("AliAnaClusterPileUp::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
 381   //Init variables
 382   TLorentzVector mom;
 383   Int_t   idMax = 0;
 384   Float_t ptMax = 0;
 385   Float_t  tMax = 0;
 386
 387   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
 388   {
 389           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
 390     //printf("calo %d, %f\n",icalo,calo->E());
 391
 392     if(!calo)  continue; // it should not happen, but just in case
 393
 394     calo->GetMomentum(mom,GetVertex(0)) ;
 395
 396     Float_t  ecluster  = mom.E();
 397     Float_t ptcluster  = mom.Pt();
 398     Float_t l0cluster  = calo->GetM02();
 399     Float_t etacluster = mom.Eta();
 400     Float_t phicluster = mom.Phi();
 401     if(phicluster < 0) phicluster+=TMath::TwoPi();
 402     Float_t tofcluster   = calo->GetTOF()*1.e9;
 403
 404     Bool_t matched = IsTrackMatched(calo,GetReader()->GetInputEvent());
 405
 406     //.......................................
 407     //If too small or big energy, skip it
 408     if(ecluster < GetMinEnergy() || ecluster > GetMaxEnergy() ) continue ;
 409
 410     //.......................................
 411     if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
 412
 413      //.......................................
 414     //Check acceptance selection
 415     if(IsFiducialCutOn())
 416     {
 417       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 418       if(! in ) continue;
 419     }
 420
 421     // Select highest pt cluster passing the cuts
 422     if(ptcluster > ptMax && tofcluster < 30)
 423     {
 424       ptMax = ptcluster;
 425                         tMax  = tofcluster;
 426       idMax = icalo;
 427     }
 428
 429     //-------------------------------------
 430     // Cluster timing for different pile-up
 431
 432     fhTimePtNoCut->Fill(ptcluster,tofcluster);
 433     if(GetReader()->IsPileUpFromSPD()) fhTimePtSPD->Fill(ptcluster,tofcluster);
 434
 435     //----------------------------------------
 436     // correlate cluster and number of vertices
 437
 438     fhPtNPileUpSPDVtx->Fill(ptcluster,nVtxSPD);
 439                 fhPtNPileUpTrkVtx->Fill(ptcluster,nVtxTrk);
 440
 441                 if(TMath::Abs(tofcluster) < 30)
 442                 {
 443                         fhPtNPileUpSPDVtxTimeCut->Fill(ptcluster,nVtxSPD);
 444                         fhPtNPileUpTrkVtxTimeCut->Fill(ptcluster,nVtxTrk);
 445                 }
 446
 447     if(tofcluster < 75 && tofcluster > -30)
 448     {
 449       fhPtNPileUpSPDVtxTimeCut2->Fill(ptcluster,nVtxSPD);
 450       fhPtNPileUpTrkVtxTimeCut2->Fill(ptcluster,nVtxTrk);
 451     }
 452
 453     // Loop on the vertices arrays, correlate with timing
 454     // only for sufficiently large cluster energy
 455     if(ecluster > 8)
 456     {
 457       fhTimeNPileUpVertSPD  ->Fill(tofcluster,nVtxSPD);
 458       fhTimeNPileUpVertTrack->Fill(tofcluster,nVtxTrk);
 459
 460       Int_t ncont = -1;
 461       Float_t z1 = -1, z2 = -1;
 462       Float_t diamZ = -1;
 463       for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
 464       {
 465         if      (esdEv)
 466         {
 467           const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
 468           ncont=pv->GetNContributors();
 469           z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
 470           z2 = pv->GetZ();
 471           diamZ = esdEv->GetDiamondZ();
 472         }//ESD
 473         else if (aodEv)
 474         {
 475           AliAODVertex *pv=aodEv->GetVertex(iVert);
 476           if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
 477           ncont=pv->GetNContributors();
 478           z1=aodEv->GetPrimaryVertexSPD()->GetZ();
 479           z2=pv->GetZ();
 480           diamZ = aodEv->GetDiamondZ();
 481         }// AOD
 482
 483         Double_t distZ  = TMath::Abs(z2-z1);
 484         diamZ  = TMath::Abs(z2-diamZ);
 485
 486         fhTimeNPileUpVertContributors  ->Fill(tofcluster,ncont);
 487         fhTimePileUpMainVertexZDistance->Fill(tofcluster,distZ);
 488         fhTimePileUpMainVertexZDiamond ->Fill(tofcluster,diamZ);
 489
 490       }// vertex loop
 491     }
 492
 493     //------------------------------------
 494     // Eta-Phi cluster position depending on timing
 495     // Continue only for BC0
 496     if      (tofcluster > 28)
 497     {
 498       fhEtaPhiBCPlus ->Fill(etacluster,phicluster);
 499       if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBCPlusPileUpSPD ->Fill(etacluster,phicluster);
 500       continue;
 501     }
 502     else if (tofcluster <-28)
 503     {
 504       fhEtaPhiBCMinus->Fill(etacluster,phicluster);
 505       if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBCMinusPileUpSPD->Fill(etacluster,phicluster);
 506       continue ;
 507     }
 508
 509     //--------------------------------------
 510     // Fill histograms for clusters in BC=0
 511
 512     fhEtaPhiBC0->Fill(etacluster,phicluster); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBC0PileUpSPD    ->Fill(etacluster,phicluster);
 513
 514     if(GetReader()->IsPileUpFromSPD())               {fhPtPileUp[0]->Fill(ptcluster); fhLambda0PileUp[0]->Fill(ptcluster,l0cluster); }
 515     if(GetReader()->IsPileUpFromEMCal())             {fhPtPileUp[1]->Fill(ptcluster); fhLambda0PileUp[1]->Fill(ptcluster,l0cluster); }
 516     if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtPileUp[2]->Fill(ptcluster); fhLambda0PileUp[2]->Fill(ptcluster,l0cluster); }
 517     if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtPileUp[3]->Fill(ptcluster); fhLambda0PileUp[3]->Fill(ptcluster,l0cluster); }
 518     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtPileUp[4]->Fill(ptcluster); fhLambda0PileUp[4]->Fill(ptcluster,l0cluster); }
 519     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtPileUp[5]->Fill(ptcluster); fhLambda0PileUp[5]->Fill(ptcluster,l0cluster); }
 520     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtPileUp[6]->Fill(ptcluster); fhLambda0PileUp[6]->Fill(ptcluster,l0cluster); }
 521
 522     if(!matched)
 523     {
 524       if(GetReader()->IsPileUpFromSPD())               {fhPtNeutralPileUp[0]->Fill(ptcluster); fhLambda0NeutralPileUp[0]->Fill(ptcluster,l0cluster); }
 525       if(GetReader()->IsPileUpFromEMCal())             {fhPtNeutralPileUp[1]->Fill(ptcluster); fhLambda0NeutralPileUp[1]->Fill(ptcluster,l0cluster); }
 526       if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtNeutralPileUp[2]->Fill(ptcluster); fhLambda0NeutralPileUp[2]->Fill(ptcluster,l0cluster); }
 527       if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtNeutralPileUp[3]->Fill(ptcluster); fhLambda0NeutralPileUp[3]->Fill(ptcluster,l0cluster); }
 528       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtNeutralPileUp[4]->Fill(ptcluster); fhLambda0NeutralPileUp[4]->Fill(ptcluster,l0cluster); }
 529       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtNeutralPileUp[5]->Fill(ptcluster); fhLambda0NeutralPileUp[5]->Fill(ptcluster,l0cluster); }
 530       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtNeutralPileUp[6]->Fill(ptcluster); fhLambda0NeutralPileUp[6]->Fill(ptcluster,l0cluster); }
 531     }
 532
 533     //----------------------------------------------------------------------------
 534     // Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
 535     // Get the fraction of the cluster energy that carries the cell with highest
 536     // energy and its absId
 537
 538     Float_t maxCellFraction = 0.;
 539     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
 540
 541     Float_t clusterLongTimePt = 0;
 542     Float_t clusterOKTimePt   = 0;
 543
 544     if(cells->GetCellAmplitude(absIdMax) < 0.1) continue ;
 545
 546     for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
 547     {
 548       Int_t absId  = calo->GetCellsAbsId()[ipos];
 549
 550       if( absId == absIdMax ) continue ;
 551
 552       Double_t time  = cells->GetCellTime(absId);
 553       Float_t  amp   = cells->GetCellAmplitude(absId);
 554       Int_t    bc    = GetReader()->GetInputEvent()->GetBunchCrossNumber();
 555       GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,time,cells);
 556       time*=1e9;
 557
 558       Float_t diff = (tofcluster-time);
 559
 560       if(GetReader()->IsInTimeWindow(time,amp)) clusterOKTimePt   += amp;
 561       else                                      clusterLongTimePt += amp;
 562
 563       if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
 564
 565       if(GetReader()->IsPileUpFromSPD())
 566       {
 567         fhClusterCellTimePileUp[0]->Fill(ptcluster, time);
 568         fhClusterTimeDiffPileUp[0]->Fill(ptcluster, diff);
 569         if(!matched) fhClusterTimeDiffNeutralPileUp[0]->Fill(ptcluster, diff);
 570       }
 571
 572       if(GetReader()->IsPileUpFromEMCal())
 573       {
 574         fhClusterCellTimePileUp[1]->Fill(ptcluster, time);
 575         fhClusterTimeDiffPileUp[1]->Fill(ptcluster, diff);
 576         if(!matched) fhClusterTimeDiffNeutralPileUp[1]->Fill(ptcluster, diff);
 577       }
 578
 579       if(GetReader()->IsPileUpFromSPDOrEMCal())
 580       {
 581         fhClusterCellTimePileUp[2]->Fill(ptcluster, time);
 582         fhClusterTimeDiffPileUp[2]->Fill(ptcluster, diff);
 583         if(!matched) fhClusterTimeDiffNeutralPileUp[2]->Fill(ptcluster, diff);
 584       }
 585
 586       if(GetReader()->IsPileUpFromSPDAndEMCal())
 587       {
 588         fhClusterCellTimePileUp[3]->Fill(ptcluster, time);
 589         fhClusterTimeDiffPileUp[3]->Fill(ptcluster, diff);
 590         if(!matched) fhClusterTimeDiffNeutralPileUp[3]->Fill(ptcluster, diff);
 591       }
 592
 593       if(GetReader()->IsPileUpFromSPDAndNotEMCal())
 594       {
 595         fhClusterCellTimePileUp[4]->Fill(ptcluster, time);
 596         fhClusterTimeDiffPileUp[4]->Fill(ptcluster, diff);
 597         if(!matched) fhClusterTimeDiffNeutralPileUp[4]->Fill(ptcluster, diff);
 598       }
 599
 600       if(GetReader()->IsPileUpFromEMCalAndNotSPD())
 601       {
 602         fhClusterCellTimePileUp[5]->Fill(ptcluster, time);
 603         fhClusterTimeDiffPileUp[5]->Fill(ptcluster, diff);
 604         if(!matched) fhClusterTimeDiffNeutralPileUp[5]->Fill(ptcluster, diff);
 605       }
 606
 607       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
 608       {
 609         fhClusterCellTimePileUp[6]->Fill(ptcluster, time);
 610         fhClusterTimeDiffPileUp[6]->Fill(ptcluster, diff);
 611         if(!matched) fhClusterTimeDiffNeutralPileUp[6]->Fill(ptcluster, diff);
 612       }
 613     }//loop
 614
 615     Float_t frac = 0;
 616     if(clusterLongTimePt+clusterOKTimePt > 0.001)
 617       frac = clusterLongTimePt/(clusterLongTimePt+clusterOKTimePt);
 618     //printf("E long %f, E OK %f, Fraction large time %f, E %f\n",clusterLongTimePt,clusterOKTimePt,frac,ptcluster);
 619
 620     if(GetReader()->IsPileUpFromSPD())               fhClusterEFracLongTimePileUp[0]->Fill(ptcluster,frac);
 621     if(GetReader()->IsPileUpFromEMCal())             fhClusterEFracLongTimePileUp[1]->Fill(ptcluster,frac);
 622     if(GetReader()->IsPileUpFromSPDOrEMCal())        fhClusterEFracLongTimePileUp[2]->Fill(ptcluster,frac);
 623     if(GetReader()->IsPileUpFromSPDAndEMCal())       fhClusterEFracLongTimePileUp[3]->Fill(ptcluster,frac);
 624     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhClusterEFracLongTimePileUp[4]->Fill(ptcluster,frac);
 625     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhClusterEFracLongTimePileUp[5]->Fill(ptcluster,frac);
 626     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhClusterEFracLongTimePileUp[6]->Fill(ptcluster,frac);
 627   }//loop
 628
 629   //----------------------------------------------------------------------------------------------
 630   // Loop again on clusters to compare this max cluster t and the rest of the clusters, if E > 0.3
 631   Int_t n20  = 0;
 632   Int_t n40  = 0;
 633   Int_t n    = 0;
 634   Int_t nOK  = 0;
 635
 636   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
 637   {
 638           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
 639
 640     if(!calo || calo->E() < 0.3 || icalo == idMax) continue;
 641
 642     Float_t tdiff = TMath::Abs(tMax-calo->GetTOF()*1e9);
 643     n++;
 644     if(tdiff < 25) nOK++;
 645     else
 646     {
 647       n20++;
 648       if(tdiff > 40 ) n40++;
 649     }
 650   }
 651
 652   // Check pile-up and fill histograms depending on the different cluster multiplicities
 653   if(GetReader()->IsPileUpFromSPD())
 654   {
 655     fhClusterMultSPDPileUp[0]->Fill(ptMax,n  );
 656     fhClusterMultSPDPileUp[1]->Fill(ptMax,nOK);
 657     fhClusterMultSPDPileUp[2]->Fill(ptMax,n20);
 658     fhClusterMultSPDPileUp[3]->Fill(ptMax,n40);
 659   }
 660   else
 661   {
 662     fhClusterMultNoPileUp[0]->Fill(ptMax,n  );
 663     fhClusterMultNoPileUp[1]->Fill(ptMax,nOK);
 664     fhClusterMultNoPileUp[2]->Fill(ptMax,n20);
 665     fhClusterMultNoPileUp[3]->Fill(ptMax,n40);
 666   }
 667
 668
 669   if(GetDebug() > 1) printf("AliAnaClusterPileUp::MakeAnalysisFillHistograms()  End fill histograms\n");
 670
 671 }
 672
 673
 674
 675 //__________________________________________________________________
 676 void AliAnaClusterPileUp::Print(const Option_t * opt) const
 677 {
 678   //Print some relevant parameters set for the analysis
 679
 680   if(! opt)
 681     return;
 682
 683   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 684   AliAnaCaloTrackCorrBaseClass::Print(" ");
 685
 686   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
 687   printf("    \n") ;
 688
 689 }