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