]> 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
Secure codin - overruns corrected.
[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//______________________________
9eeae5d5 233Bool_t AliCFParticleGenCuts::IsSelected(TObject* obj) {
563113d0 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
34911148 272 Double32_t prodVtxXYmin = 0.;
55003684 273 if (fProdVtxXMin>0 && fProdVtxYMin>0)
34911148 274 prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
275 Double32_t prodVtxXYmax = 0.;
55003684 276 if(fProdVtxXMax>0 && fProdVtxYMax>0)
34911148 277 prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
004f223f 278
107a3100 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
b95f6a36 285 TParticle* part = mcPart->Particle();
286 AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
287 TParticle* daughter=0x0;
107a3100 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
c7803356 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
107a3100 316 // cut on charge
317 if ( fRequireIsCharged || fRequireIsNeutral ) {
c7803356 318 if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
319 if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
563113d0 320 }
c7803356 321 else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
563113d0 322
107a3100 323 // cut on primary/secondary
324 if ( fRequireIsPrimary || fRequireIsSecondary) {
b95f6a36 325 if (fRequireIsPrimary && IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
326 if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
563113d0 327 }
c7803356 328 else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
329
107a3100 330 // cut on PDG code
c7803356 331 if ( fRequirePdgCode ) {
3ff4a092 332 if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
563113d0 333 }
c7803356 334 else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
107a3100 335
8fe1a43d 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 ;
2f8e8cad 341 //now number of daughters is OK
8fe1a43d 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 }
2f8e8cad 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
8fe1a43d 375 fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
2f8e8cad 376 } //end require decay channel
8fe1a43d 377 else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
378
107a3100 379
c7803356 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);
004f223f 385
386 ++iCutBit;
387 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
55003684 388 || ( fProdVtxRange2D && (fProdVtxXMin>0 && fProdVtxYMin>0) && prodVtxXYmin >= 1)
389 || ( fProdVtxRange2D && (fProdVtxXMin<=0 || fProdVtxYMin<=0) ) )
390 fBitmap->SetBitNumber(iCutBit,kTRUE);
004f223f 391
392 ++iCutBit;
393 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
55003684 394 || ( fProdVtxRange2D && (fProdVtxXMax>0 && fProdVtxYMax>0) && prodVtxXYmax <= 1)
395 || ( fProdVtxRange2D && (fProdVtxXMax<=0 || fProdVtxYMax<=0) ) )
396 fBitmap->SetBitNumber(iCutBit,kTRUE);
004f223f 397
398 ++iCutBit;
399 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
55003684 400 || ( fProdVtxRange2D && (fProdVtxXMin>0 && fProdVtxYMin>0) && prodVtxXYmin >= 1)
401 || ( fProdVtxRange2D && (fProdVtxXMin<=0 || fProdVtxYMin<=0) ) )
402 fBitmap->SetBitNumber(iCutBit,kTRUE);
004f223f 403
404 ++iCutBit;
405 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
55003684 406 || ( fProdVtxRange2D && (fProdVtxXMax>0 && fProdVtxYMax>0) && prodVtxXYmax <= 1)
407 || ( fProdVtxRange2D && (fProdVtxXMax<=0 || fProdVtxYMax<=0) ) )
408 fBitmap->SetBitNumber(iCutBit,kTRUE);
004f223f 409
c7803356 410 if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
411 if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
412 if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
413 if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
414 if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
415 if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
416 if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
417 if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
418 if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
419 if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
420 if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
421 if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
8fe1a43d 422 if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
107a3100 423}
424
b95f6a36 425//__________________________________________________________________________________
426void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
427{
428 //
429 // test if the track passes the single cuts
430 // and store the information in a bitmap
431 //
432
433 for (UInt_t i=0; i<kNCuts; i++) {
434 fBitmap->SetBitNumber(i,kFALSE);
435 fCutValues->SetAt((Double32_t)0,i) ;
436 }
437
438 // fill the cut array
439 Double32_t partVx=(Double32_t)mcPart->Xv();
440 Double32_t partVy=(Double32_t)mcPart->Yv();
441 Double32_t partVz=(Double32_t)mcPart->Zv();
442
004f223f 443 // calculate the production vertex ellipse
34911148 444 Double32_t prodVtxXYmin = 0.;
445 if (fProdVtxXMin!=0 && fProdVtxYMin!=0)
446 prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
447 Double32_t prodVtxXYmax = 0.;
448 if(fProdVtxXMax!=0 && fProdVtxYMax!=0)
449 prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
004f223f 450
b95f6a36 451 Double32_t decayVx=0.;
452 Double32_t decayVy=0.;
453 Double32_t decayVz=0.;
454 Double32_t decayL=0.;
455 Double32_t decayRxy=0.;
456
457 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fMCInfo);
4bc0c5f9 458
459 if (!aod) {
460 AliError("AOD event casting failed");
461 return;
462 }
463
b95f6a36 464 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
4bc0c5f9 465 if (!mcArray) {
466 AliError("array casting failed");
467 return;
468 }
b95f6a36 469 AliAODMCParticle* daughter = 0x0 ;
470
471 if (mcPart->GetDaughter(0)>0) {
472 daughter = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
4bc0c5f9 473 if (!daughter) {
474 AliError("daughter casting failed");
475 return;
476 }
b95f6a36 477
478 decayVx=(Double32_t)daughter->Xv();
479 decayVy=(Double32_t)daughter->Yv();
480 decayVz=(Double32_t)daughter->Zv();
481 decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) +
482 TMath::Power(partVy-decayVy,2) +
483 TMath::Power(partVz-decayVz,2) ) ;
484 decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
485
486 }
487
488 fCutValues->SetAt(partVx ,kCutProdVtxXMin) ;
489 fCutValues->SetAt(partVy ,kCutProdVtxYMin) ;
490 fCutValues->SetAt(partVz ,kCutProdVtxZMin) ;
491 fCutValues->SetAt(partVx ,kCutProdVtxXMax) ;
492 fCutValues->SetAt(partVy ,kCutProdVtxYMax) ;
493 fCutValues->SetAt(partVz ,kCutProdVtxZMax) ;
494 fCutValues->SetAt(decayVx ,kCutDecVtxXMin) ;
495 fCutValues->SetAt(decayVy ,kCutDecVtxYMin) ;
496 fCutValues->SetAt(decayVz ,kCutDecVtxZMin) ;
497 fCutValues->SetAt(decayVx ,kCutDecVtxXMax) ;
498 fCutValues->SetAt(decayVy ,kCutDecVtxYMax) ;
499 fCutValues->SetAt(decayVz ,kCutDecVtxZMax) ;
500 fCutValues->SetAt(decayL ,kCutDecLgthMin) ;
501 fCutValues->SetAt(decayL ,kCutDecLgthMax) ;
502 fCutValues->SetAt(decayRxy,kCutDecRxyMin) ;
503 fCutValues->SetAt(decayRxy,kCutDecRxyMax) ;
504
505 // cut on charge
506 if ( fRequireIsCharged || fRequireIsNeutral ) {
507 if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
508 if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
509 }
510 else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
511
512 // cut on primary/secondary
513 if ( fRequireIsPrimary || fRequireIsSecondary) {
514 if (fRequireIsPrimary && IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
515 if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
516 }
517 else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
518
519 // cut on PDG code
520 if ( fRequirePdgCode ) {
3ff4a092 521 if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
b95f6a36 522 }
523 else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
524
525 // cut on decay channel
526 if ( fDecayChannel ) {
527 Bool_t goodDecay = kTRUE ;
528 Short_t nDaughters = 0 ;
529 if (mcPart->GetDaughter(0) >=0) nDaughters = 1 + mcPart->GetDaughter(1) - mcPart->GetDaughter(0);
530
531 if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
532 if (goodDecay) {
533 // now check if decay channel is respected
534 // first try
535 for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
536 AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)+iDaughter));
b7ccf481 537 if (!daug) {
538 AliError("daughter: casting failed");
539 continue;
540 }
b95f6a36 541 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
542 }
543 if (!goodDecay) {
544 //second try inverting the order of the daughters
545 goodDecay = kTRUE ;
546 for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
547 AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter));
548 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
549 }
550 }
551 }
552 fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
553 }
554 else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
555
556
557 // now array of cut is build, fill the bitmap consequently
558 Int_t iCutBit = -1;
559 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
560 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
561 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
004f223f 562
563 ++iCutBit;
564 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
565 || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
566
567 ++iCutBit;
568 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
569 || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
570
571 ++iCutBit;
572 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
573 || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
574
575 ++iCutBit;
576 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
577 || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
578
b95f6a36 579 if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
580 if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
581 if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
582 if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
583 if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
584 if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
585 if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
586 if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
587 if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
588 if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
589 if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
590 if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
591 if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
592}
593
c7803356 594
107a3100 595//__________________________________________________________________________________
033789c9 596void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
107a3100 597{
598 //
599 // fill the QA histograms
600 //
601
602 for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++)
c7803356 603 fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
107a3100 604
004f223f 605 fhProdVtxXY[afterCuts]->Fill(fCutValues->At(4),fCutValues->At(5));
606
107a3100 607 // fill cut statistics and cut correlation histograms with information from the bitmap
608 if (afterCuts) return;
609
610 // Number of single cuts in this class
611 UInt_t ncuts = fBitmap->GetNbits();
612 for(UInt_t bit=0; bit<ncuts;bit++) {
c7803356 613 if (!fBitmap->TestBitNumber(bit)) {
107a3100 614 fhCutStatistics->Fill(bit+1);
615 for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
616 if (!fBitmap->TestBitNumber(bit2))
617 fhCutCorrelation->Fill(bit+1,bit2+1);
618 }
619 }
563113d0 620 }
107a3100 621}
622
623//__________________________________________________________________________________
624void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
625 //
626 // saves the histograms in a TList
627 //
563113d0 628
107a3100 629 DefineHistograms();
563113d0 630
107a3100 631 qaList->Add(fhCutStatistics);
632 qaList->Add(fhCutCorrelation);
633
634 for (Int_t j=0; j<kNStepQA; j++) {
004f223f 635 qaList->Add(fhProdVtxXY[j]);
107a3100 636 for(Int_t i=0; i<kNCuts; i++)
637 qaList->Add(fhQA[i][j]);
638 }
639}
640
641//__________________________________________________________________________________
642void AliCFParticleGenCuts::DefineHistograms() {
643 //
644 // histograms for cut variables, cut statistics and cut correlations
645 //
646 Int_t color = 2;
647
648 // book cut statistics and cut correlation histograms
649 fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
650 fhCutStatistics->SetLineWidth(2);
651 int k = 1;
652 fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge") ; k++;
653 fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec") ; k++;
654 fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG") ; k++;
655 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin") ; k++;
656 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax") ; k++;
657 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin") ; k++;
658 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax") ; k++;
659 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin") ; k++;
660 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++;
661 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin") ; k++;
662 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax") ; k++;
663 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin") ; k++;
664 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax") ; k++;
665 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin") ; k++;
666 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++;
667 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
668 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
669 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin") ; k++;
670 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax") ; k++;
8fe1a43d 671 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
107a3100 672
673
674 fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
675 fhCutCorrelation->SetLineWidth(2);
676 for (k=1; k<=kNCuts; k++) {
677 fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
678 fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
679 }
680
4bc0c5f9 681 Char_t str[5];
107a3100 682 for (int i=0; i<kNStepQA; i++) {
4bc0c5f9 683 if (i==0) snprintf(str,5," ");
684 else snprintf(str,5,"_cut");
8fe1a43d 685 fhQA[kCutCharge] [i] = new TH1F(Form("%s_charge%s" ,GetName(),str),"",2,0,2);
686 fhQA[kCutPrimSec] [i] = new TH1F(Form("%s_primSec%s" ,GetName(),str),"",2,0,2);
687 fhQA[kCutPDGCode] [i] = new TH1F(Form("%s_pdgCode%s" ,GetName(),str),"",2,0,2);
34911148 688 fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,-10,10);
689 fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,-10,10);
690 fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,-10,10);
691 fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,-10,10);
692 fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,-10,10);
693 fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,-10,10);
8fe1a43d 694 fhQA[kCutDecVtxXMin] [i] = new TH1F(Form("%s_decVtxXMin%s" ,GetName(),str),"",100,0,10);
695 fhQA[kCutDecVtxXMax] [i] = new TH1F(Form("%s_decVtxXMax%s" ,GetName(),str),"",100,0,10);
696 fhQA[kCutDecVtxYMin] [i] = new TH1F(Form("%s_decVtxYMin%s" ,GetName(),str),"",100,0,10);
697 fhQA[kCutDecVtxYMax] [i] = new TH1F(Form("%s_decVtxYMax%s" ,GetName(),str),"",100,0,10);
698 fhQA[kCutDecVtxZMin] [i] = new TH1F(Form("%s_decVtxZMin%s" ,GetName(),str),"",100,0,10);
699 fhQA[kCutDecVtxZMax] [i] = new TH1F(Form("%s_decVtxZMax%s" ,GetName(),str),"",100,0,10);
700 fhQA[kCutDecLgthMin] [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
701 fhQA[kCutDecLgthMax] [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
702 fhQA[kCutDecRxyMin] [i] = new TH1F(Form("%s_decRxyMin%s" ,GetName(),str),"",100,0,10);
703 fhQA[kCutDecRxyMax] [i] = new TH1F(Form("%s_decRxyMax%s" ,GetName(),str),"",100,0,10);
704 fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
004f223f 705 fhProdVtxXY [i] = new TH2F(Form("%s_prodVtxXY%s" ,GetName(),str),"",100,0,10,100,0,10);
706 fhProdVtxXY [i] ->GetXaxis()->SetTitle("x_{production vertex}");
707 fhProdVtxXY [i] ->GetYaxis()->SetTitle("y_{production vertex}");
708 fhQA[kCutProdVtxXMax] [i] ->GetXaxis()->SetTitle("x_{production vertex}");
709 fhQA[kCutProdVtxYMax] [i] ->GetXaxis()->SetTitle("y_{production vertex}");
107a3100 710 }
711 for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
563113d0 712}
107a3100 713
714
563113d0 715//______________________________
b95f6a36 716Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) {
563113d0 717 //
718 //check if particle is charged.
719 //
b95f6a36 720 if (mcPart->Charge()==0) return kFALSE;
563113d0 721 return kTRUE;
722}
723//______________________________
b95f6a36 724Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) {
563113d0 725 //
726 //check if particle is primary (standard definition)
727 //
b95f6a36 728
729 AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
730
731 if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE;
563113d0 732 return kTRUE;
733}
734//______________________________
b95f6a36 735Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) {
736 //
737 //check if particle is primary (standard definition)
738 //
739
740 if (!mcPart->IsPhysicalPrimary()) return kFALSE;
741 return kTRUE;
742}
743//______________________________
744Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
563113d0 745 //
746 //check if a charged particle is primary (standard definition)
747 //
b95f6a36 748
749 if (!fIsAODMC) {
750 if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
751 }
752 else {
753 if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
754 }
563113d0 755 return kTRUE;
756}
757//______________________________
758Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
759 //
760 //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
3ff4a092 761 //absolute value.
563113d0 762 //
763 TParticle* part = mcPart->Particle();
764 Int_t pdgCode = part->GetPdgCode();
b95f6a36 765
766 if (abs) {
767 pdgCode = TMath::Abs(pdgCode);
768 pdg = TMath::Abs(pdg);
769 }
770 if (pdgCode != pdg ) return kFALSE;
771 return kTRUE;
772}
773//______________________________
774Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
775 //
776 //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
3ff4a092 777 //absolute value.
b95f6a36 778 //
779 Int_t pdgCode = mcPart->GetPdgCode();
780
fc01457a 781 if (abs) {
782 pdgCode = TMath::Abs(pdgCode);
783 pdg = TMath::Abs(pdg);
784 }
107a3100 785 if (pdgCode != pdg ) return kFALSE;
563113d0 786 return kTRUE;
787}
788//______________________________
0c01ae65 789void AliCFParticleGenCuts::SetMCEventInfo(const TObject* mcEvent) {
563113d0 790 //
107a3100 791 // Sets pointer to MC event information (AliMCEvent)
563113d0 792 //
793
107a3100 794 if (!mcEvent) {
795 AliError("Pointer to MC Event is null !");
563113d0 796 return;
797 }
798
107a3100 799 TString className(mcEvent->ClassName());
b95f6a36 800 if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) {
801 AliError("argument must point to an AliMCEvent or an AliAODEvent !");
563113d0 802 return ;
803 }
0c01ae65 804
805 fMCInfo = (AliVEvent*)mcEvent ;
563113d0 806}