1)Terminate() method implemented in the frame. Simple examples on what to do with...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliMCAnalysisUtils.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: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
16
17 //_________________________________________________________________________
18 // Class for analysis utils for MC data
19 // stored in stack or event header.
20 // Contains:
21 //  - method to check the origin of a given track/cluster
22 //  - method to obtain the generated jets
23 //                
24 //*-- Author: Gustavo Conesa (LNF-INFN) 
25 //////////////////////////////////////////////////////////////////////////////
26
27
28 // --- ROOT system ---
29 #include <TMath.h>
30 #include <TString.h>
31 #include <TList.h>
32
33 //---- ANALYSIS system ----
34 #include "AliLog.h"
35 #include "AliMCAnalysisUtils.h"
36 #include "AliStack.h"
37 #include "TParticle.h"
38 #include "AliGenPythiaEventHeader.h"
39
40 ClassImp(AliMCAnalysisUtils)
41
42
43 //________________________________________________
44 AliMCAnalysisUtils::AliMCAnalysisUtils() : 
45 TObject(), fCurrentEvent(-1), fDebug(-1), 
46 fJetsList(new TList), fMCGenerator("PYTHIA")
47 {
48         //Ctor
49 }
50
51 //____________________________________________________________________________
52 AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :   
53 TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
54 fJetsList(mcutils.fJetsList), fMCGenerator(mcutils.fMCGenerator)
55 {
56         // cpy ctor
57         
58 }
59
60 //_________________________________________________________________________
61 AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
62 {
63         // assignment operator
64         
65         if(&mcutils == this) return *this;
66         fCurrentEvent = mcutils.fCurrentEvent ;
67         fDebug        = mcutils.fDebug;
68         fJetsList     = mcutils.fJetsList;
69         fMCGenerator  = mcutils.fMCGenerator;
70
71         return *this;
72         
73 }
74
75 //____________________________________________________________________________
76 AliMCAnalysisUtils::~AliMCAnalysisUtils() 
77 {
78         // Remove all pointers.
79         
80         if (fJetsList) {
81                 fJetsList->Clear();
82                 delete fJetsList ;
83         }   
84
85 }
86
87 //_________________________________________________________________________
88 Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliStack * stack) const {
89   //Play with the MC stack if available
90   //Check origin of the candidates, good for PYTHIA
91   
92   if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
93   //  printf("label %d, ntrack %d, nprim %d\n",label, stack->GetNtrack(), stack->GetNprimary());
94 //   for(Int_t i = 0; i< stack->GetNprimary(); i++){
95 //      TParticle *particle =   stack->Particle(i);
96 //                      //particle->Print();
97 //   }
98   if(label >= 0 && label <  stack->GetNtrack()){
99     //Mother
100     TParticle * mom = stack->Particle(label);
101     Int_t mPdg = TMath::Abs(mom->GetPdgCode());
102     Int_t mStatus =  mom->GetStatusCode() ;
103     Int_t iParent =  mom->GetFirstMother() ;
104     if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOrigin: Mother is parton %d\n",iParent);
105     
106     //GrandParent
107     TParticle * parent = new TParticle ;
108     Int_t pPdg = -1;
109     Int_t pStatus =-1;
110     if(iParent > 0){
111       parent = stack->Particle(iParent);
112       pPdg = TMath::Abs(parent->GetPdgCode());
113       pStatus = parent->GetStatusCode();  
114     }
115     else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOrigin: Parent with label %d\n",iParent);
116     
117     //return tag
118     if(mPdg == 22){
119       if(mStatus == 1){
120         if(fMCGenerator == "PYTHIA"){
121           if(iParent < 8 && iParent > 5) {//outgoing partons
122             if(pPdg == 22) return kMCPrompt;
123             else  return kMCFragmentation;
124           }//Outgoing partons
125           else if(pStatus == 11){//Decay
126             if(pPdg == 111) return kMCPi0Decay ;
127             else if (pPdg == 321)  return kMCEtaDecay ;
128             else  return kMCOtherDecay ;
129           }//Decay
130           else return kMCISR; //Initial state radiation
131         }//PYTHIA
132
133         else if(fMCGenerator == "HERWIG"){        
134           if(pStatus < 197){//Not decay
135             while(1){
136               if(parent->GetFirstMother()<=5) break;
137               iParent = parent->GetFirstMother();
138               parent=stack->Particle(iParent);
139               pStatus= parent->GetStatusCode();
140               pPdg = parent->GetPdgCode();
141             }//Look for the parton
142             
143             if(iParent < 8 && iParent > 5) {
144               if(pPdg == 22) return kMCPrompt;
145               else  return kMCFragmentation;
146             }
147             return kMCISR;//Initial state radiation
148           }//Not decay
149           else{//Decay
150             if(pPdg == 111) return kMCPi0Decay ;
151             else if (pPdg == 321)  return kMCEtaDecay ;
152             else  return kMCOtherDecay ;
153           }//Decay
154         }//HERWIG
155         else return  kMCUnknown;
156       }//Status 1 : Pythia generated
157       else if(mStatus == 0){
158         if(pPdg ==22 || pPdg ==11|| pPdg == 2112 ||  pPdg == 211 ||  
159            pPdg == 321 ||  pPdg == 2212  ||  pPdg == 130  ||  pPdg == 13 ) 
160           return kMCConversion ;
161         if(pPdg == 111) return kMCPi0Decay ;
162         else if (pPdg == 221)  return kMCEtaDecay ;
163         else  return kMCOtherDecay ;
164       }//status 0 : geant generated
165     }//Mother Photon
166     else if(mPdg == 111)  return kMCPi0 ;
167     else if(mPdg == 221)  return kMCEta ;
168     else if(mPdg ==11){
169 //       if(pPdg !=22 &&mStatus == 0) {
170 //      printf("Origin electron, pT %f, status %d, parent %s, pstatus %d, vx %f, vy %f, vz %f\n",
171 //             mom->Pt(),mStatus, parent->GetName(),pStatus,mom->Vx(),mom->Vy(), mom->Vz());
172
173 //       }
174
175       if(mStatus == 0) {
176         if(pPdg ==22) return kMCConversion ;
177         else if (pStatus == 0) return kMCConversion;
178         else return kMCElectron ;
179       }
180       else return kMCElectron ;
181     }
182     else return kMCUnknown;
183   }//Good label value
184   else{
185     if(label < 0 ) printf("AliMCAnalysisUtils::CheckOrigin: *** bad label or no stack ***:  label %d \n", label);
186     if(label >=  stack->GetNtrack()) printf("AliMCAnalysisUtils::CheckOrigin: *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
187     return kMCUnknown;
188   }//Bad label
189         
190   return kMCUnknown;
191   
192 }
193
194 //_________________________________________________________________________
195 TList * AliMCAnalysisUtils::GetJets(Int_t iEvent, AliStack * stack, AliGenEventHeader * geh) {
196  //Return list of jets (TParticles) and index of most likely parton that originated it.
197         
198         if(fCurrentEvent!=iEvent){
199                 fCurrentEvent = iEvent;
200                 fJetsList = new TList;
201                 Int_t nTriggerJets = 0;
202                 Float_t tmpjet[]={0,0,0,0};
203                 
204                 //printf("Event %d %d\n",fCurrentEvent,iEvent);
205                 //Get outgoing partons
206                 if(stack->GetNtrack() < 8) return fJetsList;
207                 TParticle * parton1 =  stack->Particle(6);
208                 TParticle * parton2 =  stack->Particle(7);
209                 if(fDebug > 2){
210                         printf("parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
211                                 parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
212                         printf("parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
213                                 parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
214                 }
215 //              //Trace the jet from the mother parton
216 //              Float_t pt  = 0;
217 //              Float_t pt1 = 0;
218 //              Float_t pt2 = 0;
219 //              Float_t e   = 0;
220 //              Float_t e1  = 0;
221 //              Float_t e2  = 0;
222 //              TParticle * tmptmp = new TParticle;
223 //              for(Int_t i = 0; i< stack->GetNprimary(); i++){
224 //                      tmptmp = stack->Particle(i);
225                 
226 //                      if(tmptmp->GetStatusCode() == 1){
227 //                              pt = tmptmp->Pt();
228 //                              e =  tmptmp->Energy();                  
229 //                              Int_t imom = tmptmp->GetFirstMother();
230 //                              Int_t imom1 = 0;
231 //                              //printf("1st imom %d\n",imom);
232 //                              while(imom > 5){
233 //                                      imom1=imom;
234 //                                      tmptmp = stack->Particle(imom);
235 //                                      imom = tmptmp->GetFirstMother();
236 //                                      //printf("imom %d       \n",imom);
237 //                              }
238 //                              //printf("Last imom %d %d\n",imom1, imom);
239 //                              if(imom1 == 6) {
240 //                                      pt1+=pt;
241 //                                      e1+=e;                          
242 //                              }
243 //                              else if (imom1 == 7){
244 //                                      pt2+=pt;
245 //                                      e2+=e;                                  }
246 //                      }// status 1
247                                 
248 //              }// for
249                 
250 //              printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
251                 
252                 //Get the jet, different way for different generator
253                 //PYTHIA
254                 if(fMCGenerator == "PYTHIA"){
255                         TParticle * jet =  new TParticle;
256                         AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
257                         nTriggerJets =  pygeh->NTriggerJets();
258                         //if(fDebug > 1)
259                                 printf("PythiaEventHeader: Njets: %d\n",nTriggerJets);
260                         Int_t iparton = -1;
261                         for(Int_t i = 0; i< nTriggerJets; i++){
262                                 iparton=-1;
263                                 pygeh->TriggerJet(i, tmpjet);
264                                 jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
265                                 //Assign an outgoing parton as mother
266                                 Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());               
267                                 Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
268                                 if(phidiff1 > phidiff2) jet->SetFirstMother(7);
269                                 else  jet->SetFirstMother(6);
270                                 //jet->Print();
271                                 if(fDebug > 1)
272                                   printf("PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
273                                          i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
274                                 fJetsList->Add(jet);                    
275                         }
276                 }//Pythia triggered jets
277                 //HERWIG
278                 else if (fMCGenerator=="HERWIG"){
279                         Int_t pdg = -1;         
280                         //Check parton 1
281                         TParticle * tmp = parton1;
282                         if(parton1->GetPdgCode()!=22){
283                                 while(pdg != 94){
284                                         if(tmp->GetFirstDaughter()==-1) return fJetsList;
285                                         tmp = stack->Particle(tmp->GetFirstDaughter());
286                                         pdg = tmp->GetPdgCode();
287                                 }//while
288                                 
289                                 //Add found jet to list
290                                 TParticle *jet1 = new TParticle(*tmp);
291                                 jet1->SetFirstMother(6);
292                                 fJetsList->Add(jet1);
293                                 //printf("jet 1:  first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
294                                 //tmp = stack->Particle(tmp->GetFirstDaughter());
295                                 //tmp->Print();
296                                 //jet1->Print();
297                                 if(fDebug > 1)                  
298                                         printf("HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
299                                                 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
300                         }//not photon
301                                 
302                         //Check parton 2
303                         pdg = -1;
304                         tmp = parton2;
305                         Int_t i = -1;
306                         if(parton2->GetPdgCode()!=22){
307                                 while(pdg != 94){
308                                         if(tmp->GetFirstDaughter()==-1) return fJetsList;
309                                         i = tmp->GetFirstDaughter();
310                                         tmp = stack->Particle(tmp->GetFirstDaughter());
311                                         pdg = tmp->GetPdgCode();
312                                 }//while
313                                 //Add found jet to list
314                                 TParticle *jet2 = new TParticle(*tmp);
315                                 jet2->SetFirstMother(7);
316                                 fJetsList->Add(jet2);
317                                 //jet2->Print();
318                                 if(fDebug > 1)
319                                         printf("HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
320                                         jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
321                                         Int_t first =  tmp->GetFirstDaughter();
322                                         Int_t last  =  tmp->GetLastDaughter();
323                                         printf("jet 2:  first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
324                                 //      for(Int_t d = first ; d < last+1; d++){
325 //                                              tmp = stack->Particle(d);
326 //                                              if(i == tmp->GetFirstMother())
327 //                                                      printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
328 //                                                      d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());                    
329 //                         }
330                                                    //tmp->Print();
331                            }//not photon
332                 }//Herwig generated jets
333         }
334         
335         return fJetsList;
336 }
337
338 //________________________________________________________________
339 void AliMCAnalysisUtils::Print(const Option_t * opt) const
340 {
341         
342         //Print some relevant parameters set for the analysis
343         if(! opt)
344                 return;
345         
346         printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
347         
348         printf("Debug level    = %d\n",fDebug);
349         printf("MC Generator   = %s\n",fMCGenerator.Data());
350
351         printf(" \n");
352         
353
354
355