New class to generate randomly trigger particles for isolation and correlation analysis.
[u/mrichter/AliRoot.git] / 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),
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(const Float_t eta, const 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,"Min E   = %3.2f - Max E   = %3.2f\n", GetMinPt(), GetMaxPt()) ;
148   parList+=onePar ;
149   snprintf(onePar,buffersize,"Min Eta = %3.2f - Max Eta = %3.2f\n", fEtaCut[0], fEtaCut[1]) ;
150   parList+=onePar ;
151   snprintf(onePar,buffersize,"Min Phi = %3.2f - Max Phi = %3.2f\n", fPhiCut[0], fPhiCut[1]) ;
152   parList+=onePar ;
153    
154   return new TObjString(parList) ;
155 }
156
157 //_______________________________________________________
158 TList *  AliAnaRandomTrigger::GetCreateOutputObjects()
159 {  
160   // Create histograms to be saved in output file and 
161   // store them in fOutputContainer
162   
163   
164   TList * outputContainer = new TList() ; 
165   outputContainer->SetName("RandomTrigger") ; 
166   
167   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins(); Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
168   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
169   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();   
170
171   fhE  = new TH1F ("hE","E distribution", nptbins,ptmin,ptmax); 
172   fhE->SetXTitle("E (GeV)");
173   outputContainer->Add(fhE);
174   
175   fhPt  = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax); 
176   fhPt->SetXTitle("p_{T} (GeV/c)");
177   outputContainer->Add(fhPt);
178   
179   fhPhi  = new TH2F ("hPhi","#phi distribution",
180                         nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
181   fhPhi->SetYTitle("#phi (rad)");
182   fhPhi->SetXTitle("p_{T} (GeV/c)");
183   outputContainer->Add(fhPhi);
184   
185   fhEta  = new TH2F ("hEta","#eta distribution",
186                         nptbins,ptmin,ptmax, netabins,etamin,etamax); 
187   fhEta->SetYTitle("#eta ");
188   fhEta->SetXTitle("p_{T} (GeV/c)");
189   outputContainer->Add(fhEta);
190   
191   fhEtaPhi  = new TH2F ("hEtaPhi","pt/eta/phi of positive charge",netabins,etamin,etamax, nphibins,phimin,phimax); 
192   fhEtaPhi->SetXTitle("#eta ");
193   fhEtaPhi->SetYTitle("#phi (rad)");  
194   outputContainer->Add(fhEtaPhi);
195     
196   return outputContainer;
197
198 }
199
200 //___________________________________________
201 void AliAnaRandomTrigger::InitParameters()
202
203   //Initialize the parameters of the analysis.
204   SetOutputAODClassName("AliAODPWG4ParticleCorrelation");
205   SetOutputAODName("RandomTrigger");
206
207   AddToHistogramsName("AnaRandomTrigger_");
208   
209   fPhiCut[0] = 0.   ;
210   fPhiCut[1] = TMath::TwoPi() ;
211   fEtaCut[0] =-1.   ;
212   fEtaCut[1] = 1.   ;
213   
214 }
215
216 //____________________________________________________________
217 void AliAnaRandomTrigger::Print(const Option_t * opt) const
218 {
219   //Print some relevant parameters set for the analysis
220   if(! opt)
221     return;
222   
223   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
224   AliAnaCaloTrackCorrBaseClass::Print(" ");     
225
226   printf("Detector = %s\n",  fDetector.Data());
227   printf("Min E   = %3.2f - Max E   = %3.2f\n", GetMinPt(), GetMaxPt());
228   printf("Min Eta = %3.2f - Max Eta = %3.2f\n", fEtaCut[0], fEtaCut[1]);
229   printf("Min Phi = %3.2f - Max Phi = %3.2f\n", fPhiCut[0], fPhiCut[1]);
230
231
232
233 //______________________________________________
234 void  AliAnaRandomTrigger::MakeAnalysisFillAOD() 
235 {
236   // Do analysis and fill aods
237   // Generate particle randomly
238     
239   // Get the random variables of the trigger
240   Float_t pt  = fRandom.Uniform(GetMinPt(), GetMaxPt());
241   Float_t eta = fRandom.Uniform(fEtaCut[0], fEtaCut[1]);
242   Float_t phi = fRandom.Uniform(fPhiCut[0], fPhiCut[1]);
243     
244   // Check if particle falls into a dead region, if inside, get new
245   Bool_t excluded =  ExcludeDeadBadRegions(eta,phi);
246   
247   // if excluded, generate a new trigger until accepted
248   while (excluded)
249   {
250     pt  = fRandom.Uniform(GetMinPt(), GetMaxPt());
251     eta = fRandom.Uniform(fEtaCut[0], fEtaCut[1]);
252     phi = fRandom.Uniform(fPhiCut[0], fPhiCut[1]);
253     
254     excluded = ExcludeDeadBadRegions(eta,phi);
255   }
256   
257   // Create the AOD trigger object
258   TLorentzVector mom;
259   mom.SetPtEtaPhiM(pt,eta,phi,0);
260   
261   AliAODPWG4Particle trigger = AliAODPWG4Particle(mom);
262   trigger.SetDetector(fDetector);
263   
264   if(GetDebug() > 1) 
265     printf("AliAnaRandomTrigger::MakeAnalysisFillAOD() - Trigger e %2.2f pt %2.2f, phi %2.2f, eta %2.2f \n",
266            trigger.E(), trigger.Pt(), trigger.Phi(), trigger.Eta());
267   
268   AddAODParticle(trigger);
269   
270   if(GetDebug() > 0)    
271     printf("AliAnaRandomTrigger::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());   
272
273
274 //_____________________________________________________
275 void  AliAnaRandomTrigger::MakeAnalysisFillHistograms() 
276 {
277   // Fill histograms
278   
279   //Loop on stored AODParticles
280   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
281   
282   if(GetDebug() > 0) 
283     printf("AliAnaRandomTrigger::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
284   
285   for(Int_t iaod = 0; iaod < naod ; iaod++){
286     AliAODPWG4Particle* trigger =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
287     
288     fhPt    ->Fill(trigger->Pt());
289     fhE     ->Fill(trigger->E());
290     fhPhi   ->Fill(trigger->Pt(), trigger->Phi());
291     fhEta   ->Fill(trigger->Pt(), trigger->Eta());
292     fhEtaPhi->Fill(trigger->Eta(),trigger->Phi());
293     
294   }// aod branch loop
295   
296 }