1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 ////////////////////////////////////////////////////////////////////////////
26 #include "AliCFParticleGenCuts.h"
27 #include "TParticle.h"
28 #include "TParticlePDG.h"
29 #include "AliMCEvent.h"
37 #include "TDecayChannel.h"
38 #include "AliAODMCParticle.h"
39 #include "AliAODEvent.h"
41 ClassImp(AliCFParticleGenCuts)
43 //______________________________
44 AliCFParticleGenCuts::AliCFParticleGenCuts() :
51 fRequireIsSecondary(0),
53 fRequireAbsolutePdg(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),
69 fDecayLengthMax(1.e+09),
74 fhCutCorrelation(0x0),
75 fCutValues(new TArrayF(kNCuts)),
81 for (int i=0; i<kNCuts; i++)
82 for (int j=0; j<kNStepQA; j++)
85 for (int j=0; j<kNStepQA; j++)
89 //______________________________
90 AliCFParticleGenCuts::AliCFParticleGenCuts(const Char_t* name, const Char_t* title) :
91 AliCFCutBase(name,title),
97 fRequireIsSecondary(0),
99 fRequireAbsolutePdg(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),
114 fDecayLengthMin(-1.),
115 fDecayLengthMax(1.e+09),
117 fDecayRxyMax(1.e+09),
119 fhCutStatistics(0x0),
120 fhCutCorrelation(0x0),
121 fCutValues(new TArrayF(kNCuts)),
122 fBitmap(new TBits(0))
127 for (int i=0; i<kNCuts; i++)
128 for (int j=0; j<kNStepQA; j++)
131 for (int j=0; j<kNStepQA; j++)
135 //______________________________
136 AliCFParticleGenCuts::AliCFParticleGenCuts(const AliCFParticleGenCuts& c) :
138 fIsAODMC(c.fIsAODMC),
140 fRequireIsCharged(c.fRequireIsCharged),
141 fRequireIsNeutral(c.fRequireIsNeutral),
142 fRequireIsPrimary(c.fRequireIsPrimary),
143 fRequireIsSecondary(c.fRequireIsSecondary),
144 fRequirePdgCode(c.fRequirePdgCode),
145 fRequireAbsolutePdg(c.fRequireAbsolutePdg),
146 fProdVtxRange2D(c.fProdVtxRange2D),
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),
163 fDecayRxyMax(c.fDecayLengthMin),
164 fDecayChannel(c.fDecayChannel),
165 fhCutStatistics(new TH1F(*c.fhCutStatistics)),
166 fhCutCorrelation(new TH2F(*c.fhCutCorrelation)),
167 fCutValues(new TArrayF(*c.fCutValues)),
168 fBitmap(new TBits(*c.fBitmap))
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();
177 for (int j=0; j<kNStepQA; j++)
178 fhProdVtxXY[j]=(TH2F*)c.fhProdVtxXY[j]->Clone();
181 //______________________________
182 AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts& c)
185 // Assignment operator
188 AliCFCutBase::operator=(c) ;
191 fRequireIsCharged=c.fRequireIsCharged;
192 fRequireIsNeutral=c.fRequireIsNeutral;
193 fRequireIsPrimary=c.fRequireIsPrimary;
194 fRequireIsSecondary=c.fRequireIsSecondary;
195 fRequirePdgCode=c.fRequirePdgCode;
196 fRequireAbsolutePdg=c.fRequireAbsolutePdg;
197 fProdVtxRange2D=c.fProdVtxRange2D;
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;
215 fDecayChannel=c.fDecayChannel;
216 fCutValues=new TArrayF(*c.fCutValues);
217 fBitmap=new TBits(*c.fBitmap);
219 if (fhCutStatistics) fhCutStatistics =new TH1F(*c.fhCutStatistics) ;
220 if (fhCutCorrelation) fhCutCorrelation=new TH2F(*c.fhCutCorrelation);
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();
226 for (int j=0; j<kNStepQA; j++)
227 fhProdVtxXY[j]=(TH2F*)c.fhProdVtxXY[j]->Clone();
232 //______________________________
233 Bool_t AliCFParticleGenCuts::IsSelected(TObject* obj) {
235 // check if selections on 'obj' are passed
236 // 'obj' must be an AliMCParticle
239 if (!obj) return kFALSE ;
241 if (!fIsAODMC) SelectionBitMap((AliMCParticle*) obj);
242 else SelectionBitMap((AliAODMCParticle*)obj);
244 if (fIsQAOn) FillHistograms(obj,0);
246 for (UInt_t icut=0; icut<fBitmap->GetNbits();icut++)
247 if (!fBitmap->TestBitNumber(icut)) return kFALSE ;
249 if (fIsQAOn) FillHistograms(obj,1);
253 //__________________________________________________________________________________
254 void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart)
257 // test if the track passes the single cuts
258 // and store the information in a bitmap
261 for (UInt_t i=0; i<kNCuts; i++) {
262 fBitmap->SetBitNumber(i,kFALSE);
263 fCutValues->SetAt((Double32_t)0,i) ;
266 // fill the cut array
267 Double32_t partVx=(Double32_t)mcPart->Xv();
268 Double32_t partVy=(Double32_t)mcPart->Yv();
269 Double32_t partVz=(Double32_t)mcPart->Zv();
271 // calculate the production vertex ellipse
272 Double32_t prodVtxXYmin = 0.;
273 if (fProdVtxXMin>0 && fProdVtxYMin>0)
274 prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
275 Double32_t prodVtxXYmax = 0.;
276 if(fProdVtxXMax>0 && fProdVtxYMax>0)
277 prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
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.;
285 TParticle* part = mcPart->Particle();
286 AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
287 TParticle* daughter=0x0;
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) ) ;
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) ;
317 if ( fRequireIsCharged || fRequireIsNeutral ) {
318 if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
319 if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
321 else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
323 // cut on primary/secondary
324 if ( fRequireIsPrimary || fRequireIsSecondary) {
325 if (fRequireIsPrimary && IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
326 if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
328 else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
331 if ( fRequirePdgCode ) {
332 if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
334 else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
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 ;
341 //now number of daughters is OK
343 // now check if decay channel is respected
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;}
350 //second try inverting the order of the daughters
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;}
357 if (!goodDecay && fRequireAbsolutePdg) {
358 //now tries inverting the sign of the daughters in case the anti-particle is also looked at
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;}
366 //fourth try inverting the order of the daughters
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;}
373 } //end check anti-particle
374 } //end # daughters OK
375 fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
376 } //end require decay channel
377 else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
380 // now array of cut is build, fill the bitmap consequently
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);
387 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
388 || ( fProdVtxRange2D && (fProdVtxXMin>0 && fProdVtxYMin>0) && prodVtxXYmin >= 1)
389 || ( fProdVtxRange2D && (fProdVtxXMin<=0 || fProdVtxYMin<=0) ) )
390 fBitmap->SetBitNumber(iCutBit,kTRUE);
393 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
394 || ( fProdVtxRange2D && (fProdVtxXMax>0 && fProdVtxYMax>0) && prodVtxXYmax <= 1)
395 || ( fProdVtxRange2D && (fProdVtxXMax<=0 || fProdVtxYMax<=0) ) )
396 fBitmap->SetBitNumber(iCutBit,kTRUE);
399 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
400 || ( fProdVtxRange2D && (fProdVtxXMin>0 && fProdVtxYMin>0) && prodVtxXYmin >= 1)
401 || ( fProdVtxRange2D && (fProdVtxXMin<=0 || fProdVtxYMin<=0) ) )
402 fBitmap->SetBitNumber(iCutBit,kTRUE);
405 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
406 || ( fProdVtxRange2D && (fProdVtxXMax>0 && fProdVtxYMax>0) && prodVtxXYmax <= 1)
407 || ( fProdVtxRange2D && (fProdVtxXMax<=0 || fProdVtxYMax<=0) ) )
408 fBitmap->SetBitNumber(iCutBit,kTRUE);
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);
422 if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
425 //__________________________________________________________________________________
426 void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
429 // test if the track passes the single cuts
430 // and store the information in a bitmap
433 for (UInt_t i=0; i<kNCuts; i++) {
434 fBitmap->SetBitNumber(i,kFALSE);
435 fCutValues->SetAt((Double32_t)0,i) ;
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();
443 // calculate the production vertex ellipse
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);
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.;
457 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fMCInfo);
458 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
459 AliAODMCParticle* daughter = 0x0 ;
461 if (mcPart->GetDaughter(0)>0) {
462 daughter = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
464 decayVx=(Double32_t)daughter->Xv();
465 decayVy=(Double32_t)daughter->Yv();
466 decayVz=(Double32_t)daughter->Zv();
467 decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) +
468 TMath::Power(partVy-decayVy,2) +
469 TMath::Power(partVz-decayVz,2) ) ;
470 decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
474 fCutValues->SetAt(partVx ,kCutProdVtxXMin) ;
475 fCutValues->SetAt(partVy ,kCutProdVtxYMin) ;
476 fCutValues->SetAt(partVz ,kCutProdVtxZMin) ;
477 fCutValues->SetAt(partVx ,kCutProdVtxXMax) ;
478 fCutValues->SetAt(partVy ,kCutProdVtxYMax) ;
479 fCutValues->SetAt(partVz ,kCutProdVtxZMax) ;
480 fCutValues->SetAt(decayVx ,kCutDecVtxXMin) ;
481 fCutValues->SetAt(decayVy ,kCutDecVtxYMin) ;
482 fCutValues->SetAt(decayVz ,kCutDecVtxZMin) ;
483 fCutValues->SetAt(decayVx ,kCutDecVtxXMax) ;
484 fCutValues->SetAt(decayVy ,kCutDecVtxYMax) ;
485 fCutValues->SetAt(decayVz ,kCutDecVtxZMax) ;
486 fCutValues->SetAt(decayL ,kCutDecLgthMin) ;
487 fCutValues->SetAt(decayL ,kCutDecLgthMax) ;
488 fCutValues->SetAt(decayRxy,kCutDecRxyMin) ;
489 fCutValues->SetAt(decayRxy,kCutDecRxyMax) ;
492 if ( fRequireIsCharged || fRequireIsNeutral ) {
493 if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
494 if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
496 else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
498 // cut on primary/secondary
499 if ( fRequireIsPrimary || fRequireIsSecondary) {
500 if (fRequireIsPrimary && IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
501 if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
503 else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
506 if ( fRequirePdgCode ) {
507 if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
509 else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
511 // cut on decay channel
512 if ( fDecayChannel ) {
513 Bool_t goodDecay = kTRUE ;
514 Short_t nDaughters = 0 ;
515 if (mcPart->GetDaughter(0) >=0) nDaughters = 1 + mcPart->GetDaughter(1) - mcPart->GetDaughter(0);
517 if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
519 // now check if decay channel is respected
521 for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
522 AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)+iDaughter));
523 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
526 //second try inverting the order of the daughters
528 for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
529 AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter));
530 if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
534 fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
536 else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
539 // now array of cut is build, fill the bitmap consequently
541 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
542 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
543 if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
546 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
547 || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
550 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
551 || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
554 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
555 || ( fProdVtxRange2D && prodVtxXYmin >= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
558 if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
559 || ( fProdVtxRange2D && prodVtxXYmax <= 1)) fBitmap->SetBitNumber(iCutBit,kTRUE);
561 if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
562 if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
563 if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
564 if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
565 if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
566 if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
567 if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
568 if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
569 if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
570 if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
571 if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
572 if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
573 if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE);
577 //__________________________________________________________________________________
578 void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
581 // fill the QA histograms
584 for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++)
585 fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
587 fhProdVtxXY[afterCuts]->Fill(fCutValues->At(4),fCutValues->At(5));
589 // fill cut statistics and cut correlation histograms with information from the bitmap
590 if (afterCuts) return;
592 // Number of single cuts in this class
593 UInt_t ncuts = fBitmap->GetNbits();
594 for(UInt_t bit=0; bit<ncuts;bit++) {
595 if (!fBitmap->TestBitNumber(bit)) {
596 fhCutStatistics->Fill(bit+1);
597 for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
598 if (!fBitmap->TestBitNumber(bit2))
599 fhCutCorrelation->Fill(bit+1,bit2+1);
605 //__________________________________________________________________________________
606 void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
608 // saves the histograms in a TList
613 qaList->Add(fhCutStatistics);
614 qaList->Add(fhCutCorrelation);
616 for (Int_t j=0; j<kNStepQA; j++) {
617 qaList->Add(fhProdVtxXY[j]);
618 for(Int_t i=0; i<kNCuts; i++)
619 qaList->Add(fhQA[i][j]);
623 //__________________________________________________________________________________
624 void AliCFParticleGenCuts::DefineHistograms() {
626 // histograms for cut variables, cut statistics and cut correlations
630 // book cut statistics and cut correlation histograms
631 fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
632 fhCutStatistics->SetLineWidth(2);
634 fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge") ; k++;
635 fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec") ; k++;
636 fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG") ; k++;
637 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin") ; k++;
638 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax") ; k++;
639 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin") ; k++;
640 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax") ; k++;
641 fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin") ; k++;
642 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++;
643 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin") ; k++;
644 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax") ; k++;
645 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin") ; k++;
646 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax") ; k++;
647 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin") ; k++;
648 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++;
649 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
650 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
651 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin") ; k++;
652 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax") ; k++;
653 fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
656 fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
657 fhCutCorrelation->SetLineWidth(2);
658 for (k=1; k<=kNCuts; k++) {
659 fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
660 fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
664 for (int i=0; i<kNStepQA; i++) {
665 if (i==0) sprintf(str," ");
666 else sprintf(str,"_cut");
667 fhQA[kCutCharge] [i] = new TH1F(Form("%s_charge%s" ,GetName(),str),"",2,0,2);
668 fhQA[kCutPrimSec] [i] = new TH1F(Form("%s_primSec%s" ,GetName(),str),"",2,0,2);
669 fhQA[kCutPDGCode] [i] = new TH1F(Form("%s_pdgCode%s" ,GetName(),str),"",2,0,2);
670 fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,-10,10);
671 fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,-10,10);
672 fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,-10,10);
673 fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,-10,10);
674 fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,-10,10);
675 fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,-10,10);
676 fhQA[kCutDecVtxXMin] [i] = new TH1F(Form("%s_decVtxXMin%s" ,GetName(),str),"",100,0,10);
677 fhQA[kCutDecVtxXMax] [i] = new TH1F(Form("%s_decVtxXMax%s" ,GetName(),str),"",100,0,10);
678 fhQA[kCutDecVtxYMin] [i] = new TH1F(Form("%s_decVtxYMin%s" ,GetName(),str),"",100,0,10);
679 fhQA[kCutDecVtxYMax] [i] = new TH1F(Form("%s_decVtxYMax%s" ,GetName(),str),"",100,0,10);
680 fhQA[kCutDecVtxZMin] [i] = new TH1F(Form("%s_decVtxZMin%s" ,GetName(),str),"",100,0,10);
681 fhQA[kCutDecVtxZMax] [i] = new TH1F(Form("%s_decVtxZMax%s" ,GetName(),str),"",100,0,10);
682 fhQA[kCutDecLgthMin] [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
683 fhQA[kCutDecLgthMax] [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
684 fhQA[kCutDecRxyMin] [i] = new TH1F(Form("%s_decRxyMin%s" ,GetName(),str),"",100,0,10);
685 fhQA[kCutDecRxyMax] [i] = new TH1F(Form("%s_decRxyMax%s" ,GetName(),str),"",100,0,10);
686 fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
687 fhProdVtxXY [i] = new TH2F(Form("%s_prodVtxXY%s" ,GetName(),str),"",100,0,10,100,0,10);
688 fhProdVtxXY [i] ->GetXaxis()->SetTitle("x_{production vertex}");
689 fhProdVtxXY [i] ->GetYaxis()->SetTitle("y_{production vertex}");
690 fhQA[kCutProdVtxXMax] [i] ->GetXaxis()->SetTitle("x_{production vertex}");
691 fhQA[kCutProdVtxYMax] [i] ->GetXaxis()->SetTitle("y_{production vertex}");
693 for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
697 //______________________________
698 Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) {
700 //check if particle is charged.
702 if (mcPart->Charge()==0) return kFALSE;
705 //______________________________
706 Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) {
708 //check if particle is primary (standard definition)
711 AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
713 if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE;
716 //______________________________
717 Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) {
719 //check if particle is primary (standard definition)
722 if (!mcPart->IsPhysicalPrimary()) return kFALSE;
725 //______________________________
726 Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
728 //check if a charged particle is primary (standard definition)
732 if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
735 if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
739 //______________________________
740 Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
742 //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
745 TParticle* part = mcPart->Particle();
746 Int_t pdgCode = part->GetPdgCode();
749 pdgCode = TMath::Abs(pdgCode);
750 pdg = TMath::Abs(pdg);
752 if (pdgCode != pdg ) return kFALSE;
755 //______________________________
756 Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
758 //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
761 Int_t pdgCode = mcPart->GetPdgCode();
764 pdgCode = TMath::Abs(pdgCode);
765 pdg = TMath::Abs(pdg);
767 if (pdgCode != pdg ) return kFALSE;
770 //______________________________
771 void AliCFParticleGenCuts::SetMCEventInfo(const TObject* mcEvent) {
773 // Sets pointer to MC event information (AliMCEvent)
777 AliError("Pointer to MC Event is null !");
781 TString className(mcEvent->ClassName());
782 if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) {
783 AliError("argument must point to an AliMCEvent or an AliAODEvent !");
787 fMCInfo = (AliVEvent*)mcEvent ;