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