AliPhysics  5364b50 (5364b50)
AliConversionSelection.cxx
Go to the documentation of this file.
2 #include "AliAODHeader.h"
3 #include "AliVVZERO.h"
4 #include "AliMultiplicity.h"
5 #include <iostream>
6 
7 
8 // Author Daniel Lohner (Daniel.Lohner@cern.ch)
9 
10 using namespace std;
11 
12 ClassImp(AliConversionSelection)
13 
14 
15 //________________________________________________________________________
17  fInputEvent(NULL),
18  fMCEvent(NULL),
19  fEventCut(evtCut),
20  fConversionCut(convCut),
21  fMesonCut(mesonCut),
22  fESDTrackCuts(NULL),
23  fGoodGammas(NULL),
24  fPi0Candidates(NULL),
25  fBGPi0s(NULL),
26  fRandomizer(NULL),
27  fBGHandler(NULL),
28  fCurrentEventNumber(-1),
29  fIsOwner(kFALSE)
30 {
31  // Default Values
32  fInvMassRange[0]=0.05;
33  fInvMassRange[1]=0.3;
34 }
35 
36 //________________________________________________________________________
38  fInputEvent(NULL),
39  fMCEvent(NULL),
40  fEventCut(NULL),
41  fConversionCut(NULL),
42  fMesonCut(NULL),
43  fESDTrackCuts(NULL),
44  fGoodGammas(NULL),
45  fPi0Candidates(NULL),
46  fBGPi0s(NULL),
47  fRandomizer(NULL),
48  fBGHandler(NULL),
49  fCurrentEventNumber(-1),
50  fIsOwner(kTRUE)
51 {
52  // Default Values
53  fInvMassRange[0]=0.05;
54  fInvMassRange[1]=0.3;
55 
57  fEventCut -> InitializeCutsFromCutString(evtCut.Data());
59  fConversionCut -> InitializeCutsFromCutString(convCut.Data());
61  fMesonCut -> InitializeCutsFromCutString(mesonCut.Data());
62 
63 }
64 
65 
66 //________________________________________________________________________
68  fInputEvent(NULL),
69  fMCEvent(NULL),
70  fEventCut(NULL),
71  fConversionCut(NULL),
72  fMesonCut(NULL),
73  fESDTrackCuts(NULL),
74  fGoodGammas(NULL),
75  fPi0Candidates(NULL),
76  fBGPi0s(NULL),
77  fRandomizer(NULL),
78  fBGHandler(NULL),
80  fIsOwner(kTRUE)
81 {
82  // Copy Constructor
86 
87  fInvMassRange[0]=ref.fInvMassRange[0];
88  fInvMassRange[1]=ref.fInvMassRange[1];
89 }
90 
91 //________________________________________________________________________
93 
94  if(fBGHandler){
95  delete fBGHandler;
96  fBGHandler=NULL;
97  }
98  if(fRandomizer){
99  delete fRandomizer;
100  fRandomizer=NULL;
101  }
102  if(fPi0Candidates){
103  delete fPi0Candidates;
104  fPi0Candidates=NULL;
105  }
106  if(fBGPi0s){
107  delete fBGPi0s;
108  fBGPi0s=NULL;
109  }
110  if(fESDTrackCuts){
111  delete fESDTrackCuts;
112  fESDTrackCuts=NULL;
113  }
114  if(fIsOwner){
115  if(fEventCut){
116  delete fEventCut;
117  fEventCut=NULL;
118  }
119  if(fConversionCut){
120  delete fConversionCut;
121  fConversionCut=NULL;
122  }
123  if(fMesonCut){
124  delete fMesonCut;
125  fMesonCut=NULL;
126  }
127  }
128 }
129 
130 //________________________________________________________________________
131 Bool_t AliConversionSelection::ProcessEvent(TClonesArray *photons,AliVEvent *inputEvent,AliMCEvent *mcEvent){
132  fInputEvent=inputEvent;
133  fMCEvent=mcEvent;
134 
135  // Protection
136  Int_t eventnumber=GetEventNumber(inputEvent);
137  if(eventnumber==fCurrentEventNumber){
138  AliWarning("Event already analyzed! Return.");
139  return kFALSE;
140  }
141  else{
142  fCurrentEventNumber=eventnumber;
143  }
144 
145  // Initialize and Reset Arrays
146  if(fGoodGammas == NULL){
147  fGoodGammas=new TObjArray(30);
148  }
149  fGoodGammas->Clear();
150 
151  if(fPi0Candidates == NULL){
152  fPi0Candidates = new TClonesArray("AliAODConversionMother",100);
153  }
154  fPi0Candidates->Delete();
155 
156  if(fBGPi0s == NULL){
157  fBGPi0s = new TClonesArray("AliAODConversionMother",100);
158  }
159  fBGPi0s->Delete();
160 
161 
162  if(!photons||!fInputEvent)return kFALSE;
163 
164  if(!fEventCut->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
165 
166  // Select photons
167  for(Int_t i = 0; i < photons->GetEntriesFast(); i++) {
168  AliAODConversionPhoton* gamma =dynamic_cast<AliAODConversionPhoton*>(photons->At(i));
169  if(!gamma) continue;
170  if(!fConversionCut->PhotonIsSelected(gamma,fInputEvent))continue;
171  fGoodGammas->Add(gamma);
172  }
173 
174  // Do MC Smearing
175  Double_t *fUnsmearedPx=NULL;
176  Double_t *fUnsmearedPy=NULL;
177  Double_t *fUnsmearedPz=NULL;
178  Double_t *fUnsmearedE=NULL;
179 
181  fUnsmearedPx = new Double_t[fGoodGammas->GetEntries()]; // Store unsmeared Momenta
182  fUnsmearedPy = new Double_t[fGoodGammas->GetEntries()];
183  fUnsmearedPz = new Double_t[fGoodGammas->GetEntries()];
184  fUnsmearedE = new Double_t[fGoodGammas->GetEntries()];
185 
186  for(Int_t gamma=0;gamma<fGoodGammas->GetEntries();gamma++){ // Smear the AODPhotons in MC
187  fUnsmearedPx[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Px();
188  fUnsmearedPy[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Py();
189  fUnsmearedPz[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Pz();
190  fUnsmearedE[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->E();
191  fMesonCut->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(gamma)));
192  }
193  }
194 
195  // Reconstruct Pi0 and BG
199 
200  // Undo MC Smearing
202  for(Int_t gamma=0;gamma<fGoodGammas->GetEntries();gamma++){ // Smear the AODPhotons in MC
203  ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPx(fUnsmearedPx[gamma]); // Reset Unsmeared Momenta
204  ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPy(fUnsmearedPy[gamma]);
205  ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPz(fUnsmearedPz[gamma]);
206  ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetE(fUnsmearedE[gamma]);
207  }
208  delete[] fUnsmearedPx; fUnsmearedPx = 0x0;
209  delete[] fUnsmearedPy; fUnsmearedPy = 0x0;
210  delete[] fUnsmearedPz; fUnsmearedPz = 0x0;
211  delete[] fUnsmearedE; fUnsmearedE = 0x0;
212  }
213  return kTRUE;
214 }
215 
216 //________________________________________________________________________
218 
219  if(index>=0&&index<GetNumberOfPi0s()){
220  return dynamic_cast<AliAODConversionMother*>(fPi0Candidates->At(index));
221  }
222  return NULL;
223 }
224 
225 //________________________________________________________________________
227 
228  if(index>=0&&index<GetNumberOfBGs()){
229  return dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(index));
230  }
231  return NULL;
232 }
233 
234 //________________________________________________________________________
236 
237  if(index>=0&&index<GetNumberOfPhotons()){
238  return dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(index));
239  }
240  return NULL;
241 }
242 
243 //________________________________________________________________________
245  // Conversion Gammas
246  if(fGoodGammas->GetEntriesFast()>1){
247  for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntriesFast()-1;firstGammaIndex++){
248  AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(firstGammaIndex));
249  if (gamma0==NULL) continue;
250  // Combine Photons
251 
252  for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntriesFast();secondGammaIndex++){
253 
254  AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));
255  if (gamma1==NULL) continue;
256  //Check for same Electron ID
257  if(gamma0->GetTrackLabelPositive()==gamma1->GetTrackLabelPositive()||gamma0->GetTrackLabelNegative()==gamma1->GetTrackLabelNegative()
258  ||gamma0->GetTrackLabelNegative()==gamma1->GetTrackLabelPositive()||gamma0->GetTrackLabelPositive()==gamma1->GetTrackLabelNegative())continue;
259 
260  AliAODConversionMother pi0cand(gamma0,gamma1);
261  pi0cand.SetLabels(firstGammaIndex,secondGammaIndex);
262 
263  // Set MC Label
264 
265  if(fMCEvent){
266 
267  TParticle *mcgam0=gamma0->GetMCParticle(fMCEvent);
268  TParticle *mcgam1=gamma1->GetMCParticle(fMCEvent);
269 
270  if(mcgam0&&mcgam1){
271  // Have same Mother?
272 
273  if(mcgam0->GetMother(0)==mcgam1->GetMother(0)){
274 
275  pi0cand.SetMCLabel(mcgam0->GetMother(0));
276  }
277  }
278  }
279 
280  if((fMesonCut->MesonIsSelected(&pi0cand))){
281  if(MesonInMassWindow(&pi0cand)){
282 
283  // Add Pi0 to Stack
284  new((*fPi0Candidates)[fPi0Candidates->GetEntriesFast()]) AliAODConversionMother(pi0cand);
285  }
286  }
287  }
288  }
289  }
290 }
291 
292 //________________________________________________________________________
294 {
295  if (pi0cand->M() > fInvMassRange[0] && pi0cand->M() < fInvMassRange[1] ){
296  return kTRUE;
297  }
298  return kFALSE;
299 }
300 
301 //________________________________________________________________________
303 
304  if(!fRandomizer){
305  fRandomizer=new TRandom3();
306  fRandomizer->SetSeed(0);
307  }
308  Double_t nRadiansPM = nDegreesPMBackground*TMath::Pi()/180;
309  Double_t rotationValue = fRandomizer->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
310  gamma->RotateZ(rotationValue);
311 }
312 
313 //________________________________________________________________________
314 
316 
317  //Rotation Method
319  // Correct for the number of rotations
320  // BG is for rotation the same, except for factor NRotations
322  for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntriesFast();firstGammaIndex++){
323  AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(firstGammaIndex));
324  if (gamma0 ==NULL) continue;
325  for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntriesFast();secondGammaIndex++){
326  AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));
327  if (gamma1==NULL) continue;
328  if(!fConversionCut->PhotonIsSelected(gamma1,fInputEvent))continue;
329  for(Int_t nRandom=0;nRandom<fMesonCut->GetNumberOfBGEvents();nRandom++){
331  AliAODConversionMother BGcandidate(gamma0,gamma1);
332  if(fMesonCut->MesonIsSelected(&BGcandidate,kFALSE)){
333  if(MesonInMassWindow(&BGcandidate)){
334  new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(BGcandidate);
335  dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
336  }
337  }
338  }
339  }
340  }
341  } else {
342  // Do Event Mixing
343  if(fBGHandler==NULL){
345  }
346 
347  for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler->GetNBGEvents(fGoodGammas,fInputEvent);nEventsInBG++){
349  if(previousEventGammas){
350  // test weighted background
351  Double_t weight=1.0;
352  // Correct for the number of eventmixing:
353  // N gammas -> (N-1) + (N-2) +(N-3) ...+ (N-(N-1)) using sum formula sum(i)=N*(N-1)/2 -> N*(N-1)/2
354  // real combinations (since you cannot combine a photon with its own)
355  // but BG leads to N_{a}*N_{b} combinations
356  weight*=0.5*(Double_t(fGoodGammas->GetEntriesFast()-1))/Double_t(previousEventGammas->size());
357  for(Int_t iCurrent=0;iCurrent<fGoodGammas->GetEntriesFast();iCurrent++){
358  AliAODConversionPhoton *gamma0 = (AliAODConversionPhoton*)(fGoodGammas->At(iCurrent));
359  for(UInt_t iPrevious=0;iPrevious<previousEventGammas->size();iPrevious++){
360  AliAODConversionPhoton *gamma1 = (AliAODConversionPhoton*)(previousEventGammas->at(iPrevious));
361  AliAODConversionMother BGcandidate(gamma0,gamma1);
362  if(fMesonCut->MesonIsSelected(&BGcandidate,kFALSE)){
363  if(MesonInMassWindow(&BGcandidate)){
364  new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(BGcandidate);
365  dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
366  }
367  }
368  }
369  }
370  }
371  }
372  }
373 }
374 //________________________________________________________________________
376 
377  switch(fEventCut->GetMultiplicityMethod()){
378  case 0:
379  return Double_t(GetNumberOfPhotons());
380  case 1:
381  return Double_t(GetNumberOfChargedTracks(inputEvent));
382  case 2:
383  return GetVZEROMult(inputEvent);
384  case 3:
385  return GetSPDMult(inputEvent);
386  case 9:
387  return 1; // if mult is used as a weight, this number can be used to switch off weighting
388  default:
389  return 0;
390  }
391 }
392 
393 //________________________________________________________________________
395 
396  Int_t ntracks = 0;
397 
398  AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
399  if(esdEvent) {
400  if(!fESDTrackCuts){
401  fESDTrackCuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
402  fESDTrackCuts->SetMaxDCAToVertexZ(2);
404  fESDTrackCuts->SetEtaRange(-etamax, etamax);
405  fESDTrackCuts->SetPtRange(0.15);
406  }
407  for(Int_t iTracks = 0; iTracks < inputEvent->GetNumberOfTracks(); iTracks++){
408  AliESDtrack* currentTrack = esdEvent->GetTrack(iTracks);
409  if(!currentTrack) continue;
410  if(fESDTrackCuts->AcceptTrack(currentTrack))ntracks++;
411  }
412  } else {
413  for(Int_t ii=0; ii<inputEvent->GetNumberOfTracks(); ii++) {
414  AliVTrack * track = dynamic_cast<AliVTrack*>(inputEvent->GetTrack(ii));
415  if (track==NULL) continue;
416  if(TMath::Abs(track->Eta())>fConversionCut->GetEtaCut())continue;
417  if(track)ntracks++;
418  }
419  }
420  return ntracks;
421 }
422 
423 //________________________________________________________________________
425 
426  AliVVZERO *vzero=inputEvent->GetVZEROData();
427  Double_t multV0A=vzero->GetMTotV0A();
428  Double_t multV0C=vzero->GetMTotV0C();
429  Double_t mult=multV0A+multV0C;
430 
431  return mult;
432 }
433 
434 //________________________________________________________________________
436 
437  AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
438  if(esdEvent) {
439  const AliMultiplicity *esdmult=esdEvent->GetMultiplicity();
440  return esdmult->GetNumberOfITSClusters(1);
441  } else {
442  // AOD implementation
443  AliAODHeader *header=(AliAODHeader*)inputEvent->GetHeader();
444  return header->GetNumberOfITSClusters(1);
445  }
446  return 0;
447 }
448 
449 //________________________________________________________________________
451 
452  AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(inputEvent);
453  if(esdEvent) {
454  return esdEvent->GetEventNumberInFile();
455  } else{
456  AliAODHeader *header=(AliAODHeader*)inputEvent->GetHeader();
457  return header->GetEventNumberESDFile();
458  }
459  return 0;
460 }
461 
462 //________________________________________________________________________
464  TString a= Form("%s_%s_%s",fEventCut->GetCutNumber().Data(), fConversionCut->GetCutNumber().Data(),fMesonCut->GetCutNumber().Data());
465  return a;
466 }
467 
TParticle * GetMCParticle(AliMCEvent *mcEvent)
Int_t GetNumberOfChargedTracks(AliVEvent *inputEvent)
AliAODConversionMother * GetPi0(Int_t index)
void SetLabels(Int_t label1, Int_t label2, Int_t label3=0)
double Double_t
Definition: External.C:58
void SetWeight(Double_t weight)
AliGammaConversionPhotonVector * GetBGGoodGammas(TObjArray *const eventGammas, AliVEvent *fInputEvent, Int_t event)
AliAODConversionMother * GetBG(Int_t index)
Double_t GetVZEROMult(AliVEvent *inputEvent)
AliConversionPhotonCuts * fConversionCut
int Int_t
Definition: External.C:63
AliConversionAODBGHandlerRP * fBGHandler
unsigned int UInt_t
Definition: External.C:33
Bool_t MesonInMassWindow(AliAODConversionMother *pi0cand)
void SmearParticle(AliAODConversionPhoton *photon)
Double_t GetMultiplicity(AliVEvent *inputEvent)
AliAODConversionPhoton * GetPhoton(Int_t index)
Class handling all kinds of selection cuts for Gamma Conversion analysis.
AliConvEventCuts * fEventCut
Bool_t EventIsSelected(AliVEvent *fInputEvent, AliMCEvent *fMCEvent)
Bool_t ProcessEvent(TClonesArray *photons, AliVEvent *inputEvent, AliMCEvent *mcEvent)
Class handling all kinds of selection cuts for Gamma Conversion analysis.
Bool_t MesonIsSelected(AliAODConversionMother *pi0, Bool_t IsSignal=kTRUE, Double_t fRapidityShift=0., Int_t leadingCellID1=0, Int_t leadingCellID2=0)
Int_t GetMultiplicityMethod()
const Double_t etamax
void AddEvent(TObjArray *const eventGammas, AliVEvent *fInputEvent)
void RotateParticle(AliAODConversionPhoton *gamma, Int_t nDegreesPMBackground)
AliESDtrackCuts * fESDTrackCuts
Class handling all kinds of selection cuts for Gamma Conversion analysis.
Double_t GetSPDMult(AliVEvent *inputEvent)
bool Bool_t
Definition: External.C:53
AliConversionSelection(AliConvEventCuts *evtCut=NULL, AliConversionPhotonCuts *convCut=NULL, AliConversionMesonCuts *mesonCut=NULL)
Int_t GetEventNumber(AliVEvent *inputEvent)
vector< AliAODConversionPhoton * > AliGammaConversionPhotonVector
Bool_t PhotonIsSelected(AliConversionPhotonBase *photon, AliVEvent *event)
AliConversionMesonCuts * fMesonCut