Moved GetClosestJets() method to helper class
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisHelperJetTasks.cxx
1
2 #include "TROOT.h"
3 #include "TList.h"
4 #include "AliMCEvent.h"
5 #include "AliAODJet.h"
6 #include "AliStack.h"
7 #include "AliGenEventHeader.h"
8 #include "AliGenCocktailEventHeader.h"
9 #include "AliGenPythiaEventHeader.h"
10 #include <fstream>
11 #include <iostream>
12 #include "AliAnalysisHelperJetTasks.h"
13
14
15 ClassImp(AliAnalysisHelperJetTasks)
16
17
18
19  
20 AliGenPythiaEventHeader*  AliAnalysisHelperJetTasks::GetPythiaEventHeader(AliMCEvent *mcEvent){
21   
22   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
23   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
24   if(!pythiaGenHeader){
25     // cocktail ??
26     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
27     
28     if (!genCocktailHeader) {
29       Printf("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__);
30       return 0;
31     }
32     TList* headerList = genCocktailHeader->GetHeaders();
33     for (Int_t i=0; i<headerList->GetEntries(); i++) {
34       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
35       if (pythiaGenHeader)
36         break;
37     }
38     if(!pythiaGenHeader){
39       Printf("%s %d: PythiaHeader not found!",(char*)__FILE__,__LINE__);
40       return 0;
41     }
42   }
43   return pythiaGenHeader;
44
45 }
46
47
48 void AliAnalysisHelperJetTasks::PrintStack(AliMCEvent *mcEvent,Int_t iFirst,Int_t iLast,Int_t iMaxPrint){
49
50   AliStack *stack = mcEvent->Stack();
51   if(!stack){
52     Printf("%s%d No Stack available",(char*)__FILE__,__LINE__);
53     return;
54   }
55
56   static Int_t iCount = 0;
57   if(iCount>iMaxPrint)return;
58   Int_t nStack = stack->GetNtrack();
59   if(iLast == 0)iLast = nStack;
60   else if(iLast > nStack)iLast = nStack;
61
62
63   Printf("####################################################################");
64   for(Int_t np = iFirst;np<iLast;++np){
65     TParticle *p = stack->Particle(np);
66     Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughter2 %d ",
67            np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter(0),p->GetDaughter(1));
68     Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi()); 
69     p->Print();    
70     Printf("---------------------------------------");
71   }
72   iCount++;
73 }
74
75
76
77
78 void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets,
79                                                AliAODJet *recJets,const Int_t &nRecJets,
80                                                Int_t *iGenIndex,Int_t *iRecIndex,
81                                                Int_t iDebug,Float_t maxDist){
82
83   //
84   // Relate the two input jet Arrays
85   //
86
87   //
88   // The association has to be unique
89   // So check in two directions
90   // find the closest rec to a gen
91   // and check if there is no other rec which is closer
92   // Caveat: Close low energy/split jets may disturb this correlation
93
94   // Idea: search in two directions generated e.g (a--e) and rec (1--3)
95   // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
96   // in the end we have something like this
97   //    1   2   3
98   // ------------
99   // a| 3   2   0
100   // b| 0   1   0
101   // c| 0   0   3
102   // d| 0   0   1
103   // e| 0   0   1
104   // Topology
105   //   1     2
106   //     a         b        
107   //
108   //  d      c
109   //        3     e
110   // Only entries with "3" match from both sides
111   
112   const int kMode = 3;
113
114   
115
116   if(nRecJets==0||nGenJets==0)return;
117
118   for(int i = 0;i < nGenJets;++i)iGenIndex[i] = -1;
119   for(int j = 0;j < nRecJets;++j)iRecIndex[j] = -1;
120
121   UShort_t *iFlag = new UShort_t[nGenJets*nRecJets];
122   for(int i = 0;i < nGenJets;++i){
123     for(int j = 0;j < nRecJets;++j){
124       iFlag[i*nGenJets+j] = 0;
125     }
126   }
127
128
129
130   // find the closest distance to the generated
131   for(int ig = 0;ig<nGenJets;++ig){
132     Float_t dist = maxDist;
133     if(iDebug>1)Printf("Gen (%d) p_T %3.3f eta %3.3f ph %3.3f ",ig,genJets[ig].Pt(),genJets[ig].Eta(),genJets[ig].Phi());
134     for(int ir = 0;ir<nRecJets;++ir){
135       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
136       if(iDebug>1)Printf("Rec (%d) p_T %3.3f eta %3.3f ph %3.3f ",ir,recJets[ir].Pt(),recJets[ir].Eta(),recJets[ir].Phi());
137       if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
138       if(dR<dist){
139         iRecIndex[ig] = ir;
140         dist = dR;
141       } 
142     }
143     if(iRecIndex[ig]>=0)iFlag[ig*nGenJets+iRecIndex[ig]]+=1;
144     // reset...
145     iRecIndex[ig] = -1;
146   }
147   // other way around
148   for(int ir = 0;ir<nRecJets;++ir){
149     Float_t dist = maxDist;
150     for(int ig = 0;ig<nGenJets;++ig){
151       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
152       if(dR<dist){
153         iGenIndex[ir] = ig;
154         dist = dR;
155       } 
156     }
157     if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]*nGenJets+ir]+=2;
158     // reset...
159     iGenIndex[ir] = -1;
160   }
161
162   // check for "true" correlations
163
164   if(iDebug>1)Printf(">>>>>> Matrix");
165
166   for(int ig = 0;ig<nGenJets;++ig){
167     for(int ir = 0;ir<nRecJets;++ir){
168       // Print
169       if(iDebug>1)printf("Flag[%d][%d] %d ",ig,ir,iFlag[ig*nGenJets+ir]);
170
171       if(kMode==3){
172         // we have a uniqie correlation
173         if(iFlag[ig*nGenJets+ir]==3){
174           iGenIndex[ir] = ig;
175           iRecIndex[ig] = ir;
176         }
177       }
178       else{
179         // we just take the correlation from on side
180         if((iFlag[ig*nGenJets+ir]&2)==2){
181           iGenIndex[ir] = ig;
182         }
183         if((iFlag[ig*nGenJets+ir]&1)==1){
184           iRecIndex[ig] = ir;
185         }
186       }
187     }
188     if(iDebug>1)printf("\n");
189   }
190
191   delete [] iFlag;
192
193 }
194
195
196
197
198