Step size region by region.
[u/mrichter/AliRoot.git] / JETAN / AliJetEventParticles.cxx
1 // $Id$
2
3 //__________________________________________________________
4 ///////////////////////////////////////////////////////////////////
5 //
6 // class AliJetEventParticles
7 //
8 // loizides@ikf.uni-frankfurt.de
9 //
10 ///////////////////////////////////////////////////////////////////
11
12 #include <Riostream.h>
13 #include <TClonesArray.h>
14 #include "AliJetParticle.h"
15 #include "AliJetEventParticles.h"
16
17 ClassImp(AliJetEventParticles)
18
19 AliJetEventParticles::AliJetEventParticles(Int_t size) :
20   fNParticles(0),
21   fParticles(new TClonesArray("AliJetParticle",size)),
22   fVertexX(0.),
23   fVertexY(0.),
24   fVertexZ(0.),
25   fTrials(0),
26   fNJets(0),
27   fNUQJets(0),
28   fXJet(-1),
29   fYJet(-1),
30   fImpact(0.),
31   fNHardScatters(0),
32   fNwNwColl(0)
33 {
34   // Default Constructor
35   for (Int_t i = 0; i < 4; i++) fZquench[i] = 0.;
36   for (Int_t i = 0; i < 10; i++) 
37     for (Int_t j = 0; j < 4; j++) {
38       fJets[j][i]=0;    // Trigger jets
39       fUQJets[j][i]=0;  // Unquenched trigger jets
40     }
41   for (Int_t i = 0; i < 5; i++){
42     fHard[i][0]=0;
43     fHard[i][1]=0;
44   }
45 }
46
47 AliJetEventParticles::AliJetEventParticles(const AliJetEventParticles& source) :
48   TObject(source), 
49   fNParticles(source.fNParticles),
50   fParticles(new TClonesArray("AliJetParticle",source.fNParticles)),
51   fVertexX(source.GetVertexX()),
52   fVertexY(source.GetVertexY()),
53   fVertexZ(source.GetVertexZ()),
54   fTrials(source.Trials()),
55   fNJets(source.NTriggerJets()),
56   fNUQJets(source.NUQTriggerJets()),
57   fXJet(source.GetXJet()),
58   fYJet(source.GetXJet()),
59   fImpact(source.GetImpact()),
60   fNHardScatters(source.GetNhard()),
61   fNwNwColl(source.GetNpart())
62 {
63   //copy constructor
64   for(Int_t i =0; i<fNParticles; i++)
65     {
66       const AliJetParticle *kjp=(const AliJetParticle *)source.fParticles->At(i);
67       new((*fParticles)[i]) AliJetParticle(*(kjp));
68     }
69   for (Int_t i = 0; i < 4; i++) fZquench[i] = 0.;
70   for (Int_t i = 0; i < 10; i++) 
71     for (Int_t j = 0; j < 4; j++) {
72       fJets[j][i]=0;    // Trigger jets
73       fUQJets[j][i]=0;  // Unquenched trigger jets
74     }
75   source.GetZQuench(fZquench);
76   for (Int_t i = 0; i < NTriggerJets(); i++){
77     source.TriggerJet(i,fJets[0][i],fJets[1][i],fJets[2][i],fJets[3][i]);
78   }
79   for (Int_t i = 0; i < NUQTriggerJets(); i++){
80     source.UQJet(i,fUQJets[0][i],fUQJets[1][i],fUQJets[2][i],fUQJets[3][i]);
81   }
82
83   source.Hard(0,fHard[0][0],fHard[1][0],fHard[2][0],fHard[3][0],fHard[4][0]);
84   source.Hard(1,fHard[0][1],fHard[1][1],fHard[2][1],fHard[3][1],fHard[4][1]);
85 }
86
87 AliJetEventParticles::~AliJetEventParticles()
88 {
89   //destructor   
90   Reset();
91   delete fParticles;
92 }
93
94 void  AliJetEventParticles::Reset(Int_t size)
95 {
96   //deletes all particles from the event
97   if(fParticles) fParticles->Clear();
98   if(size>=0) fParticles->Expand(size);
99   fNParticles = 0;
100
101   fVertexX=0.;
102   fVertexY=0.;
103   fVertexZ=0.;
104   fTrials=0;
105   fNJets=0;
106   fNUQJets=0;
107   fXJet=-1;
108   fYJet=-1;
109   fImpact=0.;
110   fNHardScatters=0;
111   fNwNwColl=0;
112   for (Int_t i = 0; i < 4; i++) fZquench[i] = 0.;
113   for (Int_t i = 0; i < 10; i++) 
114     for (Int_t j = 0; j < 4; j++) {
115       fJets[j][i]=0;    // Trigger jets
116       fUQJets[j][i]=0;  // Unquenched trigger jets
117     }
118   for (Int_t i = 0; i < 5; i++){
119     fHard[i][0]=0;
120     fHard[i][1]=0;
121   }
122
123
124 void AliJetEventParticles::AddSignal(const AliJetEventParticles& source)
125 { //mark signal particles and add them to TClonesArray
126   //note that fNParticles still keeps only background particles
127
128   Int_t nSignalParts=source.GetNParticles();
129   for(Int_t i=0; i<nSignalParts; i++)
130     {
131       const AliJetParticle *kjp=source.GetParticle(i);
132
133       AliJetParticle *ap=new((*fParticles)[fNParticles+i]) AliJetParticle(*(kjp));
134       ap->SetType(-123); //mark pythia particle
135     }
136   for(Int_t i=nSignalParts+fNParticles;i<fParticles->GetEntriesFast();i++)
137     fParticles->RemoveAt(i);
138   //cout << fParticles->GetEntries() << " " << fNParticles << " " << nSignalParts << endl;
139
140   /* should we transform the vertex???
141   fVertexX=source.GetVertexX();
142   fVertexY=source.GetVertexY();
143   fVertexZ=source.GetVertexZ();
144   */
145   fTrials=source.Trials();
146   fNJets=source.NTriggerJets();
147   fNUQJets=source.NUQTriggerJets();
148   fXJet=source.GetXJet();
149   fYJet=source.GetXJet();
150
151
152   for (Int_t i = 0; i < 4; i++) fZquench[i] = 0.;
153   for (Int_t i = 0; i < 10; i++) 
154     for (Int_t j = 0; j < 4; j++) {
155       fJets[j][i]=0;    // Trigger jets
156       fUQJets[j][i]=0;  // Unquenched trigger jets
157     }
158   source.GetZQuench(fZquench);
159   for (Int_t i = 0; i < NTriggerJets(); i++){
160     source.TriggerJet(i,fJets[0][i],fJets[1][i],fJets[2][i],fJets[3][i]);
161   }
162   for (Int_t i = 0; i < NUQTriggerJets(); i++){
163     source.UQJet(i,fUQJets[0][i],fUQJets[1][i],fUQJets[2][i],fUQJets[3][i]);
164   }
165   source.Hard(0,fHard[0][0],fHard[1][0],fHard[2][0],fHard[3][0],fHard[4][0]);
166   source.Hard(1,fHard[0][1],fHard[1][1],fHard[2][1],fHard[3][1],fHard[4][1]);
167 }
168
169 void AliJetEventParticles::AddParticle(AliJetParticle* part)
170 {
171   //Adds new particle to the event
172   fParticles->AddAt(part,fNParticles++);
173 }
174
175 void AliJetEventParticles::AddParticle(const AliJetParticle* part)
176 {
177   //Adds new particle to the event
178   new((*fParticles)[fNParticles++]) AliJetParticle(*part);
179 }
180
181 void AliJetEventParticles::AddParticle(const TParticle* part,Int_t idx, Int_t l, Int_t ncl)
182 {
183   //Adds new particle to the event
184   new((*fParticles)[fNParticles++]) AliJetParticle(part,idx,l,ncl);
185 }
186
187 void AliJetEventParticles::AddParticle(Float_t px, Float_t py, Float_t pz, 
188                               Float_t etot, Int_t idx, Int_t l, Int_t ncl)
189 {
190   //Adds new particle to the event
191   new((*fParticles)[fNParticles++]) AliJetParticle(px,py,pz,etot,idx,l,ncl); 
192 }
193
194 void AliJetEventParticles::AddParticle(Float_t px, Float_t py, Float_t pz, Float_t etot, Int_t idx, Int_t l,
195                               Int_t ncl, Float_t pt, Float_t phi, Float_t eta)
196 {
197   //Adds new particle to the event
198   new((*fParticles)[fNParticles++]) AliJetParticle(px,py,pz,etot,idx,l,ncl,pt,phi,eta); 
199 }
200
201 const AliJetParticle* AliJetEventParticles::GetParticleSafely(Int_t n)
202 {
203   //returns nth particle with range check
204   if( (n<0) || (fNParticles<=n) ) return 0;
205   return (const AliJetParticle*)fParticles->At(n);
206 }
207
208 void AliJetEventParticles::AddJet(Float_t px, Float_t py, Float_t pz, Float_t e)
209 {
210   //
211   //  Add a jet 
212   //
213   if (fNJets < 10) {
214     fJets[0][fNJets] = px;
215     fJets[1][fNJets] = py;
216     fJets[2][fNJets] = pz;
217     fJets[3][fNJets] = e;
218     fNJets++;
219   } else {
220     printf("\nWarning: More than 10 jets triggered !!\n");
221   }
222 }
223
224 void AliJetEventParticles::AddJet(Float_t p[4])
225 {
226   //
227   //  Add a jet 
228   //
229   if (fNJets < 10) {
230     fJets[0][fNJets] = p[0];
231     fJets[1][fNJets] = p[1];
232     fJets[2][fNJets] = p[2];
233     fJets[3][fNJets] = p[3];
234     fNJets++;
235   } else {
236     printf("\nWarning: More than 10 jets triggered !!\n");
237   }
238 }
239
240 void AliJetEventParticles::AddUQJet(Float_t px, Float_t py, Float_t pz, Float_t e)
241 {
242   //
243   //  Add a jet 
244   //
245   if (fNUQJets < 10) {
246     fUQJets[0][fNUQJets] = px;
247     fUQJets[1][fNUQJets] = py;
248     fUQJets[2][fNUQJets] = pz;
249     fUQJets[3][fNUQJets] = e;
250     fNUQJets++;
251   } else {
252     printf("\nWarning: More than 10 jets triggered !!\n");
253   }
254 }
255
256 void AliJetEventParticles::AddUQJet(Float_t p[4])
257 {
258   //
259   //  Add a jet 
260   //
261   if (fNUQJets < 10) {
262     fUQJets[0][fNUQJets] = p[0];
263     fUQJets[1][fNUQJets] = p[1];
264     fUQJets[2][fNUQJets] = p[2];
265     fUQJets[3][fNUQJets] = p[3];
266     fNUQJets++;
267   } else {
268     printf("\nWarning: More than 10 jets triggered !!\n");
269   }
270 }
271
272 void AliJetEventParticles::AddHard(Int_t i,Float_t px, Float_t py, Float_t pz, Float_t e, Float_t type)
273 {
274   //
275   //  Add a had parton
276   //
277   if (i < 2) {
278     fHard[0][i] = px;
279     fHard[1][i] = py;
280     fHard[2][i] = pz;
281     fHard[3][i] = e;
282     fHard[4][i] = type;
283   } else {
284     printf("\nWarning: More than 2 partons !!\n");
285   }
286 }
287
288 void AliJetEventParticles::SetZQuench(Double_t z[4])
289 {
290   //
291   // Set quenching fraction
292   //
293   for (Int_t i = 0; i < 4; i++) fZquench[i] = z[i];
294 }
295
296 void AliJetEventParticles::GetZQuench(Double_t z[4]) const
297 {
298   //
299   // Get quenching fraction
300   //
301   for (Int_t i = 0; i < 4; i++) z[i] = fZquench[i];
302 }
303
304 void AliJetEventParticles::TriggerJet(Int_t i, Float_t p[4]) const
305 {
306   //
307   // Give back jet #i
308   //
309   if (i >= fNJets) {
310     printf("\nWarning: TriggerJet, index out of Range!!\n");
311   } else {
312     p[0] = fJets[0][i];
313     p[1] = fJets[1][i];
314     p[2] = fJets[2][i];
315     p[3] = fJets[3][i];
316   }
317 }
318
319 void AliJetEventParticles::TriggerJet(Int_t i, Float_t &p1, Float_t &p2, Float_t &p3, Float_t &E) const
320 {
321   //
322   // Give back jet #i
323   //
324   if (i >= fNJets) {
325     printf("\nWarning: TriggerJet, index out of Range!!\n");
326   } else {
327     p1   = fJets[0][i];
328     p2   = fJets[1][i];
329     p3   = fJets[2][i];
330     E    = fJets[3][i];
331   }
332 }
333
334 void AliJetEventParticles::UQJet(Int_t i, Float_t p[4]) const
335 {
336   //
337   // Give back jet #i
338   //
339   if (i >= fNUQJets) {
340     printf("\nWarning: Unquenched Jets, index out of Range!!\n");
341   } else {
342     p[0] = fUQJets[0][i];
343     p[1] = fUQJets[1][i];
344     p[2] = fUQJets[2][i];
345     p[3] = fUQJets[3][i];
346   }
347 }
348
349 void AliJetEventParticles::UQJet(Int_t i, Float_t &p1, Float_t &p2, Float_t &p3, Float_t &E) const
350 {
351   //
352   // Give back jet #i
353   //
354   if (i >= fNUQJets) {
355     printf("\nWarning: Unquenched Jets, index out of Range!!\n");
356   } else {
357     p1   = fUQJets[0][i];
358     p2   = fUQJets[1][i];
359     p3   = fUQJets[2][i];
360     E    = fUQJets[3][i];
361   }
362 }
363
364 void AliJetEventParticles::Hard(Int_t i, Float_t &p1, Float_t &p2, Float_t &p3, Float_t &E, Float_t &type) const
365 {
366   //
367   // Give back jet #i
368   //
369   if (i >= 2) {
370     printf("\nWarning: Hard partons, index out of Range!!\n");
371   } else {
372     p1   = fHard[0][i];
373     p2   = fHard[1][i];
374     p3   = fHard[2][i];
375     E    = fHard[3][i];
376     type = fHard[4][i];
377   }
378 }
379
380 void AliJetEventParticles::Hard(Int_t i, Float_t p[4], Float_t &type) const
381 {
382   //
383   // Give back jet #i
384   //
385   if (i >= 2) {
386     printf("\nWarning: Hard partons, index out of Range!!\n");
387   } else {
388     p[0]   = fHard[0][i];
389     p[1]   = fHard[1][i];
390     p[2]   = fHard[2][i];
391     p[3]   = fHard[3][i];
392     type = fHard[4][i];
393   }
394 }
395
396 void AliJetEventParticles::SetXYJet(Double_t x, Double_t y)
397 {
398   //
399   //  Add jet production point
400   //
401   fXJet = x; 
402   fYJet = y; 
403 }
404
405 void AliJetEventParticles::Print(Option_t* /*t*/) const
406 {
407   cout << "--- AliJetEventParticles ---" << endl;
408   if(fHeader.Length()) cout << fHeader.Data() << endl;
409   cout << "Particles in Event: " << fNParticles << endl;
410   if(fNUQJets){
411     cout << "Unquenched Jets: " << fNUQJets << endl;
412     for(Int_t i = 0; i<fNUQJets; i++){
413       Float_t x=fUQJets[0][i];
414       Float_t y=fUQJets[1][i];
415       Float_t z=fUQJets[2][i];
416       Float_t e=fUQJets[3][i];
417       Float_t ptj=TMath::Sqrt(x*x+y*y);
418       Float_t thj=TMath::ATan2(ptj,z);
419       Float_t etaj=-TMath::Log(TMath::Tan(thj/2));
420       Float_t phj=TMath::Pi()+TMath::ATan2(-y,-x);
421       Float_t et=e*TMath::Sin(thj);
422       cout << i << " " << et << " " << etaj << " " << phj << endl;
423     }
424   }
425   if(fNJets){
426     cout << "Triggered Jets: " << fNJets << endl;
427     for(Int_t i = 0; i<fNJets; i++){
428       Float_t x=fJets[0][i];
429       Float_t y=fJets[1][i];
430       Float_t z=fJets[2][i];
431       Float_t e=fJets[3][i];
432       Float_t ptj=TMath::Sqrt(x*x+y*y);
433       Float_t thj=TMath::ATan2(ptj,z);
434       Float_t etaj=-TMath::Log(TMath::Tan(thj/2));
435       Float_t phj=TMath::Pi()+TMath::ATan2(-y,-x);
436       Float_t et=e*TMath::Sin(thj);
437       cout << i << " " << et << " " << etaj << " " << phj << endl;
438     }
439   }
440
441 }