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