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