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