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