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