5c3c2d33011e7c4ded3e44cd36e73d1577b810c6
[u/mrichter/AliRoot.git] / PWG4 / 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 /* History of cvs commits:
18  *
19  * $Log$
20  *
21  *
22  */
23 //_________________________________________________________________________
24 // Class for the analysis of particle-parton correlations
25 // Particle (for example direct gamma) must be found in a previous analysis 
26 // -- Author: Gustavo Conesa (LNF-INFN) 
27 //////////////////////////////////////////////////////////////////////////////
28
29 // --- ROOT system ---
30 #include "Riostream.h"
31 #include "TH2F.h"
32 #include "TParticle.h"
33
34 //---- ANALYSIS system ----
35 #include "AliAnaParticlePartonCorrelation.h" 
36 #include "AliLog.h"
37 #include "AliStack.h"  
38
39   ClassImp(AliAnaParticlePartonCorrelation)
40   
41
42 //____________________________________________________________________________
43   AliAnaParticlePartonCorrelation::AliAnaParticlePartonCorrelation() : 
44     AliAnaBaseClass(),   
45     fhDeltaEtaNearParton(0), fhDeltaPhiNearParton(0), 
46     fhDeltaPtNearParton(0), fhPtRatNearParton(0),
47     fhDeltaEtaAwayParton(0), fhDeltaPhiAwayParton(0), 
48     fhDeltaPtAwayParton(0), fhPtRatAwayParton(0)
49 {
50   //Default Ctor
51
52   //Initialize parameters
53   InitParameters();
54 }
55
56 //____________________________________________________________________________
57 AliAnaParticlePartonCorrelation::AliAnaParticlePartonCorrelation(const AliAnaParticlePartonCorrelation & g) :   
58   AliAnaBaseClass(g),   
59   fhDeltaEtaNearParton(g.fhDeltaEtaNearParton), fhDeltaPhiNearParton(g.fhDeltaPhiNearParton), 
60   fhDeltaPtNearParton(g.fhDeltaPtNearParton), fhPtRatNearParton(g.fhPtRatNearParton),
61   fhDeltaEtaAwayParton(g.fhDeltaEtaAwayParton), fhDeltaPhiAwayParton(g.fhDeltaPhiAwayParton), 
62   fhDeltaPtAwayParton(g.fhDeltaPtAwayParton), fhPtRatAwayParton(g.fhPtRatAwayParton)
63 {
64   // cpy ctor
65
66 }
67
68 //_________________________________________________________________________
69 AliAnaParticlePartonCorrelation & AliAnaParticlePartonCorrelation::operator = (const AliAnaParticlePartonCorrelation & source)
70 {
71   // assignment operator
72
73   if(this == &source)return *this;
74   ((AliAnaBaseClass *)this)->operator=(source);
75
76   fhDeltaEtaAwayParton = source.fhDeltaEtaAwayParton;
77   fhDeltaPhiAwayParton = source.fhDeltaPhiAwayParton;
78   fhDeltaPtAwayParton = source.fhDeltaPtAwayParton;
79   fhPtRatAwayParton = source.fhPtRatAwayParton;
80   fhDeltaEtaNearParton = source.fhDeltaEtaNearParton;
81   fhDeltaPhiNearParton = source.fhDeltaPhiNearParton;
82   fhDeltaPtNearParton = source.fhDeltaPtNearParton;
83   fhPtRatNearParton = source.fhPtRatNearParton;
84
85   return *this;
86
87 }
88
89
90 //________________________________________________________________________
91 TList *  AliAnaParticlePartonCorrelation::GetCreateOutputObjects()
92 {  
93   // Create histograms to be saved in output file 
94
95   AliDebug(1,"Init parton histograms");
96
97   TList * outputContainer = new TList() ; 
98   outputContainer->SetName("ParticlePartonHistos") ; 
99   
100   fhDeltaPhiNearParton  = new TH2F
101     ("DeltaPhiNearParton","#phi_{particle} - #phi_{parton} vs p_{T particle}",
102      200,0,120,200,0,6.4); 
103   fhDeltaPhiNearParton->SetYTitle("#Delta #phi");
104   fhDeltaPhiNearParton->SetXTitle("p_{T particle} (GeV/c)");
105   outputContainer->Add(fhDeltaPhiNearParton);
106   
107   fhDeltaEtaNearParton  = new TH2F
108     ("DeltaEtaNearParton","#eta_{particle} - #eta_{parton} vs p_{T particle}",
109      200,0,120,200,-2,2); 
110   fhDeltaEtaNearParton->SetYTitle("#Delta #eta");
111   fhDeltaEtaNearParton->SetXTitle("p_{T particle} (GeV/c)");
112   outputContainer->Add(fhDeltaEtaNearParton);
113   
114   fhDeltaPtNearParton  = new TH2F
115     ("DeltaPtNearParton","#p_{T particle} - #p_{T parton} vs p_{T particle}",
116      200,0,120,100,-10,10); 
117   fhDeltaPtNearParton->SetYTitle("#Delta #p_{T}");
118   fhDeltaPtNearParton->SetXTitle("p_{T particle} (GeV/c)"); 
119   outputContainer->Add(fhDeltaPtNearParton);
120   
121   fhPtRatNearParton  = new TH2F
122     ("PtRatNearParton","#p_{T parton} / #p_{T particle} vs p_{T particle}",
123      200,0,120,200,0,5); 
124   fhPtRatNearParton->SetYTitle("ratio");
125   fhPtRatNearParton->SetXTitle("p_{T particle} (GeV/c)");
126   outputContainer->Add(fhPtRatNearParton);
127   
128   fhDeltaPhiAwayParton  = new TH2F
129     ("DeltaPhiAwayParton","#phi_{particle} - #phi_{parton} vs p_{T particle}",
130      200,0,120,200,0,6.4); 
131   fhDeltaPhiAwayParton->SetYTitle("#Delta #phi");
132   fhDeltaPhiAwayParton->SetXTitle("p_{T particle} (GeV/c)");
133   outputContainer->Add(fhDeltaPhiAwayParton);
134   
135   fhDeltaEtaAwayParton  = new TH2F
136     ("DeltaEtaAwayParton","#eta_{particle} - #eta_{parton} vs p_{T particle}",
137      200,0,120,200,-2,2); 
138   fhDeltaEtaAwayParton->SetYTitle("#Delta #eta");
139   fhDeltaEtaAwayParton->SetXTitle("p_{T particle} (GeV/c)");
140   outputContainer->Add(fhDeltaEtaAwayParton);
141   
142   fhDeltaPtAwayParton  = new TH2F
143     ("DeltaPtAwayParton","#p_{T particle} - #p_{T parton} vs p_{T particle}",
144      200,0,120,100,-10,10); 
145   fhDeltaPtAwayParton->SetYTitle("#Delta #p_{T}");
146   fhDeltaPtAwayParton->SetXTitle("p_{T particle} (GeV/c)"); 
147   outputContainer->Add(fhDeltaPtAwayParton);
148   
149   fhPtRatAwayParton  = new TH2F
150     ("PtRatAwayParton","#p_{T parton} / #p_{T particle} vs p_{T particle}",
151      200,0,120,200,0,5); 
152   fhPtRatAwayParton->SetYTitle("ratio");
153   fhPtRatAwayParton->SetXTitle("p_{T particle} (GeV/c)");
154   outputContainer->Add(fhPtRatAwayParton);
155   
156   return outputContainer;
157 }
158
159 //____________________________________________________________________________
160 void AliAnaParticlePartonCorrelation::InitParameters()
161 {
162   
163   //Initialize the parameters of the analysis.
164   ;
165   
166 }
167
168 //__________________________________________________________________
169 void AliAnaParticlePartonCorrelation::Print(const Option_t * opt) const
170 {
171   
172   //Print some relevant parameters set for the analysis
173   if(! opt)
174     return;
175   
176
177
178 //__________________________________________________________________
179 void  AliAnaParticlePartonCorrelation::MakeAnalysisFillAOD()  
180 {
181   //Particle-Parton Correlation Analysis, create AODs
182   //Add partons to the reference list of the trigger particle
183   //Partons are considered those in the first eight possitions in the stack
184   //being 0, and 1 the 2 protons, and 6 and 7 the outgoing final partons.
185
186   if(GetDebug() > 1){
187     printf("Begin parton correlation analysis, fill AODs \n");
188     printf("In particle branch aod entries %d\n", GetAODBranch()->GetEntries());
189   }
190   
191   //Loop on stored AOD particles
192   Int_t naod = GetAODBranch()->GetEntriesFast();
193   for(Int_t iaod = 0; iaod < naod ; iaod++){
194     AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
195     
196     AliStack * stack =  GetMCStack() ;
197     if(!stack) AliFatal("No Stack available, STOP");
198     
199     if(stack->GetNtrack() < 8) {
200       printf("*** small number of particles, not a PYTHIA simulation? ***:  n tracks %d \n", stack->GetNprimary());
201       continue ;
202     }
203     
204     //Fill AOD reference only with partons
205     TParticle * parton = new TParticle ;
206     
207     for(Int_t ipr = 0;ipr < 8; ipr ++ ){
208       parton = stack->Particle(ipr) ;
209       particle->AddTrack(parton);
210       //parton->Print();
211     }
212
213   }//Aod branch loop
214
215  if(GetDebug() > 1) printf("End parton correlation analysis, fill AODs \n");
216 }
217  
218 //__________________________________________________________________
219 void  AliAnaParticlePartonCorrelation::MakeAnalysisFillHistograms()  
220 {
221   //Particle-Parton Correlation Analysis, fill histograms
222   if(GetDebug() > 1){
223     printf("Begin parton correlation analysis, fill histograms \n");
224     printf("In particle branch aod entries %d\n", GetAODBranch()->GetEntries());
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 = GetAODBranch()->GetEntriesFast();
232   TParticle *  mom =new TParticle ;
233
234   for(Int_t iaod = 0; iaod < naod ; iaod++){
235     AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->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())->GetEntries() < 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())->GetEntries(), 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