]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - CORRFW/AliCFParticleGenCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Coding rules
[u/mrichter/AliRoot.git] / CORRFW / AliCFParticleGenCuts.cxx
... / ...
CommitLineData
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// ---- CORRECTION FRAMEWORK ----
18// class AliCFParticleGenCuts implementation
19// Using this class a user may define selections relative to
20// MC particle (AliMCParticle) using generation-level information.
21////////////////////////////////////////////////////////////////////////////
22// author : R. Vernet (renaud.vernet@cern.ch)
23////////////////////////////////////////////////////////////////////////////
24
25#include "AliLog.h"
26#include "AliCFParticleGenCuts.h"
27#include "TParticle.h"
28#include "TParticlePDG.h"
29#include "AliMCEvent.h"
30#include "TObject.h"
31#include "AliStack.h"
32#include "TH1F.h"
33#include "TH2F.h"
34#include "TBits.h"
35#include "TList.h"
36#include "TArrayF.h"
37#include "TDecayChannel.h"
38#include "AliAODMCParticle.h"
39#include "AliAODEvent.h"
40
41ClassImp(AliCFParticleGenCuts)
42
43//______________________________
44AliCFParticleGenCuts::AliCFParticleGenCuts() :
45 AliCFCutBase(),
46 fIsAODMC(0),
47 fMCInfo(0x0),
48 fRequireIsCharged(0),
49 fRequireIsNeutral(0),
50 fRequireIsPrimary(0),
51 fRequireIsSecondary(0),
52 fRequirePdgCode(0),
53 fRequireAbsolutePdg(0),
54 fProdVtxRange2D(0),
55 fPdgCode(0),
56 fProdVtxXMin (-1.e+09),
57 fProdVtxYMin (-1.e+09),
58 fProdVtxZMin (-1.e+09),
59 fProdVtxXMax ( 1.e+09),
60 fProdVtxYMax ( 1.e+09),
61 fProdVtxZMax ( 1.e+09),
62 fDecayVtxXMin(-1.e+09),
63 fDecayVtxYMin(-1.e+09),
64 fDecayVtxZMin(-1.e+09),
65 fDecayVtxXMax( 1.e+09),
66 fDecayVtxYMax( 1.e+09),
67 fDecayVtxZMax( 1.e+09),
68 fDecayLengthMin(-1.),
69 fDecayLengthMax(1.e+09),
70 fDecayRxyMin(-1),
71 fDecayRxyMax(1.e+09),
72 fDecayChannel(0x0),
73 fhCutStatistics(0x0),
74 fhCutCorrelation(0x0),
75 fCutValues(new TArrayF(kNCuts)),
76 fBitmap(new TBits(0))
77{
78 //
79 //ctor
80 //
81 for (int i=0; i<kNCuts; i++)
82 for (int j=0; j<kNStepQA; j++)
83 fhQA[i][j]=0x0;
84
85 for (int j=0; j<kNStepQA; j++)
86 fhProdVtxXY[j]=0x0;
87}
88
89//______________________________
90AliCFParticleGenCuts::AliCFParticleGenCuts(const Char_t* name, const Char_t* title) :
91 AliCFCutBase(name,title),
92 fIsAODMC(0),
93 fMCInfo(0x0),
94 fRequireIsCharged(0),
95 fRequireIsNeutral(0),
96 fRequireIsPrimary(0),
97 fRequireIsSecondary(0),
98 fRequirePdgCode(0),
99 fRequireAbsolutePdg(0),
100 fProdVtxRange2D(0),
101 fPdgCode(0),
102 fProdVtxXMin (-1.e+09),
103 fProdVtxYMin (-1.e+09),
104 fProdVtxZMin (-1.e+09),
105 fProdVtxXMax ( 1.e+09),
106 fProdVtxYMax ( 1.e+09),
107 fProdVtxZMax ( 1.e+09),
108 fDecayVtxXMin(-1.e+09),
109 fDecayVtxYMin(-1.e+09),
110 fDecayVtxZMin(-1.e+09),
111 fDecayVtxXMax( 1.e+09),
112 fDecayVtxYMax( 1.e+09),
113 fDecayVtxZMax( 1.e+09),
114 fDecayLengthMin(-1.),
115 fDecayLengthMax(1.e+09),
116 fDecayRxyMin(-1.),
117 fDecayRxyMax(1.e+09),
118 fDecayChannel(0x0),
119 fhCutStatistics(0x0),
120 fhCutCorrelation(0x0),
121 fCutValues(new TArrayF(kNCuts)),
122 fBitmap(new TBits(0))
123{
124 //
125 //ctor
126 //
127 for (int i=0; i<kNCuts; i++)
128 for (int j=0; j<kNStepQA; j++)
129 fhQA[i][j]=0x0;
130
131 for (int j=0; j<kNStepQA; j++)
132 fhProdVtxXY[j]=0x0;
133}
134
135//______________________________
136AliCFParticleGenCuts::AliCFParticleGenCuts(const AliCFParticleGenCuts& c) :
137 AliCFCutBase(c),
138 fIsAODMC(c.fIsAODMC),
139 fMCInfo(c.fMCInfo),
140 fRequireIsCharged(c.fRequireIsCharged),
141 fRequireIsNeutral(c.fRequireIsNeutral),
142 fRequireIsPrimary(c.fRequireIsPrimary),
143 fRequireIsSecondary(c.fRequireIsSecondary),
144 fRequirePdgCode(c.fRequirePdgCode),
145 fRequireAbsolutePdg(c.fRequireAbsolutePdg),
146 fProdVtxRange2D(c.fProdVtxRange2D),
147 fPdgCode(c.fPdgCode),
148 fProdVtxXMin (c.fProdVtxXMin),
149 fProdVtxYMin (c.fProdVtxYMin),
150 fProdVtxZMin (c.fProdVtxZMin),
151 fProdVtxXMax (c.fProdVtxXMax),
152 fProdVtxYMax (c.fProdVtxYMax),
153 fProdVtxZMax (c.fProdVtxZMax),
154 fDecayVtxXMin(c.fDecayVtxXMin),
155 fDecayVtxYMin(c.fDecayVtxYMin),
156 fDecayVtxZMin(c.fDecayVtxZMin),
157 fDecayVtxXMax(c.fDecayVtxXMax),
158 fDecayVtxYMax(c.fDecayVtxYMax),
159 fDecayVtxZMax(c.fDecayVtxZMax),
160 fDecayLengthMin(c.fDecayLengthMin),
161 fDecayLengthMax(c.fDecayLengthMin),
162 fDecayRxyMin(c.fDecayLengthMin),
163 fDecayRxyMax(c.fDecayLengthMin),
164 fDecayChannel(c.fDecayChannel),
165 fhCutStatistics(new TH1F(*c.fhCutStatistics)),
166 fhCutCorrelation(new TH2F(*c.fhCutCorrelation)),
167 fCutValues(new TArrayF(*c.fCutValues)),
168 fBitmap(new TBits(*c.fBitmap))
169{
170 //
171 //copy ctor
172 //
173 for (int i=0; i<kNCuts; i++)
174 for (int j=0; j<kNStepQA; j++)
175 fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
176
177 for (int j=0; j<kNStepQA; j++)
178 fhProdVtxXY[j]=(TH2F*)c.fhProdVtxXY[j]->Clone();
179}
180
181//______________________________
182AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts& c)
183{
184 //
185 // Assignment operator
186 //
187 if (this != &c) {
188 AliCFCutBase::operator=(c) ;
189 fIsAODMC=c.fIsAODMC;
190 fMCInfo=c.fMCInfo;
191 fRequireIsCharged=c.fRequireIsCharged;
192 fRequireIsNeutral=c.fRequireIsNeutral;
193 fRequireIsPrimary=c.fRequireIsPrimary;
194 fRequireIsSecondary=c.fRequireIsSecondary;
195 fRequirePdgCode=c.fRequirePdgCode;
196 fRequireAbsolutePdg=c.fRequireAbsolutePdg;
197 fProdVtxRange2D=c.fProdVtxRange2D;
198 fPdgCode=c.fPdgCode;
199 fProdVtxXMin=c.fProdVtxXMin;
200 fProdVtxYMin=c.fProdVtxYMin;
201 fProdVtxZMin=c.fProdVtxZMin;
202 fProdVtxXMax=c.fProdVtxXMax;
203 fProdVtxYMax=c.fProdVtxYMax;
204 fProdVtxZMax=c.fProdVtxZMax;
205 fDecayVtxXMin=c.fDecayVtxXMin;
206 fDecayVtxYMin=c.fDecayVtxYMin;
207 fDecayVtxZMin=c.fDecayVtxZMin;
208 fDecayVtxXMax=c.fDecayVtxXMax;
209 fDecayVtxYMax=c.fDecayVtxYMax;
210 fDecayVtxZMax=c.fDecayVtxZMax;
211 fDecayLengthMin=c.fDecayVtxZMax;
212 fDecayLengthMax=c.fDecayLengthMax;
213 fDecayRxyMin=c.fDecayRxyMin;
214 fDecayRxyMax=c.fDecayRxyMax;
215 fDecayChannel=c.fDecayChannel;
216 fCutValues=new TArrayF(*c.fCutValues);
217 fBitmap=new TBits(*c.fBitmap);
218
219 if (fhCutStatistics) fhCutStatistics =new TH1F(*c.fhCutStatistics) ;
220 if (fhCutCorrelation) fhCutCorrelation=new TH2F(*c.fhCutCorrelation);
221
222 for (int i=0; i<kNCuts; i++)
223 for (int j=0; j<kNStepQA; j++)
224 fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
225
226 for (int j=0; j<kNStepQA; j++)
227 fhProdVtxXY[j]=(TH2F*)c.fhProdVtxXY[j]->Clone();
228 }
229 return *this ;
230}
231
232//______________________________
233Bool_t AliCFParticleGenCuts::IsSelected(TObject* obj) {
234 //
235 // check if selections on 'obj' are passed
236 // 'obj' must be an AliMCParticle
237 //
238
239 if (!obj) return kFALSE ;
240
241 if (!fIsAODMC) SelectionBitMap((AliMCParticle*) obj);
242 else SelectionBitMap((AliAODMCParticle*)obj);
243
244 if (fIsQAOn) FillHistograms(obj,0);
245
246 for (UInt_t icut=0; icut<fBitmap->GetNbits();icut++)
247 if (!fBitmap->TestBitNumber(icut)) return kFALSE ;
248
249 if (fIsQAOn) FillHistograms(obj,1);
250 return kTRUE;
251}
252
253//__________________________________________________________________________________
254void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart)
255{
256 //
257 // test if the track passes the single cuts
258 // and store the information in a bitmap
259 //
260
261 for (UInt_t i=0; i<kNCuts; i++) {
262 fBitmap->SetBitNumber(i,kFALSE);
263 fCutValues->SetAt((Double32_t)0,i) ;
264 }
265
266 // fill the cut array
267 Double32_t partVx=(Double32_t)mcPart->Xv();
268 Double32_t partVy=(Double32_t)mcPart->Yv();
269 Double32_t partVz=(Double32_t)mcPart->Zv();
270
271 // calculate the production vertex ellipse
272 Double32_t prodVtxXYmin = 0.;
273 if (fProdVtxXMin!=0 && fProdVtxYMin!=0)
274 prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
275 Double32_t prodVtxXYmax = 0.;
276 if(fProdVtxXMax!=0 && fProdVtxYMax!=0)
277 prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
278
279 Double32_t decayVx=0.;
280 Double32_t decayVy=0.;
281 Double32_t decayVz=0.;
282 Double32_t decayL=0.;
283 Double32_t decayRxy=0.;
284
285 TParticle* part = mcPart->Particle();
286 AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
287 TParticle* daughter=0x0;
288 if ( part->GetNDaughters() > 0 ) {
289 daughter = stack->Particle(part->GetFirstDaughter()) ;
290 decayVx=(Double32_t)daughter->Vx();
291 decayVy=(Double32_t)daughter->Vy();
292 decayVz=(Double32_t)daughter->Vz();
293 decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) +
294 TMath::Power(partVy-decayVy,2) +
295 TMath::Power(partVz-decayVz,2) ) ;
296 decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
297 }
298
299 fCutValues->SetAt(partVx ,kCutProdVtxXMin) ;
300 fCutValues->SetAt(partVy ,kCutProdVtxYMin) ;
301 fCutValues->SetAt(partVz ,kCutProdVtxZMin) ;
302 fCutValues->SetAt(partVx ,kCutProdVtxXMax) ;
303 fCutValues->SetAt(partVy ,kCutProdVtxYMax) ;
304 fCutValues->SetAt(partVz ,kCutProdVtxZMax) ;
305 fCutValues->SetAt(decayVx ,kCutDecVtxXMin) ;
306 fCutValues->SetAt(decayVy ,kCutDecVtxYMin) ;
307 fCutValues->SetAt(decayVz ,kCutDecVtxZMin) ;
308 fCutValues->SetAt(decayVx ,kCutDecVtxXMax) ;
309 fCutValues->SetAt(decayVy ,kCutDecVtxYMax) ;
310 fCutValues->SetAt(decayVz ,kCutDecVtxZMax) ;
311 fCutValues->SetAt(decayL ,kCutDecLgthMin) ;
312 fCutValues->SetAt(decayL ,kCutDecLgthMax) ;
313 fCutValues->SetAt(decayRxy,kCutDecRxyMin) ;
314 fCutValues->SetAt(decayRxy,kCutDecRxyMax) ;
315
316 // cut on charge
317 if ( fRequireIsCharged || fRequireIsNeutral ) {
318 if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
319 if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
320 }
321 else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
322
323 // cut on primary/secondary
324 if ( fRequireIsPrimary || fRequireIsSecondary) {
325 if (fRequireIsPrimary && IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
326 if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
327 }
328 else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
329
330 // cut on PDG code
331 if ( fRequirePdgCode ) {
332 if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
333 }
334 else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
335
336 // cut on decay channel
337 if ( fDecayChannel ) {
338 Bool_t goodDecay = kTRUE ;
339 Short_t nDaughters = mcPart->Particle()->GetNDaughters() ;
340 if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
341 //now number of daughters is OK
342 if (goodDecay) {
343 // now check if decay channel is respected
344 // first try
345 for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
346 TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
347 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
348 }
349 if (!goodDecay) {
350 //second try inverting the order of the daughters
351 goodDecay = kTRUE ;
352 for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
353 TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
354 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
355 }
356 }
357 if (!goodDecay && fRequireAbsolutePdg) {
358 //now tries inverting the sign of the daughters in case the anti-particle is also looked at
359 // third try
360 goodDecay = kTRUE ;
361 for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
362 TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
363 if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
364 }
365 if (!goodDecay) {
366 //fourth try inverting the order of the daughters
367 goodDecay = kTRUE ;
368 for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
369 TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
370 if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
371 }
372 }
373 } //end check anti-particle
374 } //end # daughters OK
375 fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
376 } //end require decay channel
377 else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
378
379
380 // now array of cut is build, fill the bitmap consequently
381 Int_t iCutBit = -1;
382 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
383 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
384 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
385
386 ++iCutBit;
387 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
388 || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
389
390 ++iCutBit;
391 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
392 || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
393
394 ++iCutBit;
395 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
396 || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
397
398 ++iCutBit;
399 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
400 || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
401
402 if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
403 if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
404 if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
405 if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
406 if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
407 if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
408 if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
409 if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
410 if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
411 if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
412 if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
413 if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
414 if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
415}
416
417//__________________________________________________________________________________
418void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
419{
420 //
421 // test if the track passes the single cuts
422 // and store the information in a bitmap
423 //
424
425 for (UInt_t i=0; i<kNCuts; i++) {
426 fBitmap->SetBitNumber(i,kFALSE);
427 fCutValues->SetAt((Double32_t)0,i) ;
428 }
429
430 // fill the cut array
431 Double32_t partVx=(Double32_t)mcPart->Xv();
432 Double32_t partVy=(Double32_t)mcPart->Yv();
433 Double32_t partVz=(Double32_t)mcPart->Zv();
434
435 // calculate the production vertex ellipse
436 Double32_t prodVtxXYmin = 0.;
437 if (fProdVtxXMin!=0 && fProdVtxYMin!=0)
438 prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
439 Double32_t prodVtxXYmax = 0.;
440 if(fProdVtxXMax!=0 && fProdVtxYMax!=0)
441 prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
442
443 Double32_t decayVx=0.;
444 Double32_t decayVy=0.;
445 Double32_t decayVz=0.;
446 Double32_t decayL=0.;
447 Double32_t decayRxy=0.;
448
449 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fMCInfo);
450 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
451 AliAODMCParticle* daughter = 0x0 ;
452
453 if (mcPart->GetDaughter(0)>0) {
454 daughter = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
455
456 decayVx=(Double32_t)daughter->Xv();
457 decayVy=(Double32_t)daughter->Yv();
458 decayVz=(Double32_t)daughter->Zv();
459 decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) +
460 TMath::Power(partVy-decayVy,2) +
461 TMath::Power(partVz-decayVz,2) ) ;
462 decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
463
464 }
465
466 fCutValues->SetAt(partVx ,kCutProdVtxXMin) ;
467 fCutValues->SetAt(partVy ,kCutProdVtxYMin) ;
468 fCutValues->SetAt(partVz ,kCutProdVtxZMin) ;
469 fCutValues->SetAt(partVx ,kCutProdVtxXMax) ;
470 fCutValues->SetAt(partVy ,kCutProdVtxYMax) ;
471 fCutValues->SetAt(partVz ,kCutProdVtxZMax) ;
472 fCutValues->SetAt(decayVx ,kCutDecVtxXMin) ;
473 fCutValues->SetAt(decayVy ,kCutDecVtxYMin) ;
474 fCutValues->SetAt(decayVz ,kCutDecVtxZMin) ;
475 fCutValues->SetAt(decayVx ,kCutDecVtxXMax) ;
476 fCutValues->SetAt(decayVy ,kCutDecVtxYMax) ;
477 fCutValues->SetAt(decayVz ,kCutDecVtxZMax) ;
478 fCutValues->SetAt(decayL ,kCutDecLgthMin) ;
479 fCutValues->SetAt(decayL ,kCutDecLgthMax) ;
480 fCutValues->SetAt(decayRxy,kCutDecRxyMin) ;
481 fCutValues->SetAt(decayRxy,kCutDecRxyMax) ;
482
483 // cut on charge
484 if ( fRequireIsCharged || fRequireIsNeutral ) {
485 if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
486 if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
487 }
488 else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
489
490 // cut on primary/secondary
491 if ( fRequireIsPrimary || fRequireIsSecondary) {
492 if (fRequireIsPrimary && IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
493 if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
494 }
495 else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
496
497 // cut on PDG code
498 if ( fRequirePdgCode ) {
499 if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
500 }
501 else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
502
503 // cut on decay channel
504 if ( fDecayChannel ) {
505 Bool_t goodDecay = kTRUE ;
506 Short_t nDaughters = 0 ;
507 if (mcPart->GetDaughter(0) >=0) nDaughters = 1 + mcPart->GetDaughter(1) - mcPart->GetDaughter(0);
508
509 if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
510 if (goodDecay) {
511 // now check if decay channel is respected
512 // first try
513 for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
514 AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)+iDaughter));
515 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
516 }
517 if (!goodDecay) {
518 //second try inverting the order of the daughters
519 goodDecay = kTRUE ;
520 for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
521 AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter));
522 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
523 }
524 }
525 }
526 fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
527 }
528 else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
529
530
531 // now array of cut is build, fill the bitmap consequently
532 Int_t iCutBit = -1;
533 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
534 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
535 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
536
537 ++iCutBit;
538 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
539 || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
540
541 ++iCutBit;
542 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
543 || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
544
545 ++iCutBit;
546 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
547 || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
548
549 ++iCutBit;
550 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
551 || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
552
553 if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
554 if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
555 if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
556 if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
557 if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
558 if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
559 if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
560 if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
561 if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
562 if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
563 if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
564 if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
565 if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
566}
567
568
569//__________________________________________________________________________________
570void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
571{
572 //
573 // fill the QA histograms
574 //
575
576 for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++)
577 fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
578
579 fhProdVtxXY[afterCuts]->Fill(fCutValues->At(4),fCutValues->At(5));
580
581 // fill cut statistics and cut correlation histograms with information from the bitmap
582 if (afterCuts) return;
583
584 // Number of single cuts in this class
585 UInt_t ncuts = fBitmap->GetNbits();
586 for(UInt_t bit=0; bit<ncuts;bit++) {
587 if (!fBitmap->TestBitNumber(bit)) {
588 fhCutStatistics->Fill(bit+1);
589 for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
590 if (!fBitmap->TestBitNumber(bit2))
591 fhCutCorrelation->Fill(bit+1,bit2+1);
592 }
593 }
594 }
595}
596
597//__________________________________________________________________________________
598void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
599 //
600 // saves the histograms in a TList
601 //
602
603 DefineHistograms();
604
605 qaList->Add(fhCutStatistics);
606 qaList->Add(fhCutCorrelation);
607
608 for (Int_t j=0; j<kNStepQA; j++) {
609 qaList->Add(fhProdVtxXY[j]);
610 for(Int_t i=0; i<kNCuts; i++)
611 qaList->Add(fhQA[i][j]);
612 }
613}
614
615//__________________________________________________________________________________
616void AliCFParticleGenCuts::DefineHistograms() {
617 //
618 // histograms for cut variables, cut statistics and cut correlations
619 //
620 Int_t color = 2;
621
622 // book cut statistics and cut correlation histograms
623 fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
624 fhCutStatistics->SetLineWidth(2);
625 int k = 1;
626 fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge") ; k++;
627 fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec") ; k++;
628 fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG") ; k++;
629 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin") ; k++;
630 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax") ; k++;
631 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin") ; k++;
632 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax") ; k++;
633 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin") ; k++;
634 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++;
635 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin") ; k++;
636 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax") ; k++;
637 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin") ; k++;
638 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax") ; k++;
639 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin") ; k++;
640 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++;
641 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
642 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
643 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin") ; k++;
644 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax") ; k++;
645 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
646
647
648 fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
649 fhCutCorrelation->SetLineWidth(2);
650 for (k=1; k<=kNCuts; k++) {
651 fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
652 fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
653 }
654
655 Char_t str[256];
656 for (int i=0; i<kNStepQA; i++) {
657 if (i==0) sprintf(str," ");
658 else sprintf(str,"_cut");
659 fhQA[kCutCharge] [i] = new TH1F(Form("%s_charge%s" ,GetName(),str),"",2,0,2);
660 fhQA[kCutPrimSec] [i] = new TH1F(Form("%s_primSec%s" ,GetName(),str),"",2,0,2);
661 fhQA[kCutPDGCode] [i] = new TH1F(Form("%s_pdgCode%s" ,GetName(),str),"",2,0,2);
662 fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,-10,10);
663 fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,-10,10);
664 fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,-10,10);
665 fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,-10,10);
666 fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,-10,10);
667 fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,-10,10);
668 fhQA[kCutDecVtxXMin] [i] = new TH1F(Form("%s_decVtxXMin%s" ,GetName(),str),"",100,0,10);
669 fhQA[kCutDecVtxXMax] [i] = new TH1F(Form("%s_decVtxXMax%s" ,GetName(),str),"",100,0,10);
670 fhQA[kCutDecVtxYMin] [i] = new TH1F(Form("%s_decVtxYMin%s" ,GetName(),str),"",100,0,10);
671 fhQA[kCutDecVtxYMax] [i] = new TH1F(Form("%s_decVtxYMax%s" ,GetName(),str),"",100,0,10);
672 fhQA[kCutDecVtxZMin] [i] = new TH1F(Form("%s_decVtxZMin%s" ,GetName(),str),"",100,0,10);
673 fhQA[kCutDecVtxZMax] [i] = new TH1F(Form("%s_decVtxZMax%s" ,GetName(),str),"",100,0,10);
674 fhQA[kCutDecLgthMin] [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
675 fhQA[kCutDecLgthMax] [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
676 fhQA[kCutDecRxyMin] [i] = new TH1F(Form("%s_decRxyMin%s" ,GetName(),str),"",100,0,10);
677 fhQA[kCutDecRxyMax] [i] = new TH1F(Form("%s_decRxyMax%s" ,GetName(),str),"",100,0,10);
678 fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
679 fhProdVtxXY [i] = new TH2F(Form("%s_prodVtxXY%s" ,GetName(),str),"",100,0,10,100,0,10);
680 fhProdVtxXY [i] ->GetXaxis()->SetTitle("x_{production vertex}");
681 fhProdVtxXY [i] ->GetYaxis()->SetTitle("y_{production vertex}");
682 fhQA[kCutProdVtxXMax] [i] ->GetXaxis()->SetTitle("x_{production vertex}");
683 fhQA[kCutProdVtxYMax] [i] ->GetXaxis()->SetTitle("y_{production vertex}");
684 }
685 for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
686}
687
688
689//______________________________
690Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) {
691 //
692 //check if particle is charged.
693 //
694 if (mcPart->Charge()==0) return kFALSE;
695 return kTRUE;
696}
697//______________________________
698Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) {
699 //
700 //check if particle is primary (standard definition)
701 //
702
703 AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
704
705 if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE;
706 return kTRUE;
707}
708//______________________________
709Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) {
710 //
711 //check if particle is primary (standard definition)
712 //
713
714 if (!mcPart->IsPhysicalPrimary()) return kFALSE;
715 return kTRUE;
716}
717//______________________________
718Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
719 //
720 //check if a charged particle is primary (standard definition)
721 //
722
723 if (!fIsAODMC) {
724 if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
725 }
726 else {
727 if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
728 }
729 return kTRUE;
730}
731//______________________________
732Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
733 //
734 //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
735 //absolute value.
736 //
737 TParticle* part = mcPart->Particle();
738 Int_t pdgCode = part->GetPdgCode();
739
740 if (abs) {
741 pdgCode = TMath::Abs(pdgCode);
742 pdg = TMath::Abs(pdg);
743 }
744 if (pdgCode != pdg ) return kFALSE;
745 return kTRUE;
746}
747//______________________________
748Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
749 //
750 //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
751 //absolute value.
752 //
753 Int_t pdgCode = mcPart->GetPdgCode();
754
755 if (abs) {
756 pdgCode = TMath::Abs(pdgCode);
757 pdg = TMath::Abs(pdg);
758 }
759 if (pdgCode != pdg ) return kFALSE;
760 return kTRUE;
761}
762//______________________________
763void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) {
764 //
765 // Sets pointer to MC event information (AliMCEvent)
766 //
767
768 if (!mcEvent) {
769 AliError("Pointer to MC Event is null !");
770 return;
771 }
772
773 TString className(mcEvent->ClassName());
774 if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) {
775 AliError("argument must point to an AliMCEvent or an AliAODEvent !");
776 return ;
777 }
778
779 if (fIsAODMC) fMCInfo = dynamic_cast<AliAODEvent*>(mcEvent) ;
780 else fMCInfo = dynamic_cast<AliMCEvent*> (mcEvent) ;
781}