TENDER becomes Tender, removing .so
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliConversionSelection.cxx
1 #include "AliConversionSelection.h"
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 //________________________________________________________________________
16 AliConversionSelection::AliConversionSelection(AliConvEventCuts *evtCut, AliConversionPhotonCuts *convCut, AliConversionMesonCuts *mesonCut) : TObject(),
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 //________________________________________________________________________
37 AliConversionSelection::AliConversionSelection(TString evtCut, TString convCut, TString mesonCut) : TObject(),
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
56         fEventCut = new AliConvEventCuts();
57         fEventCut -> InitializeCutsFromCutString(evtCut.Data());
58         fConversionCut = new AliConversionPhotonCuts();
59         fConversionCut -> InitializeCutsFromCutString(convCut.Data());
60         fMesonCut = new AliConversionMesonCuts();
61         fMesonCut -> InitializeCutsFromCutString(mesonCut.Data());
62
63 }
64
65
66 //________________________________________________________________________
67 AliConversionSelection::AliConversionSelection(const AliConversionSelection &ref) : TObject(ref),
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),
79         fCurrentEventNumber(-1),
80         fIsOwner(kTRUE)
81 {
82         // Copy Constructor
83         fEventCut=new AliConvEventCuts(*ref.fEventCut);
84         fConversionCut=new AliConversionPhotonCuts(*ref.fConversionCut);
85         fMesonCut=new AliConversionMesonCuts(*ref.fMesonCut);
86
87         fInvMassRange[0]=ref.fInvMassRange[0];
88         fInvMassRange[1]=ref.fInvMassRange[1];
89 }
90
91 //________________________________________________________________________
92 AliConversionSelection::~AliConversionSelection(){
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
180         if(fMesonCut->UseMCPSmearing() && fMCEvent){
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
196         CalculatePi0Candidates();
197         CalculateBackground();
198         if(fBGHandler)fBGHandler->AddEvent(fGoodGammas,fInputEvent);
199
200         // Undo MC Smearing
201         if(fMesonCut->UseMCPSmearing() && fMCEvent){
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 //________________________________________________________________________
217 AliAODConversionMother* AliConversionSelection::GetPi0(Int_t index){
218
219         if(index>=0&&index<GetNumberOfPi0s()){
220                 return dynamic_cast<AliAODConversionMother*>(fPi0Candidates->At(index));
221         }
222         return NULL;
223 }
224
225 //________________________________________________________________________
226 AliAODConversionMother* AliConversionSelection::GetBG(Int_t index){
227
228         if(index>=0&&index<GetNumberOfBGs()){
229                 return dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(index));
230         }
231         return NULL;
232 }
233
234 //________________________________________________________________________
235 AliAODConversionPhoton* AliConversionSelection::GetPhoton(Int_t index){
236
237         if(index>=0&&index<GetNumberOfPhotons()){
238                 return dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(index));
239         }
240         return NULL;
241 }
242
243 //________________________________________________________________________
244 void AliConversionSelection::CalculatePi0Candidates(){
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->Stack());
268                                         TParticle *mcgam1=gamma1->GetMCParticle(fMCEvent->Stack());
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 //________________________________________________________________________
293 Bool_t AliConversionSelection::MesonInMassWindow(AliAODConversionMother *pi0cand)
294 {
295         if (pi0cand->M() > fInvMassRange[0] && pi0cand->M() < fInvMassRange[1] ){
296                 return kTRUE;
297         }
298         return kFALSE;
299 }
300
301 //________________________________________________________________________
302 void AliConversionSelection::RotateParticle(AliAODConversionPhoton *gamma,Int_t nDegreesPMBackground){
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
315 void AliConversionSelection::CalculateBackground(){
316
317         //Rotation Method
318         if(fMesonCut->UseRotationMethod()){
319                 // Correct for the number of rotations
320                 // BG is for rotation the same, except for factor NRotations
321                 Double_t weight=1./Double_t(fMesonCut->GetNumberOfBGEvents());
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++){
330                                         RotateParticle(gamma1,fMesonCut->NDegreesRotation());
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){
344                         fBGHandler=new AliConversionAODBGHandlerRP(fEventCut->IsHeavyIon(),fMesonCut->UseTrackMultiplicity());
345                 }
346
347                 for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler->GetNBGEvents(fGoodGammas,fInputEvent);nEventsInBG++){
348                         AliGammaConversionPhotonVector *previousEventGammas = fBGHandler->GetBGGoodGammas(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 //________________________________________________________________________
375 Double_t AliConversionSelection::GetMultiplicity(AliVEvent *inputEvent){
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 //________________________________________________________________________
394 Int_t AliConversionSelection::GetNumberOfChargedTracks(AliVEvent *inputEvent){
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);
403                         Double_t etamax=fConversionCut->GetEtaCut();
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 //________________________________________________________________________
424 Double_t AliConversionSelection::GetVZEROMult(AliVEvent *inputEvent){
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 //________________________________________________________________________
435 Double_t AliConversionSelection::GetSPDMult(AliVEvent *inputEvent){
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 //________________________________________________________________________
450 Int_t AliConversionSelection::GetEventNumber(AliVEvent *inputEvent){
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 //________________________________________________________________________
463 TString AliConversionSelection::GetCutString(){
464         TString a= Form("%s_%s_%s",fEventCut->GetCutNumber().Data(), fConversionCut->GetCutNumber().Data(),fMesonCut->GetCutNumber().Data());
465         return a;
466 }
467