1)Terminate() method implemented in the frame. Simple examples on what to do with...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaCaloTriggerMC.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 //_________________________________________________________________________
18 // An analysis task to check the trigger data in ESD with MC data
19 // Creates an ntuple for 2x2 and NxN triggers
20 // Each ntuple connects the maximum trigger amplitudes 
21 // and its positions with reconstructed clusters and MC 
22 //
23 //*-- Yves Schutz (CERN) & Gustavo Conesa Balbastre (INFN-LNF)
24 //////////////////////////////////////////////////////////////////////////////
25
26 #include <TChain.h>
27 #include <TFile.h> 
28 #include <TNtuple.h>
29 #include <TVector3.h> 
30
31 #include "AliAnalysisManager.h"
32 #include "AliMCEventHandler.h"
33 #include "AliMCEvent.h"
34 #include "AliAnaCaloTriggerMC.h" 
35 #include "AliESDEvent.h" 
36 #include "AliESDCaloCluster.h"
37 #include "AliLog.h"
38 #include "AliStack.h"
39 #include "TParticle.h"
40
41 //______________________________________________________________________________
42 AliAnaCaloTriggerMC::AliAnaCaloTriggerMC() : 
43   fChain(0),
44   fESD(0), 
45   fOutputContainer(0),
46   fCalorimeter("PHOS"),
47   fNtTrigger22(0), 
48   fNtTriggerNN(0)
49
50 {
51   // Default constructor.
52
53 }
54 //______________________________________________________________________________
55 AliAnaCaloTriggerMC::AliAnaCaloTriggerMC(const char *name) : 
56   AliAnalysisTask(name, "AnaCaloTriggerMC"),
57   fChain(0),
58   fESD(0), 
59   fOutputContainer(0),
60   fCalorimeter("PHOS"),
61   fNtTrigger22(0), 
62   fNtTriggerNN(0)
63
64 {
65   // Constructor.
66   // Input slot #0 works with an Ntuple
67   DefineInput(0, TChain::Class());
68   // Output slot #0 writes into a TH1 container
69   DefineOutput(0,  TObjArray::Class()) ; 
70 }
71 //____________________________________________________________________________
72 AliAnaCaloTriggerMC::AliAnaCaloTriggerMC(const AliAnaCaloTriggerMC & ct) : 
73   AliAnalysisTask(ct),fChain(ct.fChain), fESD(ct.fESD),
74   fOutputContainer(ct.fOutputContainer), fCalorimeter(ct. fCalorimeter),
75   fNtTrigger22(ct.fNtTrigger22), fNtTriggerNN(ct.fNtTriggerNN)
76 {
77
78   // cpy ctor
79   SetName (ct.GetName()) ; 
80   SetTitle(ct.GetTitle()) ;
81  
82 }
83
84 //_________________________________________________________________________
85 AliAnaCaloTriggerMC & AliAnaCaloTriggerMC::operator = (const AliAnaCaloTriggerMC & source)
86 {
87   // assignment operator
88   
89   if(&source == this) return *this;
90
91   fChain = source.fChain ; 
92   fESD = source.fESD ;
93   fOutputContainer = source.fOutputContainer ;
94   fCalorimeter = source. fCalorimeter ;
95   fNtTrigger22 = source.fNtTrigger22 ;
96   fNtTriggerNN = source.fNtTriggerNN ;
97
98   return *this;
99   
100 }
101
102 //______________________________________________________________________________
103 AliAnaCaloTriggerMC::~AliAnaCaloTriggerMC()
104 {
105   // dtor
106   //  fOutputContainer->Clear() ; 
107   //  delete fOutputContainer ;
108
109 }
110
111
112 //______________________________________________________________________________
113 void AliAnaCaloTriggerMC::ConnectInputData(const Option_t*)
114 {
115   // Initialisation of branch container and histograms 
116     
117   AliInfo(Form("*** Initialization of %s", GetName())) ; 
118   
119   // Get input data
120   fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
121   if (!fChain) {
122     AliError(Form("Input 0 for %s not found\n", GetName()));
123     return ;
124   }
125   
126   fESD = new AliESDEvent();
127   fESD->ReadFromTree(fChain);
128
129 }
130
131 //________________________________________________________________________
132
133 void AliAnaCaloTriggerMC::CreateOutputObjects()
134
135  
136   // Create the output container
137   OpenFile(0);
138  
139   // create histograms 
140   fNtTrigger22 = new TNtuple(fCalorimeter+"trigger22", "Trigger data 2x2 patch", "a22:a220:enMax:phEnMax:eta22:phi22:etaMax:phiMax:phEtaMax:phPhiMax");
141   fNtTriggerNN = new TNtuple(fCalorimeter+"triggerNN", "Trigger data NxN patch", "aNN:aNN0:enMax:phEnMax:etaNN:phiNN:etaMax:phiMax:phEtaMax:phPhiMax");
142   
143   // create output container
144   
145   fOutputContainer = new TObjArray(2) ; 
146   fOutputContainer->SetName(GetName()) ; 
147
148   fOutputContainer->AddAt(fNtTrigger22,             0) ; 
149   fOutputContainer->AddAt(fNtTriggerNN,             1) ; 
150
151 }
152
153 //______________________________________________________________________________
154 void AliAnaCaloTriggerMC::Exec(Option_t *) 
155 {
156
157   // Processing of one event
158    
159   Long64_t entry = fChain->GetReadEntry() ;
160  
161   if (!fESD) {
162     AliError("fESD is not connected to the input!") ; 
163     return ; 
164   }
165
166   if ( !((entry-1)%100) ) 
167     AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ; 
168
169   //Get MC data
170   AliStack* stack = 0x0; 
171   AliMCEventHandler*    mctruth = (AliMCEventHandler*) 
172     ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
173   
174   if(mctruth)
175     stack = mctruth->MCEvent()->Stack();
176   
177   if (!stack) {
178     AliError("Stack not found") ; 
179     return ; 
180   }
181   
182   // Get trigger information of fCalorimeter 
183   TArrayF * triggerAmplitudes = 0x0 ;
184   TArrayF * triggerPosition   = 0x0 ;
185   Int_t numberOfCaloClusters  = 0 ;
186   
187   if(fCalorimeter == "PHOS"){
188     triggerAmplitudes      = fESD->GetPHOSTriggerAmplitudes();
189     triggerPosition        = fESD->GetPHOSTriggerPosition();
190   }
191   else if(fCalorimeter == "EMCAL"){
192     triggerAmplitudes    = fESD->GetEMCALTriggerAmplitudes();
193     triggerPosition      = fESD->GetEMCALTriggerPosition();
194   }
195
196   if( triggerAmplitudes && triggerPosition ){
197     // trigger amplitudes
198     const Float_t ka22    = static_cast<Float_t>(triggerAmplitudes->At(0)) ; 
199     const Float_t ka22O   = static_cast<Float_t>(triggerAmplitudes->At(1)) ; 
200     const Float_t kaNN    = static_cast<Float_t>(triggerAmplitudes->At(2)) ; 
201     const Float_t kaNNO   = static_cast<Float_t>(triggerAmplitudes->At(3)) ; 
202     
203     // trigger position
204     const Float_t kx22  =  static_cast<Float_t>(triggerPosition->At(0)) ; 
205     const Float_t ky22  =  static_cast<Float_t>(triggerPosition->At(1)) ;
206     const Float_t kz22  =  static_cast<Float_t>(triggerPosition->At(2)) ;
207     const Float_t kxNN  =  static_cast<Float_t>(triggerPosition->At(3)) ; 
208     const Float_t kyNN  =  static_cast<Float_t>(triggerPosition->At(4)) ;
209     const Float_t kzNN  =  static_cast<Float_t>(triggerPosition->At(5)) ; 
210     
211     Float_t enMax       = 0. ;
212     Float_t phEnMax     = 0. ;
213     Float_t etaMax      = 0.5 ;
214     Float_t phiMax      = 0. ; 
215     Float_t phEtaMax    = 0.5 ;
216     Float_t phPhiMax    = 0. ; 
217     
218     TVector3 vpos22(kx22, ky22, kz22) ;
219     TVector3 vposNN(kxNN, kyNN, kzNN) ;
220     Float_t eta22 = vpos22.Eta() ; 
221     Float_t phi22 = vpos22.Phi() * TMath::RadToDeg() + 360. ; 
222     Float_t etaNN = vposNN.Eta() ; 
223     Float_t phiNN = vposNN.Phi() * TMath::RadToDeg() + 360. ; 
224     
225     Int_t      icaloCluster ; 
226     
227     // loop over the Calorimeters Clusters
228     Float_t cluEnergy = 0;
229     Int_t labelmax = -5;
230     for(icaloCluster = 0 ; icaloCluster < numberOfCaloClusters ; icaloCluster++) {
231
232       AliESDCaloCluster * cluster = fESD->GetCaloCluster(icaloCluster) ;
233
234       if (cluster && ( (fCalorimeter == "PHOS" && cluster->IsPHOS())  ||  
235                        (fCalorimeter == "EMCAL" && cluster->IsEMCAL()))) {
236
237         cluEnergy = cluster->E() ; 
238         Float_t pos[3] ;
239         TVector3 vpos ;
240         
241         cluster->GetPosition( pos ) ;
242         
243         if ( cluEnergy > enMax) { 
244           enMax = cluEnergy ; 
245           vpos.SetXYZ(pos[0], pos[1], pos[2]) ; 
246           etaMax = vpos.Eta() ; 
247           phiMax = vpos.Phi() ; 
248           labelmax = cluster->GetLabel();
249         }
250         
251         Double_t * pid = cluster->GetPid() ;
252         
253         if(pid[AliPID::kPhoton] > 0.9) {
254           if ( cluEnergy > phEnMax) { 
255             phEnMax = cluEnergy ; 
256             vpos.SetXYZ(pos[0], pos[1], pos[2]) ; 
257             phEtaMax = vpos.Eta() ; 
258             phPhiMax = vpos.Phi() ; 
259           }
260         }
261       }//if cluster
262     }//CaloCluster loop
263
264     if(labelmax < stack->GetNtrack() && labelmax >= 0 ){
265       TParticle * particle = stack->Particle(labelmax); 
266       Float_t ptgen = particle->Energy();
267       fNtTrigger22->Fill(ka22, ka22O, ptgen, enMax, phEnMax, eta22, phi22, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
268       fNtTriggerNN->Fill(kaNN, kaNNO, ptgen, enMax, phEnMax, etaNN, phiNN, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
269     }
270     else AliDebug(1, Form("Wrong label %d, ntrack %d, Emax %f ",labelmax, stack->GetNtrack(), phEnMax));
271   }//If trigger arrays filled
272     
273   PostData(0, fOutputContainer);
274   
275 }
276
277
278 //______________________________________________________________________________
279 //void AliAnaCaloTriggerMC::Terminate(Option_t *) const
280 //{
281 //  // Processing when the event loop is ended
282 //
283 //}