PWG4/PartCorrDep/AliAnaShowerParameter.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 /* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
  16
  17 //_________________________________________________________________________
  18 //
  19 // Class cloned from AliAnaPhoton, main aim is shower shape studies
  20 //
  21 //
  22 //
  23 //-- Author: Jocelyn Mlynarz (WSU) and Gustavo Conesa (LPSC-CNRS)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 #include <TH3F.h>
  29 #include <TH2F.h>
  30 #include <TClonesArray.h>
  31 //#include <TObjString.h>
  32 #include <Riostream.h>
  33 #include "TParticle.h"
  34 //#include <fstream>
  35
  36 // --- Analysis system ---
  37 #include "AliAnaShowerParameter.h"
  38 #include "AliCaloTrackReader.h"
  39 #include "AliStack.h"
  40 #include "AliCaloPID.h"
  41 #include "AliMCAnalysisUtils.h"
  42 #include "AliFiducialCut.h"
  43 #include "AliAODCaloCluster.h"
  44 #include "AliAODMCParticle.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliAnalysisTaskParticleCorrelation.h"
  47 #include "AliAODEvent.h"
  48
  49
  50 ClassImp(AliAnaShowerParameter)
  51
  52 //____________________________________________________________________________
  53 AliAnaShowerParameter::AliAnaShowerParameter() :
  54 AliAnaPartCorrBaseClass(), fCalorimeter(""), fNCellsCutMin(0),
  55 fNCellsCutMax(0), fLambdaCut(0), fTimeCutMin(-1), fTimeCutMax(9999999),
  56 fhNClusters(0), fhNCellCluster(0), fhEtaPhiPtCluster(0),
  57 fhLambdaPtCluster(0),
  58
  59 //MC
  60
  61 fhLambdaPtPhoton(0), fhLambdaPtPi0(0), fhLambdaPtPion(0), fhPtTruthPi0(0)
  62
  63 {
  64   //default ctor
  65
  66   //Initialize parameters
  67   InitParameters();
  68
  69 }
  70
  71 //____________________________________________________________________________
  72 AliAnaShowerParameter::~AliAnaShowerParameter()
  73 {
  74   //dtor
  75
  76 }
  77
  78 //________________________________________________________________________
  79 TObjString *  AliAnaShowerParameter::GetAnalysisCuts()
  80 {
  81   //Save parameters used for analysis
  82   TString parList ; //this will be list of parameters used for this analysis.
  83   const Int_t buffersize = 255;
  84   char onePar[buffersize] ;
  85
  86   snprintf(onePar,buffersize,"--- AliAnaShowerParameter ---\n") ;
  87   parList+=onePar ;
  88   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
  89   parList+=onePar ;
  90   snprintf(onePar,buffersize,"fNCellsCutMin: (Cut in the minimum number of cells) %e\n",fNCellsCutMin) ;
  91   parList+=onePar ;
  92   snprintf(onePar,buffersize,"fNCellsCutMax: (Cut in the maximum number of cells) %e\n",fNCellsCutMax) ;
  93   parList+=onePar ;
  94   snprintf(onePar,buffersize,"fLambdaCut: (Cut in the minimum lambda) %e\n",fLambdaCut) ;
  95   parList+=onePar ;
  96
  97   //Get parameters set in base class.
  98   parList += GetBaseParametersList() ;
  99
 100   //Get parameters set in PID class.
 101   parList += GetCaloPID()->GetPIDParametersList() ;
 102
 103   //Get parameters set in FiducialCut class (not available yet)
 104   //parlist += GetFidCut()->GetFidCutParametersList()
 105
 106   return new TObjString(parList) ;
 107 }
 108
 109
 110 //________________________________________________________________________
 111 TList *  AliAnaShowerParameter::GetCreateOutputObjects()
 112 {
 113   // Create histograms to be saved in output file and
 114   // store them in outputContainer
 115   TList * outputContainer = new TList() ;
 116   outputContainer->SetName("PhotonHistos") ;
 117
 118   Int_t nptbins  = GetHistoPtBins();
 119   Int_t nphibins = GetHistoPhiBins();
 120   Int_t netabins = GetHistoEtaBins();
 121   Float_t ptmax  = GetHistoPtMax();
 122   Float_t phimax = GetHistoPhiMax();
 123   Float_t etamax = GetHistoEtaMax();
 124   Float_t ptmin  = GetHistoPtMin();
 125   Float_t phimin = GetHistoPhiMin();
 126   Float_t etamin = GetHistoEtaMin();
 127
 128   //General non-MC Cluster histograms
 129   fhNClusters = new TH1F ("hNClusters","NClusters",21,-0.5,20.5);
 130   fhNClusters->SetXTitle("N_{Clusters}");
 131   outputContainer->Add(fhNClusters) ;
 132
 133   fhNCellCluster = new TH2F ("hNCellCluster","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5);
 134   fhNCellCluster->SetYTitle("N_{Cells}");
 135   fhNCellCluster->SetXTitle("p_{T}");
 136   outputContainer->Add(fhNCellCluster) ;
 137
 138   fhEtaPhiPtCluster = new TH3F
 139   ("hEtaPhiPtCluster","#phi_{Cluster} and #eta_{Cluster}",nptbins,ptmin,ptmax,nphibins,phimin,phimax,netabins,etamin,etamax);
 140   fhEtaPhiPtCluster->SetZTitle("#eta");
 141   fhEtaPhiPtCluster->SetYTitle("#phi");
 142   fhEtaPhiPtCluster->SetXTitle("pT_{Cluster} (GeV/c)");
 143   outputContainer->Add(fhEtaPhiPtCluster) ;
 144
 145   fhLambdaPtCluster  = new TH2F
 146   ("hLambdaCluster","#lambda_{Cluster}",nptbins,ptmin,ptmax,300,0,3);
 147   fhLambdaPtCluster->SetYTitle("#lambda_{0}^{2}");
 148   fhLambdaPtCluster->SetXTitle("p_{T Cluster} (GeV/c)");
 149   outputContainer->Add(fhLambdaPtCluster) ;
 150
 151   if(IsDataMC()){
 152
 153     fhLambdaPtPhoton  = new TH2F
 154     ("hLambdaPtPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2);
 155     fhLambdaPtPhoton->SetYTitle("#lambda_{0}^{2}");
 156     fhLambdaPtPhoton->SetXTitle("pT_{#gamma, Reco} (GeV)");
 157     outputContainer->Add(fhLambdaPtPhoton) ;
 158
 159     fhLambdaPtPi0  = new TH2F
 160     ("hLambdaPtPi0","#lambda_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,2);
 161     fhLambdaPtPi0->SetYTitle("#lambda_{0}^{2}");
 162     fhLambdaPtPi0->SetXTitle("pT_{#pi^{0}, Reco} (GeV)");
 163     outputContainer->Add(fhLambdaPtPi0) ;
 164
 165     fhLambdaPtPion  = new TH2F
 166     ("hLambdaPtPion","#lambda_{#pi^{+}}",nptbins,ptmin,ptmax,200,0,2);
 167     fhLambdaPtPion->SetYTitle("#lambda_{0}^{2}");
 168     fhLambdaPtPion->SetXTitle("pT_{#pi^{+}, Reco} (GeV)");
 169     outputContainer->Add(fhLambdaPtPion) ;
 170
 171     fhPtTruthPi0  = new TH1D("hPtTruthPi0","#pi^{0} MC truth pT",nptbins,ptmin,ptmax) ;
 172     fhPtTruthPi0->SetXTitle("pT_{#pi^{0}, Reco} (GeV)");
 173     outputContainer->Add(fhPtTruthPi0) ;
 174
 175   }//Histos with MC
 176
 177   return outputContainer ;
 178
 179 }
 180
 181 //____________________________________________________________________________
 182 void AliAnaShowerParameter::Init()
 183 {
 184
 185   //Init
 186   //Do some checks
 187   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 188     printf("AliAnaShowerParameter::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 189     abort();
 190   }
 191   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 192     printf("AliAnaShowerParameter::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 193     abort();
 194   }
 195
 196 }
 197
 198 //____________________________________________________________________________
 199 void AliAnaShowerParameter::InitParameters()
 200 {
 201
 202   //Initialize the parameters of the analysis.
 203   AddToHistogramsName("AnaLambda_");
 204
 205   fCalorimeter = "EMCAL" ;
 206
 207   fTimeCutMin  = -1;
 208   fTimeCutMax  = 9999999;
 209   fNCellsCutMin = 0 ;
 210   fNCellsCutMax = 0 ;
 211   fLambdaCut = 0.01 ;
 212
 213 }
 214
 215 //__________________________________________________________________
 216 void  AliAnaShowerParameter::MakeAnalysisFillHistograms()
 217 {
 218
 219   //Do analysis and fill histograms
 220   if ( GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Starting analysis.\n");
 221
 222   // Access MC information in stack if requested, check that it exists.
 223   AliStack * stack = 0x0;
 224   TClonesArray * mcparticles0 = 0x0;
 225   Int_t NClusters = 0 ;
 226   TLorentzVector momCluster ;
 227
 228   //Check if the stack is available when analysing MC data.
 229   if(IsDataMC()){
 230
 231     if(GetReader()->ReadStack()){
 232       stack =  GetMCStack() ;
 233       if(!stack) {
 234         printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
 235                                 abort();
 236       }
 237
 238     }
 239     else if(GetReader()->ReadAODMCParticles()){
 240
 241       //Get the list of MC particles
 242       mcparticles0 = GetReader()->GetAODMCParticles(0);
 243       if(!mcparticles0 && GetDebug() > 0)       {
 244         printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() -  Standard MCParticles not available !\n");
 245       }
 246     }// is data and MC
 247   }
 248   //Loop on stored AOD photons
 249
 250   TClonesArray *  clustArray = 0x0 ;
 251   clustArray = GetAODCaloClusters() ;
 252   NClusters = clustArray->GetEntriesFast() ;
 253   fhNClusters->Fill(NClusters) ;
 254
 255   if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - aod branch entries %d\n", NClusters);
 256
 257   for(Int_t iCluster = 0 ; iCluster < NClusters ; iCluster++){
 258     Int_t input = 0 ;
 259
 260     AliAODCaloCluster * caloCluster =  (AliAODCaloCluster*) (clustArray->At(iCluster)) ;
 261     caloCluster->GetMomentum(momCluster,GetVertex(0)) ;
 262     AliAODPWG4Particle * aodCluster = new AliAODPWG4Particle(momCluster) ;
 263     Int_t labelCluster = caloCluster->GetLabel() ;
 264     aodCluster->SetLabel(labelCluster) ;;
 265     aodCluster->SetInputFileIndex(input) ;
 266     aodCluster->SetCaloLabel(caloCluster->GetID(),-1) ;
 267     aodCluster->SetDetector(fCalorimeter) ;
 268     aodCluster->SetTag(GetMCAnalysisUtils()->CheckOrigin(caloCluster->GetLabels(),caloCluster->GetNLabels(),GetReader(), aodCluster->GetInputFileIndex())) ;
 269
 270     if (GetDebug() > 3) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster particle label = %d\n",labelCluster) ;
 271
 272     //Fill Cluster histograms
 273     Float_t ptcluster  = aodCluster->Pt() ;
 274     Float_t phicluster = aodCluster->Phi() ;
 275     Float_t etacluster = aodCluster->Eta() ;
 276     Float_t lambdaMainCluster = caloCluster->GetM02() ;
 277     Int_t iNumCell = caloCluster->GetNCells() ;
 278
 279     if (iNumCell>=fNCellsCutMin&&iNumCell<=fNCellsCutMax&&lambdaMainCluster>=fLambdaCut){
 280
 281       //Fill the basic non-MC information about the cluster.
 282       fhNCellCluster->Fill(ptcluster,iNumCell) ;
 283       fhEtaPhiPtCluster->Fill(ptcluster,phicluster,etacluster) ;
 284       fhLambdaPtCluster->Fill(ptcluster,lambdaMainCluster) ;
 285
 286       //Play with the MC data if available
 287       if(IsDataMC()){
 288
 289         //Get the tag from AliMCAnalysisUtils for PID
 290         Int_t tag = aodCluster->GetTag();
 291         if (GetDebug() > 3) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster particle tag= %d\n",tag) ;
 292         if ( aodCluster->GetLabel() < 0){
 293           if(GetDebug() > 0)
 294             printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Label is -1; problem in the MC ESD? ");
 295           continue;
 296         }
 297
 298         //Check if the tag matches one of the different particle types and fill the corresponding histograms
 299         if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)&&!(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))) {
 300           fhLambdaPtPhoton->Fill(ptcluster,lambdaMainCluster) ;
 301         }//kMCPhoton
 302
 303         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)) {
 304           fhLambdaPtPi0->Fill(ptcluster,lambdaMainCluster) ;
 305         }
 306         if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion)) {
 307           fhLambdaPtPion->Fill(ptcluster,lambdaMainCluster) ;
 308         }
 309
 310         // Access MC information in stack if requested, check that it exists.
 311         Int_t label = aodCluster->GetLabel();
 312         if(label < 0) {
 313           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
 314           continue;
 315         }
 316       }
 317     }
 318   }
 319
 320   Float_t fPtGener = 0 ;
 321   if(IsDataMC() && stack){
 322     TParticle* pGener = stack->Particle(0) ;
 323     fPtGener = pGener->Pt() ;
 324   }
 325
 326
 327 }
 328
 329
 330 //__________________________________________________________________
 331 void AliAnaShowerParameter::Print(const Option_t * opt) const
 332 {
 333   //Print some relevant parameters set for the analysis
 334
 335   if(! opt)
 336     return;
 337
 338   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 339   AliAnaPartCorrBaseClass::Print(" ");
 340   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
 341   printf("Min number of cells in cluster is        > %f \n", fNCellsCutMin);
 342   printf("Max number of cells in cluster is        > %f \n", fNCellsCutMax);
 343   printf("Min lambda of cluster is        > %f \n", fLambdaCut);
 344   printf("    \n") ;
 345
 346 }