PWGPP/AliFilteredTreeEventCuts.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 #include <iostream>
  17 #include <TList.h>
  18
  19 #include "AliLog.h"
  20 #include "AliESDEvent.h"
  21 #include "AliESDVertex.h"
  22 #include "AliMCEvent.h"
  23 #include "AliHeader.h"
  24 #include "AliGenEventHeader.h"
  25 #include "AliGenPythiaEventHeader.h"
  26 #include "AliGenCocktailEventHeader.h"
  27 #include "AliGenDPMjetEventHeader.h"
  28 #include "AliStack.h"
  29
  30
  31 #include "AliFilteredTreeEventCuts.h"
  32
  33 using namespace std;
  34
  35 ClassImp(AliFilteredTreeEventCuts)
  36
  37 Int_t AliFilteredTreeEventCuts::fgLastProcessType = -1;
  38
  39 //_____________________________________________________________________________
  40 AliFilteredTreeEventCuts::AliFilteredTreeEventCuts(const Char_t* name,const Char_t *title) :
  41 AliAnalysisCuts(name, title)
  42 , fTriggerRequired(kTRUE)
  43 , fRecVertexRequired(kTRUE)
  44 , fEventProcessType(kInvalidProcess)
  45 , fMinNContributors(0)
  46 , fMaxNContributors(0)
  47 , fMaxR(0)
  48 , fMinZv(0)
  49 , fMaxZv(0)
  50 , fMeanXv(0)
  51 , fMeanYv(0)
  52 , fMeanZv(0)
  53 , fSigmaMeanXv(0)
  54 , fSigmaMeanYv(0)
  55 , fSigmaMeanZv(0)
  56 , fRedoTPCVertex(kTRUE)
  57 , fUseBeamSpotConstraint(kTRUE)
  58 , fEventSelectedRequired(kFALSE)
  59 {
  60   // default constructor
  61
  62   // init data members with defaults
  63   Init();
  64 }
  65
  66 //_____________________________________________________________________________
  67 AliFilteredTreeEventCuts::~AliFilteredTreeEventCuts()
  68 {
  69   // destructor
  70 }
  71
  72 //_____________________________________________________________________________
  73 void AliFilteredTreeEventCuts::Init()
  74 {
  75   // set default values
  76   SetTriggerRequired();
  77   SetRecVertexRequired();
  78   SetEventProcessType();
  79   SetNContributorsRange();
  80   SetMaxR();
  81   SetZvRange();
  82   SetMeanXYZv();
  83   SetSigmaMeanXYZv();
  84   SetRedoTPCVertex();
  85   SetUseBeamSpotConstraint();
  86 }
  87
  88 //_____________________________________________________________________________
  89 Bool_t AliFilteredTreeEventCuts::AcceptEvent(AliESDEvent *esdEvent,AliMCEvent *mcEvent, const AliESDVertex *vtx)
  90 {
  91   // Check event selection cuts
  92   Bool_t retValue=kTRUE;
  93
  94   if(!esdEvent) return kFALSE;
  95   if(!IsRecVertexRequired()) return kTRUE;
  96   if(!vtx) return kFALSE;
  97   if(!vtx->GetStatus()) return kFALSE;
  98
  99   if(mcEvent) {
 100    // check MC event conditions
 101    AliHeader* header = mcEvent->Header();
 102    if(!header) return kFALSE;
 103
 104     // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive)
 105     if(fEventProcessType == kInvalidProcess) {
 106       retValue=kTRUE;
 107     }
 108     else if(fEventProcessType == kSD || fEventProcessType == kDD) {
 109       MCProcessType processType = GetEventProcessType(header);
 110       if(processType == kND) retValue=kFALSE;
 111       else retValue=kTRUE;
 112     }
 113     else if(fEventProcessType == GetEventProcessType(header)) {
 114       retValue=kTRUE;
 115     }
 116     else
 117       retValue=kFALSE;
 118   }
 119
 120   if(vtx->GetZv() < fMinZv) return kFALSE;
 121   if(vtx->GetZv() > fMaxZv) return kFALSE;
 122
 123 return retValue;
 124 }
 125
 126 //_____________________________________________________________________________
 127 Bool_t AliFilteredTreeEventCuts::AcceptMCEvent(AliMCEvent *mcEvent)
 128 {
 129   // Check event selection cuts
 130   if(!mcEvent) return kFALSE;
 131
 132   Bool_t retValue=kTRUE;
 133
 134   // check MC event conditions
 135   AliHeader* header = mcEvent->Header();
 136   if(!header) return kFALSE;
 137
 138   AliGenEventHeader* genHeader = header->GenEventHeader();
 139   if (!genHeader) {
 140     AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 141     return kFALSE;
 142   }
 143   TArrayF vtxMC(3);
 144   genHeader->PrimaryVertex(vtxMC);
 145
 146   // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive)
 147   if(fEventProcessType == kInvalidProcess) {
 148      retValue=kTRUE;
 149   } else {
 150      if(fEventProcessType == GetEventProcessType(header)) retValue=kTRUE;
 151      else retValue=kFALSE;
 152   }
 153
 154   /*
 155   Float_t R = TMath::Sqrt(vtxMC[0]*vtxMC[0]+vtxMC[1]*vtxMC[1]);
 156   if(R > fMaxR) return kFALSE;
 157   */
 158
 159   if(vtxMC[2] < fMinZv) return kFALSE;
 160   if(vtxMC[2] > fMaxZv) return kFALSE;
 161
 162 return retValue;
 163 }
 164
 165 //_____________________________________________________________________________
 166 Long64_t AliFilteredTreeEventCuts::Merge(TCollection* list)
 167 {
 168   // Merge list of objects (needed by PROOF)
 169   if (!list)
 170   return 0;
 171
 172   if (list->IsEmpty())
 173   return 1;
 174
 175   TIterator* iter = list->MakeIterator();
 176   TObject* obj = 0;
 177
 178   Int_t count=0;
 179   while((obj = iter->Next()) != 0)
 180   {
 181     AliFilteredTreeEventCuts* entry = dynamic_cast<AliFilteredTreeEventCuts*>(obj);
 182     if (entry == 0)
 183       continue;
 184
 185   count++;
 186   }
 187
 188 return count;
 189 }
 190
 191 //_____________________________________________________________________________
 192 AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
 193   //
 194   // get the process type of the event.
 195   //
 196
 197
 198   // Check for simple headers first
 199
 200   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
 201   if (pythiaGenHeader) {
 202     return GetPythiaEventProcessType(pythiaGenHeader,adebug);
 203   }
 204
 205   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader());
 206   if (dpmJetGenHeader) {
 207     return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
 208   }
 209
 210
 211   // check for cocktail
 212
 213   AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
 214   if (!genCocktailHeader) {
 215     printf("AliFilteredTreeEventCuts::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
 216     return kInvalidProcess;
 217   }
 218
 219   TList* headerList = genCocktailHeader->GetHeaders();
 220   if (!headerList) {
 221     return kInvalidProcess;
 222   }
 223
 224   for (Int_t i=0; i<headerList->GetEntries(); i++) {
 225
 226     pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 227     if (pythiaGenHeader) {
 228       return GetPythiaEventProcessType(pythiaGenHeader,adebug);
 229     }
 230
 231     dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(headerList->At(i));
 232     if (dpmJetGenHeader) {
 233       return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
 234     }
 235   }
 236   return kInvalidProcess;
 237 }
 238
 239 //____________________________________________________________________
 240 AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, DiffTreatment diffTreatment)
 241 {
 242   //
 243   // get process type
 244   //
 245   // diffTreatment:
 246   //   kMCFlags: use MC flags
 247   //   kUA5Cuts: SD events are those that have the MC SD flag and fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND
 248   //   kE710Cuts: SD events are those that have the MC SD flag and fulfill 2 < M^2 < 0.05s; DD from MC flag; Remainder is ND
 249   //   kALICEHadronLevel: SD events are those that fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND
 250   //
 251
 252   MCProcessType mcProcessType = GetEventProcessType(header);
 253
 254   if (diffTreatment == kMCFlags)
 255     return mcProcessType;
 256
 257   if (!esd)
 258   {
 259     Printf("ERROR: AliFilteredTreeEventCuts::GetEventProcessType: diffTreatment != kMCFlags and esd == 0");
 260     return kInvalidProcess;
 261   }
 262
 263   Float_t cms = esd->GetESDRun()->GetBeamEnergy();
 264   if (esd->GetESDRun()->IsBeamEnergyIsSqrtSHalfGeV())
 265     cms *= 2;
 266   //Printf("cms = %f", cms);
 267
 268   if (diffTreatment == kUA5Cuts && mcProcessType == kSD)
 269   {
 270     if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
 271       return kSD;
 272   }
 273   else if (diffTreatment == kE710Cuts && mcProcessType == kSD)
 274   {
 275     if (IsHadronLevelSingleDiffractive(stack, cms, 2. / cms / cms, 0.05))
 276       return kSD;
 277   }
 278   else if (diffTreatment == kALICEHadronLevel)
 279   {
 280     if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05))
 281       return kSD;
 282   }
 283
 284   if (mcProcessType == kSD)
 285     return kND;
 286
 287   return mcProcessType;
 288 }
 289
 290
 291 AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
 292
 293   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
 294
 295   if (!pythiaGenHeader) {
 296     printf("AliFilteredTreeEventCuts::GetProcessType : Unknown gen Header type). \n");
 297     return kInvalidProcess;
 298   }
 299
 300
 301   Int_t pythiaType = pythiaGenHeader->ProcessType();
 302   fgLastProcessType = pythiaType;
 303   MCProcessType globalType = kInvalidProcess;
 304
 305
 306   if (adebug) {
 307     printf("AliFilteredTreeEventCuts::GetProcessType : Pythia process type found: %d \n",pythiaType);
 308   }
 309
 310
 311   if(pythiaType==92||pythiaType==93){
 312     globalType = kSD;
 313   }
 314   else if(pythiaType==94){
 315     globalType = kDD;
 316   }
 317   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??}
 318   else {
 319     globalType = kND;
 320   }
 321   return globalType;
 322 }
 323
 324
 325 AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
 326   //
 327   // get the process type of the event.
 328   //
 329
 330   // can only read pythia headers, either directly or from cocktalil header
 331   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);
 332
 333   if (!dpmJetGenHeader) {
 334     printf("AliFilteredTreeEventCuts::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
 335     return kInvalidProcess;
 336   }
 337
 338   Int_t dpmJetType = dpmJetGenHeader->ProcessType();
 339   fgLastProcessType = dpmJetType;
 340   MCProcessType globalType = kInvalidProcess;
 341
 342
 343   if (adebug) {
 344     printf("AliFilteredTreeEventCuts::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
 345   }
 346
 347
 348   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
 349     globalType = kND;
 350   }
 351   else if (dpmJetType==5 || dpmJetType==6) {
 352     globalType = kSD;
 353   }
 354   else if (dpmJetType==7) {
 355     globalType = kDD;
 356   }
 357   return globalType;
 358 }
 359
 360 //____________________________________________________________________
 361 Bool_t AliFilteredTreeEventCuts::IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax)
 362 {
 363   //
 364   // return if process is single diffractive on hadron level
 365   //
 366   // xiMax and xiMin cut on M^2/s
 367   //
 368   // Based on code from Martin Poghoysan
 369   //
 370
 371   TParticle* part1 = 0;
 372   TParticle* part2 = 0;
 373
 374   Double_t smallestY = 1e10;
 375   Double_t largestY  = -1e10;
 376
 377   for (Int_t iParticle = 0; iParticle < stack->GetNprimary(); iParticle++)
 378   {
 379     TParticle* part = stack->Particle(iParticle);
 380     if (!part)
 381       continue;
 382
 383     Int_t pdg = TMath::Abs(part->GetPdgCode());
 384
 385     Int_t child1 = part->GetFirstDaughter();
 386     Int_t ist    = part->GetStatusCode();
 387
 388     Int_t mfl  = Int_t (pdg / TMath::Power(10, Int_t(TMath::Log10(pdg))));
 389     if (child1 > -1 || ist != 1)
 390       mfl = 0; // select final state charm and beauty
 391     if (!(stack->IsPhysicalPrimary(iParticle) || pdg == 111 || pdg == 3212 || pdg==3124 || mfl >= 4))
 392       continue;
 393     Int_t imother = part->GetFirstMother();
 394     if (imother>0)
 395     {
 396       TParticle *partM = stack->Particle(imother);
 397       Int_t pdgM=TMath::Abs(partM->GetPdgCode());
 398       if (pdgM==111 || pdgM==3124 || pdgM==3212)
 399         continue;
 400     }
 401
 402     Double_t y = 0;
 403
 404     // fix for problem with getting mass of particle 3124
 405     if (pdg != 3124)
 406       y = Rapidity(part->Pt(), part->Pz(), part->GetMass());
 407     else
 408       y = Rapidity(part->Pt(), part->Pz(), 1.5195);
 409
 410     if (y < smallestY)
 411     {
 412       smallestY = y;
 413       part1 = part;
 414     }
 415
 416     if (y > largestY)
 417     {
 418       largestY = y;
 419       part2 = part;
 420     }
 421   }
 422
 423   if (part1 == 0 || part2 == 0)
 424     return kFALSE;
 425
 426   Int_t pdg1 = part1->GetPdgCode();
 427   Int_t pdg2 = part2->GetPdgCode();
 428
 429   Double_t pt1 = part1->Pt();
 430   Double_t pt2 = part2->Pt();
 431   Double_t pz1 = part1->Pz();
 432   Double_t pz2 = part2->Pz();
 433
 434   Double_t y1 = TMath::Abs(Rapidity(pt1, pz1, 0.938));
 435   Double_t y2 = TMath::Abs(Rapidity(pt2, pz2, 0.938));
 436
 437   Int_t arm = -99999;
 438   if (pdg1 == 2212 && pdg2 == 2212)
 439   {
 440     if (y1 > y2)
 441       arm = 0;
 442     else
 443       arm = 1;
 444   }
 445   else if (pdg1 == 2212)
 446     arm = 0;
 447   else if (pdg2 == 2212)
 448     arm = 1;
 449
 450   Double_t M02s = 1. - 2 * part1->Energy() / cms;
 451   Double_t M12s = 1. - 2 * part2->Energy() / cms;
 452
 453   if (arm == 0 && M02s > xiMin && M02s < xiMax)
 454     return kTRUE;
 455   else if (arm == 1 && M12s > xiMin && M12s < xiMax)
 456     return kTRUE;
 457
 458   return kFALSE;
 459 }
 460
 461 //____________________________________________________________________
 462 Double_t AliFilteredTreeEventCuts::Rapidity(Double_t pt, Double_t pz, Double_t m)
 463 {
 464   //
 465   // calculates rapidity keeping the sign in case E == pz
 466   //
 467
 468   Double_t energy = TMath::Sqrt(pt*pt+pz*pz+m*m);
 469   if (energy != TMath::Abs(pz))
 470     return 0.5*TMath::Log((energy+pz)/(energy-pz));
 471
 472   Printf("W- mt=0");
 473   return TMath::Sign(1.e30,pz);
 474 }
 475