PartCorr split in 2 Base and Dep; coding violations corrected; PHOS geometry can...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaParticlePartonCorrelation.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 /* $Id:  $ */
16
17 //_________________________________________________________________________
18 // Class for the analysis of particle-parton correlations
19 // Particle (for example direct gamma) must be found in a previous analysis 
20 // -- Author: Gustavo Conesa (LNF-INFN) 
21 //////////////////////////////////////////////////////////////////////////////
22
23 // --- ROOT system ---
24 //#include "Riostream.h"
25 #include "TH2F.h"
26 #include "TParticle.h"
27
28 //---- ANALYSIS system ----
29 #include "AliAnaParticlePartonCorrelation.h" 
30 #include "AliLog.h"
31 #include "AliStack.h"  
32 #include "AliAODPWG4ParticleCorrelation.h"
33
34   ClassImp(AliAnaParticlePartonCorrelation)
35   
36
37 //____________________________________________________________________________
38   AliAnaParticlePartonCorrelation::AliAnaParticlePartonCorrelation() : 
39     AliAnaPartCorrBaseClass(),   
40     fhDeltaEtaNearParton(0), fhDeltaPhiNearParton(0), 
41     fhDeltaPtNearParton(0), fhPtRatNearParton(0),
42     fhDeltaEtaAwayParton(0), fhDeltaPhiAwayParton(0), 
43     fhDeltaPtAwayParton(0), fhPtRatAwayParton(0)
44 {
45   //Default Ctor
46
47   //Initialize parameters
48   InitParameters();
49 }
50
51 //____________________________________________________________________________
52 AliAnaParticlePartonCorrelation::AliAnaParticlePartonCorrelation(const AliAnaParticlePartonCorrelation & g) :   
53   AliAnaPartCorrBaseClass(g),   
54   fhDeltaEtaNearParton(g.fhDeltaEtaNearParton), fhDeltaPhiNearParton(g.fhDeltaPhiNearParton), 
55   fhDeltaPtNearParton(g.fhDeltaPtNearParton), fhPtRatNearParton(g.fhPtRatNearParton),
56   fhDeltaEtaAwayParton(g.fhDeltaEtaAwayParton), fhDeltaPhiAwayParton(g.fhDeltaPhiAwayParton), 
57   fhDeltaPtAwayParton(g.fhDeltaPtAwayParton), fhPtRatAwayParton(g.fhPtRatAwayParton)
58 {
59   // cpy ctor
60
61 }
62
63 //_________________________________________________________________________
64 AliAnaParticlePartonCorrelation & AliAnaParticlePartonCorrelation::operator = (const AliAnaParticlePartonCorrelation & source)
65 {
66   // assignment operator
67
68   if(this == &source)return *this;
69   ((AliAnaPartCorrBaseClass *)this)->operator=(source);
70
71   fhDeltaEtaAwayParton = source.fhDeltaEtaAwayParton;
72   fhDeltaPhiAwayParton = source.fhDeltaPhiAwayParton;
73   fhDeltaPtAwayParton = source.fhDeltaPtAwayParton;
74   fhPtRatAwayParton = source.fhPtRatAwayParton;
75   fhDeltaEtaNearParton = source.fhDeltaEtaNearParton;
76   fhDeltaPhiNearParton = source.fhDeltaPhiNearParton;
77   fhDeltaPtNearParton = source.fhDeltaPtNearParton;
78   fhPtRatNearParton = source.fhPtRatNearParton;
79
80   return *this;
81
82 }
83
84
85 //________________________________________________________________________
86 TList *  AliAnaParticlePartonCorrelation::GetCreateOutputObjects()
87 {  
88   // Create histograms to be saved in output file 
89
90   AliDebug(1,"Init parton histograms");
91
92   TList * outputContainer = new TList() ; 
93   outputContainer->SetName("ParticlePartonHistos") ; 
94   
95   fhDeltaPhiNearParton  = new TH2F
96     ("DeltaPhiNearParton","#phi_{particle} - #phi_{parton} vs p_{T particle}",
97      200,0,120,200,0,6.4); 
98   fhDeltaPhiNearParton->SetYTitle("#Delta #phi");
99   fhDeltaPhiNearParton->SetXTitle("p_{T particle} (GeV/c)");
100   outputContainer->Add(fhDeltaPhiNearParton);
101   
102   fhDeltaEtaNearParton  = new TH2F
103     ("DeltaEtaNearParton","#eta_{particle} - #eta_{parton} vs p_{T particle}",
104      200,0,120,200,-2,2); 
105   fhDeltaEtaNearParton->SetYTitle("#Delta #eta");
106   fhDeltaEtaNearParton->SetXTitle("p_{T particle} (GeV/c)");
107   outputContainer->Add(fhDeltaEtaNearParton);
108   
109   fhDeltaPtNearParton  = new TH2F
110     ("DeltaPtNearParton","#p_{T particle} - #p_{T parton} vs p_{T particle}",
111      200,0,120,100,-10,10); 
112   fhDeltaPtNearParton->SetYTitle("#Delta #p_{T}");
113   fhDeltaPtNearParton->SetXTitle("p_{T particle} (GeV/c)"); 
114   outputContainer->Add(fhDeltaPtNearParton);
115   
116   fhPtRatNearParton  = new TH2F
117     ("PtRatNearParton","#p_{T parton} / #p_{T particle} vs p_{T particle}",
118      200,0,120,200,0,5); 
119   fhPtRatNearParton->SetYTitle("ratio");
120   fhPtRatNearParton->SetXTitle("p_{T particle} (GeV/c)");
121   outputContainer->Add(fhPtRatNearParton);
122   
123   fhDeltaPhiAwayParton  = new TH2F
124     ("DeltaPhiAwayParton","#phi_{particle} - #phi_{parton} vs p_{T particle}",
125      200,0,120,200,0,6.4); 
126   fhDeltaPhiAwayParton->SetYTitle("#Delta #phi");
127   fhDeltaPhiAwayParton->SetXTitle("p_{T particle} (GeV/c)");
128   outputContainer->Add(fhDeltaPhiAwayParton);
129   
130   fhDeltaEtaAwayParton  = new TH2F
131     ("DeltaEtaAwayParton","#eta_{particle} - #eta_{parton} vs p_{T particle}",
132      200,0,120,200,-2,2); 
133   fhDeltaEtaAwayParton->SetYTitle("#Delta #eta");
134   fhDeltaEtaAwayParton->SetXTitle("p_{T particle} (GeV/c)");
135   outputContainer->Add(fhDeltaEtaAwayParton);
136   
137   fhDeltaPtAwayParton  = new TH2F
138     ("DeltaPtAwayParton","#p_{T particle} - #p_{T parton} vs p_{T particle}",
139      200,0,120,100,-10,10); 
140   fhDeltaPtAwayParton->SetYTitle("#Delta #p_{T}");
141   fhDeltaPtAwayParton->SetXTitle("p_{T particle} (GeV/c)"); 
142   outputContainer->Add(fhDeltaPtAwayParton);
143   
144   fhPtRatAwayParton  = new TH2F
145     ("PtRatAwayParton","#p_{T parton} / #p_{T particle} vs p_{T particle}",
146      200,0,120,200,0,5); 
147   fhPtRatAwayParton->SetYTitle("ratio");
148   fhPtRatAwayParton->SetXTitle("p_{T particle} (GeV/c)");
149   outputContainer->Add(fhPtRatAwayParton);
150   
151   return outputContainer;
152 }
153
154 //____________________________________________________________________________
155 void AliAnaParticlePartonCorrelation::InitParameters()
156 {
157   
158   //Initialize the parameters of the analysis.
159   SetInputAODName("photons");
160   
161 }
162
163 //__________________________________________________________________
164 void AliAnaParticlePartonCorrelation::Print(const Option_t * opt) const
165 {
166   
167   //Print some relevant parameters set for the analysis
168   if(! opt)
169     return;
170   
171
172
173 //__________________________________________________________________
174 void  AliAnaParticlePartonCorrelation::MakeAnalysisFillAOD()  
175 {
176   //Particle-Parton Correlation Analysis, create AODs
177   //Add partons to the reference list of the trigger particle
178   //Partons are considered those in the first eight possitions in the stack
179   //being 0, and 1 the 2 protons, and 6 and 7 the outgoing final partons.
180   if(!GetInputAODBranch())
181         AliFatal(Form("ParticlePartonCorrelation::FillAOD: No input particles in AOD with name branch < %s > \n",GetInputAODName().Data()));    
182   
183   if(GetDebug() > 1){
184     printf("Begin parton correlation analysis, fill AODs \n");
185     printf("In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
186   }
187   
188   //Loop on stored AOD particles
189   Int_t naod = GetInputAODBranch()->GetEntriesFast();
190   for(Int_t iaod = 0; iaod < naod ; iaod++){
191     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
192     
193     AliStack * stack =  GetMCStack() ;
194     if(!stack) AliFatal("No Stack available, STOP");
195     
196     if(stack->GetNtrack() < 8) {
197       printf("*** small number of particles, not a PYTHIA simulation? ***:  n tracks %d \n", stack->GetNprimary());
198       continue ;
199     }
200     
201     //Fill AOD reference only with partons
202     TParticle * parton = new TParticle ;
203     
204     for(Int_t ipr = 0;ipr < 8; ipr ++ ){
205       parton = stack->Particle(ipr) ;
206       particle->AddTrack(parton);
207       //parton->Print();
208     }
209
210   }//Aod branch loop
211
212  if(GetDebug() > 1) printf("End parton correlation analysis, fill AODs \n");
213 }
214  
215 //__________________________________________________________________
216 void  AliAnaParticlePartonCorrelation::MakeAnalysisFillHistograms()  
217 {
218   //Particle-Parton Correlation Analysis, fill histograms
219   if(!GetInputAODBranch())
220         AliFatal(Form("ParticlePartonCorrelation::FillHistos: No input particles in AOD with name branch < %s > \n",GetInputAODName().Data())); 
221   
222   if(GetDebug() > 1){
223     printf("Begin parton correlation analysis, fill histograms \n");
224     printf("In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
225   }
226
227   AliStack * stack =  GetMCStack() ;
228   if(!stack) AliFatal("No Stack available, STOP");
229
230   //Loop on stored AOD particles
231   Int_t naod = GetInputAODBranch()->GetEntriesFast();
232   TParticle *  mom =new TParticle ;
233
234   for(Int_t iaod = 0; iaod < naod ; iaod++){
235     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
236     
237     Float_t ptTrigg  = particle->Pt();
238     Float_t phiTrigg = particle->Phi();
239     Float_t etaTrigg = particle->Eta(); 
240     Int_t imom = particle->GetLabel();
241     Int_t iparent  = 2000;
242     Int_t iawayparent = -1;
243
244     if(!(particle->GetRefTracks()) || (particle->GetRefTracks())->GetEntriesFast() < 7) AliFatal("Reference list with partons not filled, STOP analysis");
245
246     //Check and get indeces of mother and parton    
247     if(imom < 8 ) iparent = imom ;   //mother is already a parton
248     else if (imom <  stack->GetNtrack()) {
249       mom =  stack->Particle(imom);
250       iparent=mom->GetFirstMother();
251       //cout<<" iparent "<<iparent<<endl;
252       while(iparent > 7 ){
253         mom = stack->Particle(iparent);
254         imom = iparent ; //Mother label is of the inmediate parton daughter
255         iparent = mom->GetFirstMother();
256         //cout<<" while iparent "<<iparent<<endl;
257       }   
258     }
259     
260     if(GetDebug() > 1) printf("N reference partons %d; labels:  mother %d, parent %d \n", (particle->GetRefTracks())->GetEntriesFast(), imom, iparent);
261
262     
263     if(iparent < 0 || iparent > 8) { 
264       if(GetDebug() > 0 ) printf("Failed to find appropriate parton, index %d", iparent);
265       continue ;
266     }
267
268     //Near parton is the parton that fragmented and created the mother    
269     TParticle * nearParton = (TParticle*) (particle->GetRefTracks())->At(iparent);
270     Float_t  ptNearParton    = nearParton->Pt();
271     Float_t  phiNearParton   = nearParton->Phi() ;
272     Float_t  etaNearParton   = nearParton->Eta() ;
273     
274     fhDeltaEtaNearParton->Fill(ptTrigg,etaTrigg-etaNearParton);
275     fhDeltaPhiNearParton->Fill(ptTrigg,phiTrigg-phiNearParton);
276     fhDeltaPtNearParton->Fill(ptTrigg,ptTrigg-ptNearParton);
277     fhPtRatNearParton->Fill(ptTrigg,ptNearParton/ptTrigg);
278     
279     if(iparent == 7) iawayparent =6;
280     else if(iparent == 6) iawayparent =7;
281     else{
282       printf("Parent parton is not final state, skip \n");
283       continue;
284     }
285
286     //Away parton is the other final parton.
287     TParticle * awayParton = (TParticle*) (particle->GetRefTracks())->At(iawayparent);
288     Float_t  ptAwayParton    = awayParton->Pt();
289     Float_t  phiAwayParton   = awayParton->Phi() ;
290     Float_t  etaAwayParton   = awayParton->Eta() ;
291     fhDeltaEtaAwayParton->Fill(ptTrigg,etaTrigg-etaAwayParton);
292     fhDeltaPhiAwayParton->Fill(ptTrigg,phiTrigg-phiAwayParton);
293     fhDeltaPtAwayParton->Fill(ptTrigg,ptTrigg-ptAwayParton);
294     fhPtRatAwayParton->Fill(ptTrigg,ptAwayParton/ptTrigg);
295     
296  
297   }
298
299   if(GetDebug() > 1) printf("End parton correlation analysis, fill histograms \n");
300   
301