]> 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
move cluster position correction to the AliTRDcluster class
[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 ;
317 if (goodDecay) {
318 // now check if decay channel is respected
319 // first try
320 for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
321 TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
322 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
323 }
324 if (!goodDecay) {
325 //second try inverting the order of the daughters
326 goodDecay = kTRUE ;
327 for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
328 TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
329 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
330 }
331 }
332 }
333 fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
334 }
335 else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
336
107a3100 337
c7803356 338 // now array of cut is build, fill the bitmap consequently
339 Int_t iCutBit = -1;
340 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
341 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
342 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
343 if ( fCutValues->At(++iCutBit) > fProdVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
344 if ( fCutValues->At(++iCutBit) < fProdVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
345 if ( fCutValues->At(++iCutBit) > fProdVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
346 if ( fCutValues->At(++iCutBit) < fProdVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
347 if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
348 if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
349 if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
350 if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
351 if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
352 if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
353 if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
354 if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
355 if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
356 if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
357 if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
358 if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
8fe1a43d 359 if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
107a3100 360}
361
b95f6a36 362//__________________________________________________________________________________
363void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
364{
365 //
366 // test if the track passes the single cuts
367 // and store the information in a bitmap
368 //
369
370 for (UInt_t i=0; i<kNCuts; i++) {
371 fBitmap->SetBitNumber(i,kFALSE);
372 fCutValues->SetAt((Double32_t)0,i) ;
373 }
374
375 // fill the cut array
376 Double32_t partVx=(Double32_t)mcPart->Xv();
377 Double32_t partVy=(Double32_t)mcPart->Yv();
378 Double32_t partVz=(Double32_t)mcPart->Zv();
379
380 Double32_t decayVx=0.;
381 Double32_t decayVy=0.;
382 Double32_t decayVz=0.;
383 Double32_t decayL=0.;
384 Double32_t decayRxy=0.;
385
386 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fMCInfo);
387 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
388 AliAODMCParticle* daughter = 0x0 ;
389
390 if (mcPart->GetDaughter(0)>0) {
391 daughter = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
392
393 decayVx=(Double32_t)daughter->Xv();
394 decayVy=(Double32_t)daughter->Yv();
395 decayVz=(Double32_t)daughter->Zv();
396 decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) +
397 TMath::Power(partVy-decayVy,2) +
398 TMath::Power(partVz-decayVz,2) ) ;
399 decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
400
401 }
402
403 fCutValues->SetAt(partVx ,kCutProdVtxXMin) ;
404 fCutValues->SetAt(partVy ,kCutProdVtxYMin) ;
405 fCutValues->SetAt(partVz ,kCutProdVtxZMin) ;
406 fCutValues->SetAt(partVx ,kCutProdVtxXMax) ;
407 fCutValues->SetAt(partVy ,kCutProdVtxYMax) ;
408 fCutValues->SetAt(partVz ,kCutProdVtxZMax) ;
409 fCutValues->SetAt(decayVx ,kCutDecVtxXMin) ;
410 fCutValues->SetAt(decayVy ,kCutDecVtxYMin) ;
411 fCutValues->SetAt(decayVz ,kCutDecVtxZMin) ;
412 fCutValues->SetAt(decayVx ,kCutDecVtxXMax) ;
413 fCutValues->SetAt(decayVy ,kCutDecVtxYMax) ;
414 fCutValues->SetAt(decayVz ,kCutDecVtxZMax) ;
415 fCutValues->SetAt(decayL ,kCutDecLgthMin) ;
416 fCutValues->SetAt(decayL ,kCutDecLgthMax) ;
417 fCutValues->SetAt(decayRxy,kCutDecRxyMin) ;
418 fCutValues->SetAt(decayRxy,kCutDecRxyMax) ;
419
420 // cut on charge
421 if ( fRequireIsCharged || fRequireIsNeutral ) {
422 if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
423 if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
424 }
425 else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
426
427 // cut on primary/secondary
428 if ( fRequireIsPrimary || fRequireIsSecondary) {
429 if (fRequireIsPrimary && IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
430 if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
431 }
432 else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
433
434 // cut on PDG code
435 if ( fRequirePdgCode ) {
3ff4a092 436 if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
b95f6a36 437 }
438 else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
439
440 // cut on decay channel
441 if ( fDecayChannel ) {
442 Bool_t goodDecay = kTRUE ;
443 Short_t nDaughters = 0 ;
444 if (mcPart->GetDaughter(0) >=0) nDaughters = 1 + mcPart->GetDaughter(1) - mcPart->GetDaughter(0);
445
446 if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
447 if (goodDecay) {
448 // now check if decay channel is respected
449 // first try
450 for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
451 AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)+iDaughter));
452 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
453 }
454 if (!goodDecay) {
455 //second try inverting the order of the daughters
456 goodDecay = kTRUE ;
457 for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
458 AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter));
459 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
460 }
461 }
462 }
463 fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
464 }
465 else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
466
467
468 // now array of cut is build, fill the bitmap consequently
469 Int_t iCutBit = -1;
470 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
471 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
472 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
473 if ( fCutValues->At(++iCutBit) > fProdVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
474 if ( fCutValues->At(++iCutBit) < fProdVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
475 if ( fCutValues->At(++iCutBit) > fProdVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
476 if ( fCutValues->At(++iCutBit) < fProdVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
477 if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
478 if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
479 if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
480 if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
481 if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
482 if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
483 if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
484 if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
485 if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
486 if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
487 if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
488 if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
489 if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
490}
491
c7803356 492
107a3100 493//__________________________________________________________________________________
033789c9 494void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
107a3100 495{
496 //
497 // fill the QA histograms
498 //
499
500 for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++)
c7803356 501 fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
107a3100 502
503 // fill cut statistics and cut correlation histograms with information from the bitmap
504 if (afterCuts) return;
505
506 // Number of single cuts in this class
507 UInt_t ncuts = fBitmap->GetNbits();
508 for(UInt_t bit=0; bit<ncuts;bit++) {
c7803356 509 if (!fBitmap->TestBitNumber(bit)) {
107a3100 510 fhCutStatistics->Fill(bit+1);
511 for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
512 if (!fBitmap->TestBitNumber(bit2))
513 fhCutCorrelation->Fill(bit+1,bit2+1);
514 }
515 }
563113d0 516 }
107a3100 517}
518
519//__________________________________________________________________________________
520void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
521 //
522 // saves the histograms in a TList
523 //
563113d0 524
107a3100 525 DefineHistograms();
563113d0 526
107a3100 527 qaList->Add(fhCutStatistics);
528 qaList->Add(fhCutCorrelation);
529
530 for (Int_t j=0; j<kNStepQA; j++) {
531 for(Int_t i=0; i<kNCuts; i++)
532 qaList->Add(fhQA[i][j]);
533 }
534}
535
536//__________________________________________________________________________________
537void AliCFParticleGenCuts::DefineHistograms() {
538 //
539 // histograms for cut variables, cut statistics and cut correlations
540 //
541 Int_t color = 2;
542
543 // book cut statistics and cut correlation histograms
544 fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
545 fhCutStatistics->SetLineWidth(2);
546 int k = 1;
547 fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge") ; k++;
548 fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec") ; k++;
549 fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG") ; k++;
550 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin") ; k++;
551 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax") ; k++;
552 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin") ; k++;
553 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax") ; k++;
554 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin") ; k++;
555 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++;
556 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin") ; k++;
557 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax") ; k++;
558 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin") ; k++;
559 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax") ; k++;
560 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin") ; k++;
561 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++;
562 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
563 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
564 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin") ; k++;
565 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax") ; k++;
8fe1a43d 566 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
107a3100 567
568
569 fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
570 fhCutCorrelation->SetLineWidth(2);
571 for (k=1; k<=kNCuts; k++) {
572 fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
573 fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
574 }
575
576 Char_t str[256];
577 for (int i=0; i<kNStepQA; i++) {
578 if (i==0) sprintf(str," ");
579 else sprintf(str,"_cut");
8fe1a43d 580 fhQA[kCutCharge] [i] = new TH1F(Form("%s_charge%s" ,GetName(),str),"",2,0,2);
581 fhQA[kCutPrimSec] [i] = new TH1F(Form("%s_primSec%s" ,GetName(),str),"",2,0,2);
582 fhQA[kCutPDGCode] [i] = new TH1F(Form("%s_pdgCode%s" ,GetName(),str),"",2,0,2);
583 fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,0,10);
584 fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,0,10);
585 fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,0,10);
586 fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,0,10);
587 fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,0,10);
588 fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,0,10);
589 fhQA[kCutDecVtxXMin] [i] = new TH1F(Form("%s_decVtxXMin%s" ,GetName(),str),"",100,0,10);
590 fhQA[kCutDecVtxXMax] [i] = new TH1F(Form("%s_decVtxXMax%s" ,GetName(),str),"",100,0,10);
591 fhQA[kCutDecVtxYMin] [i] = new TH1F(Form("%s_decVtxYMin%s" ,GetName(),str),"",100,0,10);
592 fhQA[kCutDecVtxYMax] [i] = new TH1F(Form("%s_decVtxYMax%s" ,GetName(),str),"",100,0,10);
593 fhQA[kCutDecVtxZMin] [i] = new TH1F(Form("%s_decVtxZMin%s" ,GetName(),str),"",100,0,10);
594 fhQA[kCutDecVtxZMax] [i] = new TH1F(Form("%s_decVtxZMax%s" ,GetName(),str),"",100,0,10);
595 fhQA[kCutDecLgthMin] [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
596 fhQA[kCutDecLgthMax] [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
597 fhQA[kCutDecRxyMin] [i] = new TH1F(Form("%s_decRxyMin%s" ,GetName(),str),"",100,0,10);
598 fhQA[kCutDecRxyMax] [i] = new TH1F(Form("%s_decRxyMax%s" ,GetName(),str),"",100,0,10);
599 fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
107a3100 600 }
601 for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
563113d0 602}
107a3100 603
604
563113d0 605//______________________________
b95f6a36 606Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) {
563113d0 607 //
608 //check if particle is charged.
609 //
b95f6a36 610 if (mcPart->Charge()==0) return kFALSE;
563113d0 611 return kTRUE;
612}
613//______________________________
b95f6a36 614Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) {
563113d0 615 //
616 //check if particle is primary (standard definition)
617 //
b95f6a36 618
619 AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
620
621 if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE;
563113d0 622 return kTRUE;
623}
624//______________________________
b95f6a36 625Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) {
626 //
627 //check if particle is primary (standard definition)
628 //
629
630 if (!mcPart->IsPhysicalPrimary()) return kFALSE;
631 return kTRUE;
632}
633//______________________________
634Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
563113d0 635 //
636 //check if a charged particle is primary (standard definition)
637 //
b95f6a36 638
639 if (!fIsAODMC) {
640 if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
641 }
642 else {
643 if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
644 }
563113d0 645 return kTRUE;
646}
647//______________________________
648Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
649 //
650 //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
3ff4a092 651 //absolute value.
563113d0 652 //
653 TParticle* part = mcPart->Particle();
654 Int_t pdgCode = part->GetPdgCode();
b95f6a36 655
656 if (abs) {
657 pdgCode = TMath::Abs(pdgCode);
658 pdg = TMath::Abs(pdg);
659 }
660 if (pdgCode != pdg ) return kFALSE;
661 return kTRUE;
662}
663//______________________________
664Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
665 //
666 //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
3ff4a092 667 //absolute value.
b95f6a36 668 //
669 Int_t pdgCode = mcPart->GetPdgCode();
670
fc01457a 671 if (abs) {
672 pdgCode = TMath::Abs(pdgCode);
673 pdg = TMath::Abs(pdg);
674 }
107a3100 675 if (pdgCode != pdg ) return kFALSE;
563113d0 676 return kTRUE;
677}
678//______________________________
107a3100 679void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) {
563113d0 680 //
107a3100 681 // Sets pointer to MC event information (AliMCEvent)
563113d0 682 //
683
107a3100 684 if (!mcEvent) {
685 AliError("Pointer to MC Event is null !");
563113d0 686 return;
687 }
688
107a3100 689 TString className(mcEvent->ClassName());
b95f6a36 690 if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) {
691 AliError("argument must point to an AliMCEvent or an AliAODEvent !");
563113d0 692 return ;
693 }
694
b95f6a36 695 if (fIsAODMC) fMCInfo = dynamic_cast<AliAODEvent*>(mcEvent) ;
696 else fMCInfo = dynamic_cast<AliMCEvent*> (mcEvent) ;
563113d0 697}