PWGGA/CaloTrackCorrelations/AliAnaRandomTrigger.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
  16 //_________________________________________________________________________
  17 // Gerenate a random trigger, input for other analysis
  18 // Set flat energy distribution over acceptance of EMCAL, PHOS or CTS
  19 // Be careful, correlate only with Min Bias events this trigger
  20 //
  21 //
  22 //-- Author: Gustavo Conesa (LPSC-Grenoble)
  23 //_________________________________________________________________________
  24
  25
  26 // --- ROOT system ---
  27 #include <TH2F.h>
  28 #include <TClonesArray.h>
  29
  30 //---- AliRoot system ----
  31 #include "AliAnaRandomTrigger.h"
  32 #include "AliAODPWG4ParticleCorrelation.h"
  33 #include "AliEMCALGeometry.h"
  34
  35 ClassImp(AliAnaRandomTrigger)
  36
  37 //__________________________________________
  38 AliAnaRandomTrigger::AliAnaRandomTrigger() :
  39     AliAnaCaloTrackCorrBaseClass(),
  40     fDetector("EMCAL"), fRandom(0), fNRandom(0),
  41     fhE(0),             fhPt(0),
  42     fhPhi(0),           fhEta(0),
  43     fhEtaPhi(0)
  44 {
  45   //Default Ctor
  46
  47   //Initialize parameters
  48   InitParameters();
  49
  50 }
  51
  52 //_________________________________________________________________________
  53 Bool_t AliAnaRandomTrigger::ExcludeDeadBadRegions(Float_t eta, Float_t phi)
  54 {
  55   // Check if there is a dead or bad region in a detector
  56   // Now only EMCAL
  57
  58   if(fDetector!="EMCAL") return kFALSE;
  59
  60   //-------------------------------------
  61   // Get the corresponding cell in EMCAL, check if it exists in acceptance (phi gaps, borders)
  62   //-------------------------------------
  63
  64   Int_t absId = -1;
  65   if(!GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta,phi, absId)) return kTRUE; // remove if out of EMCAL acceptance, phi gaps
  66
  67   Int_t icol = -1, irow = -1, iRCU = -1;
  68   Int_t sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId,"EMCAL", icol, irow, iRCU);
  69
  70   //printf("eta %f, phi %f, ieta %d, iphi %d, sm %d\n",eta,phi,icol,irow,sm);
  71
  72   //-------------------------------------
  73   // Remove in case of close to border, by default always 1 but check what was set in reco utils
  74   //-------------------------------------
  75
  76   Bool_t okrow = kFALSE;
  77         Bool_t okcol = kFALSE;
  78   Int_t nborder = GetCaloUtils()->GetEMCALRecoUtils()->GetNumberOfCellsFromEMCALBorder();
  79   if (nborder<1) nborder = 1;
  80
  81   // Rows
  82   if(sm < 10)
  83   {
  84     if(irow >= nborder && irow < 24-nborder) okrow =kTRUE;
  85   }
  86   else
  87   {
  88     if((GetCaloUtils()->EMCALGeometryName()).Contains("12SM")) // 1/3 SM
  89     {
  90       if(irow >= nborder && irow < 8-nborder) okrow =kTRUE;
  91     }
  92     else // 1/2 SM
  93     {
  94       if(irow >= nborder && irow <12-nborder) okrow =kTRUE;
  95     }
  96   }
  97
  98   // Columns
  99   if(sm%2==0)
 100   {
 101     if(icol >= nborder)     okcol = kTRUE;
 102   }
 103   else
 104   {
 105     if(icol <  48-nborder)  okcol = kTRUE;
 106   }
 107
 108   //printf("okcol %d, okrow %d\n",okcol,okrow);
 109   if (!okcol || !okrow) return kTRUE;
 110
 111   //-------------------------------------
 112   // Check if the cell or those around are bad
 113   //-------------------------------------
 114
 115   if(GetCaloUtils()->GetEMCALChannelStatus(sm,icol, irow) > 0) return kTRUE ; // trigger falls into a bad channel
 116
 117   // Check if close there was a bad channel
 118 //  for(Int_t i = -1; i <= 1; i++)
 119 //  {
 120 //    for(Int_t j = -1; j <= 1; j++)
 121 //    {
 122 //      //printf("\t check icol %d, irow %d \n",icol+i, irow+j);
 123 //      if(GetCaloUtils()->GetEMCALChannelStatus(sm,icol+i, irow+j) > 0) return kTRUE ; // trigger falls into a bad channel
 124 //      //printf("\t ok\n");
 125 //    }
 126 //  }
 127
 128    //printf("\t OK\n");
 129
 130   return kFALSE;
 131
 132 }
 133
 134
 135 //__________________________________________________
 136 TObjString *  AliAnaRandomTrigger::GetAnalysisCuts()
 137 {
 138   //Save parameters used for analysis
 139   TString parList ; //this will be list of parameters used for this analysis.
 140   const Int_t buffersize = 255;
 141   char onePar[buffersize] ;
 142
 143   snprintf(onePar,buffersize,"--- AliAnaRandomTrigger ---\n") ;
 144   parList+=onePar ;
 145   snprintf(onePar,buffersize,"Detector: %s\n"    , fDetector.Data()) ;
 146   parList+=onePar ;
 147   snprintf(onePar,buffersize,"N per event = %d\n", fNRandom       ) ;
 148   parList+=onePar ;
 149   snprintf(onePar,buffersize,"Min E   = %3.2f - Max E   = %3.2f\n", GetMinPt(), GetMaxPt()) ;
 150   parList+=onePar ;
 151   snprintf(onePar,buffersize,"Min Eta = %3.2f - Max Eta = %3.2f\n", fEtaCut[0], fEtaCut[1]) ;
 152   parList+=onePar ;
 153   snprintf(onePar,buffersize,"Min Phi = %3.2f - Max Phi = %3.2f\n", fPhiCut[0], fPhiCut[1]) ;
 154   parList+=onePar ;
 155
 156   return new TObjString(parList) ;
 157 }
 158
 159 //_______________________________________________________
 160 TList *  AliAnaRandomTrigger::GetCreateOutputObjects()
 161 {
 162   // Create histograms to be saved in output file and
 163   // store them in fOutputContainer
 164
 165   TList * outputContainer = new TList() ;
 166   outputContainer->SetName("RandomTrigger") ;
 167
 168   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins(); Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
 169   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
 170   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 171
 172   fhE  = new TH1F ("hE","Random E distribution", nptbins,ptmin,ptmax);
 173   fhE->SetXTitle("E (GeV)");
 174   outputContainer->Add(fhE);
 175
 176   fhPt  = new TH1F ("hPt","Random p_{T} distribution", nptbins,ptmin,ptmax);
 177   fhPt->SetXTitle("p_{T} (GeV/c)");
 178   outputContainer->Add(fhPt);
 179
 180   fhPhi  = new TH2F ("hPhi","Random #phi distribution",
 181                         nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 182   fhPhi->SetYTitle("#phi (rad)");
 183   fhPhi->SetXTitle("p_{T} (GeV/c)");
 184   outputContainer->Add(fhPhi);
 185
 186   fhEta  = new TH2F ("hEta","Random #eta distribution",
 187                         nptbins,ptmin,ptmax, netabins,etamin,etamax);
 188   fhEta->SetYTitle("#eta ");
 189   fhEta->SetXTitle("p_{T} (GeV/c)");
 190   outputContainer->Add(fhEta);
 191
 192   fhEtaPhi  = new TH2F ("hEtaPhi","Random #eta vs #phi ",netabins,etamin,etamax, nphibins,phimin,phimax);
 193   fhEtaPhi->SetXTitle("#eta ");
 194   fhEtaPhi->SetYTitle("#phi (rad)");
 195   outputContainer->Add(fhEtaPhi);
 196
 197   return outputContainer;
 198
 199 }
 200
 201 //___________________________________________
 202 void AliAnaRandomTrigger::InitParameters()
 203 {
 204   //Initialize the parameters of the analysis.
 205   SetOutputAODClassName("AliAODPWG4ParticleCorrelation");
 206   SetOutputAODName("RandomTrigger");
 207
 208   AddToHistogramsName("AnaRandomTrigger_");
 209
 210   fNRandom   = 1    ;
 211   fPhiCut[0] = 0.   ;
 212   fPhiCut[1] = TMath::TwoPi() ;
 213   fEtaCut[0] =-1.   ;
 214   fEtaCut[1] = 1.   ;
 215
 216 }
 217
 218 //____________________________________________________________
 219 void AliAnaRandomTrigger::Print(const Option_t * opt) const
 220 {
 221   //Print some relevant parameters set for the analysis
 222   if(! opt)
 223     return;
 224
 225   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 226   AliAnaCaloTrackCorrBaseClass::Print(" ");
 227
 228   printf("Detector = %s\n",  fDetector.Data());
 229   printf("Min E   = %3.2f - Max E   = %3.2f\n", GetMinPt(), GetMaxPt());
 230   printf("Min Eta = %3.2f - Max Eta = %3.2f\n", fEtaCut[0], fEtaCut[1]);
 231   printf("Min Phi = %3.2f - Max Phi = %3.2f\n", fPhiCut[0], fPhiCut[1]);
 232
 233 }
 234
 235 //______________________________________________
 236 void  AliAnaRandomTrigger::MakeAnalysisFillAOD()
 237 {
 238   // Do analysis and fill aods
 239   // Generate particle randomly
 240   // fNRandom particles per event
 241
 242   for(Int_t irandom = 0; irandom < fNRandom; irandom++)
 243   {
 244     // Get the random variables of the trigger
 245     Float_t pt  = fRandom.Uniform(GetMinPt(), GetMaxPt());
 246     Float_t eta = fRandom.Uniform(fEtaCut[0], fEtaCut[1]);
 247     Float_t phi = fRandom.Uniform(fPhiCut[0], fPhiCut[1]);
 248
 249     // Check if particle falls into a dead region, if inside, get new
 250     Bool_t excluded =  ExcludeDeadBadRegions(eta,phi);
 251
 252     // if excluded, generate a new trigger until accepted
 253     while (excluded)
 254     {
 255       pt  = fRandom.Uniform(GetMinPt(), GetMaxPt());
 256       eta = fRandom.Uniform(fEtaCut[0], fEtaCut[1]);
 257       phi = fRandom.Uniform(fPhiCut[0], fPhiCut[1]);
 258
 259       excluded = ExcludeDeadBadRegions(eta,phi);
 260     }
 261
 262     // Create the AOD trigger object
 263     TLorentzVector mom;
 264     mom.SetPtEtaPhiM(pt,eta,phi,0);
 265
 266     AliAODPWG4Particle trigger = AliAODPWG4Particle(mom);
 267     trigger.SetDetector(fDetector);
 268
 269     if(GetDebug() > 1)
 270       printf("AliAnaRandomTrigger::MakeAnalysisFillAOD() -  iRandom %d, Trigger e %2.2f pt %2.2f, phi %2.2f, eta %2.2f \n",
 271              irandom, trigger.E(), trigger.Pt(), trigger.Phi(), trigger.Eta());
 272
 273     AddAODParticle(trigger);
 274   }
 275
 276   if(GetDebug() > 0)
 277     printf("AliAnaRandomTrigger::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
 278 }
 279
 280 //_____________________________________________________
 281 void  AliAnaRandomTrigger::MakeAnalysisFillHistograms()
 282 {
 283   // Fill histograms
 284
 285   //Loop on stored AODParticles
 286   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 287
 288   if(GetDebug() > 0)
 289     printf("AliAnaRandomTrigger::MakeAnalysisFillHistograms() - aod branch entries %d, fNRandom %d\n", naod, fNRandom);
 290
 291   for(Int_t iaod = 0; iaod < naod ; iaod++)
 292   {
 293     AliAODPWG4Particle* trigger =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 294
 295     fhPt    ->Fill(trigger->Pt());
 296     fhE     ->Fill(trigger->E());
 297     fhPhi   ->Fill(trigger->Pt(), trigger->Phi());
 298     fhEta   ->Fill(trigger->Pt(), trigger->Eta());
 299     fhEtaPhi->Fill(trigger->Eta(),trigger->Phi());
 300
 301   }// aod branch loop
 302
 303 }