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