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