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