PWGDQ/dielectron/AliDielectronHF.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, 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 ///////////////////////////////////////////////////////////////////////////
  17 //                Dielectron HF                                  //
  18 //                                                                       //
  19 //                                                                       //
  20 /*
  21 Detailed description
  22
  23
  24 */
  25 //                                                                       //
  26 ///////////////////////////////////////////////////////////////////////////
  27
  28 #include <TVectorD.h>
  29 #include <TH1.h>
  30 #include <TAxis.h>
  31 #include <AliVParticle.h>
  32
  33 #include <AliLog.h>
  34
  35 #include "AliDielectron.h"
  36 #include "AliDielectronHelper.h"
  37 #include "AliDielectronMC.h"
  38 #include "AliDielectronPair.h"
  39 #include "AliDielectronSignalMC.h"
  40
  41 #include "AliDielectronHistos.h"
  42 #include "AliDielectronHF.h"
  43
  44 ClassImp(AliDielectronHF)
  45
  46 AliDielectronHF::AliDielectronHF() :
  47   TNamed(),
  48   fArrPairType(0x0),
  49   fPairType(kOSonly),
  50   fSignalsMC(0x0),
  51   fAxes(kMaxCuts),
  52   fHasMC(kFALSE),
  53   fStepGenerated(kFALSE),
  54   fRefObj(0x0)
  55 {
  56   //
  57   // Default Constructor
  58   //
  59   for (Int_t i=0; i<kMaxCuts; ++i){
  60     fVarCuts[i]=0;
  61     fVarCutType[i]=0;
  62     fBinType[i]=kStdBin;
  63   }
  64   fAxes.SetOwner(kTRUE);
  65 }
  66
  67 //______________________________________________
  68 AliDielectronHF::AliDielectronHF(const char* name, const char* title) :
  69   TNamed(name, title),
  70   fArrPairType(0x0),
  71   fPairType(kOSonly),
  72   fSignalsMC(0x0),
  73   fAxes(kMaxCuts),
  74   fHasMC(kFALSE),
  75   fStepGenerated(kFALSE),
  76   fRefObj(0x0)
  77 {
  78   //
  79   // Named Constructor
  80   //
  81   for (Int_t i=0; i<kMaxCuts; ++i){
  82     fVarCuts[i]=0;
  83     fVarCutType[i]=0;
  84     fBinType[i]=kStdBin;
  85   }
  86   fAxes.SetOwner(kTRUE);
  87 }
  88
  89 //______________________________________________
  90 AliDielectronHF::~AliDielectronHF()
  91 {
  92   //
  93   // Default Destructor
  94   //
  95   fAxes.Delete();
  96 }
  97
  98 //________________________________________________________________
  99 void AliDielectronHF::SetRefHist(TH1 *obj, UInt_t vars[4])
 100 {
 101   //
 102   // store reference object and its varaibles
 103   //
 104
 105   //  UInt_t val[2]={AliDielectronVarManager::kM,AliDielectronVarManager::kPt};
 106   AliDielectronHistos::StoreVariables(obj,vars);
 107   AliDielectronHistos::AdaptNameTitle(obj,"Pair");
 108   obj->SetName("");
 109   fRefObj     = obj;
 110
 111 }
 112
 113 //________________________________________________________________
 114 void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
 115                                              Int_t nbins, Double_t min, Double_t max, Bool_t log, Bool_t leg, EBinType btype)
 116 {
 117   //
 118   // Add a variable to the mixing handler
 119   //
 120
 121   // limit number of variables to kMaxCuts
 122   if (fAxes.GetEntriesFast()>=kMaxCuts) return;
 123
 124   TVectorD *binLimits=0x0;
 125   if (!log) binLimits=AliDielectronHelper::MakeLinBinning(nbins,min,max);
 126   else binLimits=AliDielectronHelper::MakeLogBinning(nbins,min,max);
 127   if (!binLimits) return;
 128
 129   Int_t size=fAxes.GetEntriesFast();
 130   fVarCuts[size]=(UShort_t)type;
 131   fVarCutType[size]=leg;
 132   fAxes.Add(binLimits->Clone());
 133   fBinType[size]=btype;
 134 }
 135
 136 //________________________________________________________________
 137 void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
 138                                              const char* binLimitStr, Bool_t leg, EBinType btype)
 139 {
 140   //
 141   // Add a variable to the mixing handler with arbitrary binning
 142   //
 143
 144   // limit number of variables to kMaxCuts
 145   if (fAxes.GetEntriesFast()>=kMaxCuts) return;
 146
 147   TVectorD *binLimits=AliDielectronHelper::MakeArbitraryBinning(binLimitStr);
 148   if (!binLimits) return;
 149
 150   Int_t size=fAxes.GetEntriesFast();
 151   fVarCuts[size]=(UShort_t)type;
 152   fVarCutType[size]=leg;
 153   fAxes.Add(binLimits);
 154   fBinType[size]=btype;
 155 }
 156
 157 //________________________________________________________________
 158 void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
 159                                              TVectorD * binLimits, Bool_t leg, EBinType btype)
 160 {
 161   //
 162   // Add a variable to the mixing handler with a vector
 163   // the TVectorD is assumed to be surplus after the creation and will be deleted!!!
 164   //
 165
 166   // limit number of variables to kMaxCuts
 167   if (fAxes.GetEntriesFast()>=kMaxCuts) return;
 168
 169   if (!binLimits) return;
 170
 171   Int_t size=fAxes.GetEntriesFast();
 172   fVarCuts[size]=(UShort_t)type;
 173   fVarCutType[size]=leg;
 174   fAxes.Add(binLimits);
 175   fBinType[size]=btype;
 176 }
 177
 178 //______________________________________________
 179 void AliDielectronHF::Fill(Int_t label1, Int_t label2, Int_t nSignal)
 180 {
 181   //
 182   // fill the pure MC part of the container starting from a pair of 2 particles (part1 and part2 are legs)
 183   //
 184   // fill only if we have asked for these steps
 185   if(!fStepGenerated) return;
 186
 187   AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
 188   AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
 189
 190   AliDielectronMC* dieMC = AliDielectronMC::Instance();
 191
 192   Int_t mLabel1 = dieMC->GetMothersLabel(label1);    // should work for both ESD and AOD
 193   Int_t mLabel2 = dieMC->GetMothersLabel(label2);
 194
 195   // check the same mother option
 196   AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
 197   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
 198   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
 199
 200   // fill the leg variables
 201   Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
 202   Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
 203   AliDielectronVarManager::Fill(part1,valuesLeg1);
 204   AliDielectronVarManager::Fill(part2,valuesLeg2);
 205
 206   // fill the pair and event variables
 207   Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
 208   AliDielectronVarManager::Fill(dieMC->GetMCEvent(), valuesPair);
 209   AliDielectronVarManager::FillVarMCParticle2(part1,part2,valuesPair);
 210
 211   if(part1->Charge()*part2->Charge()<0) {
 212     //valuesPair[AliDielectronVarManager::kPairType]=1;
 213 //     for(Int_t i=0; i<fAxes.GetEntriesFast(); i++) {
 214 //       printf("[D]: %s %f %f %f \n",
 215 //           AliDielectronVarManager::GetValueName(fVarCuts[i]),
 216 //           valuesLeg1[fVarCuts[i]], valuesLeg2[fVarCuts[i]], valuesPair[fVarCuts[i]]);
 217 //     }
 218     Fill(nSignal, valuesPair,  valuesLeg1, valuesLeg2);
 219   }
 220   // on OS at the moment
 221   // else if(part1->Charge()>0)
 222   //   valuesPair[AliDielectronVarManager::kPairType]=0;
 223   // else
 224   //   valuesPair[AliDielectronVarManager::kPairType]=2; // if one of the two particles is neutral, the pair will go here
 225
 226   return;
 227 }
 228 //______________________________________________
 229 void AliDielectronHF::Fill(Int_t pairIndex, const AliDielectronPair *particle)
 230 {
 231   //
 232   // fill histograms for event, pair and daughter cuts and pair types
 233   //
 234
 235   // only OS pairs in case of MC
 236   if(fHasMC && pairIndex!=AliDielectron::kEv1PM) return;
 237
 238   // only selected pair types in case of data
 239   if(!IsPairTypeSelected(pairIndex)) return;
 240
 241   // get event and pair variables
 242   Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
 243   AliDielectronVarManager::Fill(particle,valuesPair);
 244
 245   // get leg variables
 246   Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]={0};
 247   AliDielectronVarManager::Fill(particle->GetFirstDaughter(),valuesLeg1);
 248   Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]={0};
 249   AliDielectronVarManager::Fill(particle->GetSecondDaughter(),valuesLeg2);
 250
 251   // fill
 252   if(!fHasMC) { Fill(pairIndex, valuesPair,  valuesLeg1, valuesLeg2); }
 253   if(fHasMC && fSignalsMC) {
 254     for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 255       if(AliDielectronMC::Instance()->IsMCTruth(particle, (AliDielectronSignalMC*)fSignalsMC->At(i)))
 256         Fill(i, valuesPair,  valuesLeg1, valuesLeg2);
 257     }
 258   }
 259
 260   return;
 261
 262 }
 263
 264 //______________________________________________
 265 void AliDielectronHF::Fill(Int_t Index, Double_t * const valuesPair, Double_t * const valuesLeg1, Double_t * const valuesLeg2)
 266 {
 267   //
 268   // main fill function using index and values as input
 269   //
 270
 271   TObjArray *histArr = static_cast<TObjArray*>(fArrPairType.At(Index));
 272   if(!histArr) return;
 273
 274   Int_t size  = GetNumberOfBins();
 275   // loop over all histograms
 276   for(Int_t ihist=0; ihist<size; ihist++) {
 277
 278     Int_t sizeAdd   = 1;
 279     Bool_t selected = kTRUE;
 280
 281     // loop over all cut variables
 282     Int_t nvars = fAxes.GetEntriesFast();
 283     for(Int_t ivar=0; ivar<nvars; ivar++) {
 284
 285       // get bin limits
 286       TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
 287       Int_t nbins    = bins->GetNrows()-1;
 288
 289       // bin limits for current ivar bin
 290       Int_t ibin   = (ihist/sizeAdd)%nbins;
 291       Double_t lowEdge = (*bins)[ibin];
 292       Double_t upEdge  = (*bins)[ibin+1];
 293       switch(fBinType[ivar]) {
 294       case kStdBin:     upEdge=(*bins)[ibin+1];     break;
 295       case kBinToMax:   upEdge=(*bins)[nbins];      break;
 296       case kBinFromMin: lowEdge=(*bins)[0];         break;
 297       case kSymBin:     upEdge=(*bins)[nbins-ibin];
 298         if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
 299         break;
 300       }
 301
 302       // leg variable
 303       if(fVarCutType[ivar]) {
 304         if( (valuesLeg1[fVarCuts[ivar]] < lowEdge || valuesLeg1[fVarCuts[ivar]] >= upEdge) ||
 305             (valuesLeg2[fVarCuts[ivar]] < lowEdge || valuesLeg2[fVarCuts[ivar]] >= upEdge) ) {
 306           selected=kFALSE;
 307           break;
 308         }
 309       }
 310       else { // pair and event variables
 311         if( (valuesPair[fVarCuts[ivar]] < lowEdge || valuesPair[fVarCuts[ivar]] >= upEdge) ) {
 312           selected=kFALSE;
 313           break;
 314         }
 315       }
 316
 317       sizeAdd*=nbins;
 318     } //end of var cut loop
 319
 320     // do not fill the histogram
 321     if(!selected) continue;
 322
 323     // fill the object with Pair and event values (TODO: this needs to be changed)
 324     TH1 *tmp=static_cast<TH1*>(histArr->At(ihist));
 325     TString title = tmp->GetName();
 326     AliDebug(10,title.Data());
 327
 328     AliDielectronHistos::FillValues(tmp, valuesPair);
 329     //    AliDebug(10,Form("Fill var %d %s value %f in %s \n",fVar,AliDielectronVarManager::GetValueName(fVar),valuesPair[fVar],tmp->GetName()));
 330   } //end of hist loop
 331
 332 }
 333
 334 //______________________________________________
 335 void AliDielectronHF::Init()
 336 {
 337   //
 338   // initialise event buffers
 339   //
 340
 341   // has MC signals
 342   fHasMC=AliDielectronMC::Instance()->HasMC();
 343   Int_t steps = 0;
 344   if(fHasMC) steps=fSignalsMC->GetEntries();
 345   if(fStepGenerated) steps*=2;
 346
 347   // init pair type array
 348   fArrPairType.SetName(Form("%s_HF",GetName()));
 349   if(fHasMC) fArrPairType.Expand(steps);
 350   else fArrPairType.Expand(AliDielectron::kEv1PMRot+1);
 351
 352   Int_t size  = GetNumberOfBins();
 353   AliDebug(10,Form("Creating a histo array with size %d \n",size));
 354
 355   Int_t sizeAdd  = 1;
 356
 357   // fill object array with the histograms
 358   TObjArray *histArr = new TObjArray();
 359   histArr->Expand(size);
 360
 361   //  printf("fRefObj %p \n",fRefObj);
 362   for(Int_t ihist=0; ihist<size; ihist++) {
 363     histArr->AddAt(fRefObj->Clone(""), ihist);
 364     //histArr->AddAt(fRefObj->Clone(Form("h%04d",ihist)), ihist);
 365   }
 366
 367   // loop over all cut variables
 368   Int_t nvars = fAxes.GetEntriesFast();
 369   for(Int_t ivar=0; ivar<nvars; ivar++) {
 370
 371     // get bin limits
 372     TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
 373     Int_t nbins    = bins->GetNrows()-1;
 374
 375
 376     // loop over all histograms an set unique titles
 377     for(Int_t ihist=0; ihist<size; ihist++) {
 378
 379       // get the lower limit for current ivar bin
 380       Int_t ibin   = (ihist/sizeAdd)%nbins;
 381       Double_t lowEdge = (*bins)[ibin];
 382       Double_t upEdge  = (*bins)[ibin+1];
 383       switch(fBinType[ivar]) {
 384       case kStdBin:     upEdge=(*bins)[ibin+1];     break;
 385       case kBinToMax:   upEdge=(*bins)[nbins];      break;
 386       case kBinFromMin: lowEdge=(*bins)[0];         break;
 387       case kSymBin:     upEdge=(*bins)[nbins-ibin];
 388         if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
 389         break;
 390       }
 391
 392       TH1 *tmp=static_cast<TH1*>(histArr->At(ihist));
 393       TString title = tmp->GetName();
 394       if(ivar!=0)           title+=":";
 395       if(fVarCutType[ivar]) title+="Leg";
 396       title+=AliDielectronVarManager::GetValueName(fVarCuts[ivar]);
 397       title+=Form("#%.2f#%.2f",lowEdge,upEdge);
 398       tmp->SetName(title.Data());
 399       AliDebug(10,title.Data());
 400     }
 401     sizeAdd*=nbins;
 402   }
 403
 404   // copy array to the selected pair types/ MC sources
 405   if(fHasMC) {
 406     for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 407         TString title = Form("(Signal: %s)",fSignalsMC->At(i)->GetTitle());
 408         fArrPairType[i]=(TObjArray*)histArr->Clone(title.Data());
 409         if(fStepGenerated)  {
 410           title+=" MC truth";
 411           fArrPairType[i+fSignalsMC->GetEntries()]=(TObjArray*)histArr->Clone(title.Data());
 412         }
 413       }
 414    }
 415    else {
 416     for(Int_t i=0; i<AliDielectron::kEv1PMRot+1; i++) {
 417       if(IsPairTypeSelected(i)) fArrPairType[i]=(TObjArray*)histArr->Clone(Form("%s",AliDielectron::PairClassName(i)));
 418       else fArrPairType[i]=0x0;
 419     }
 420   }
 421
 422   // clean up
 423   if(histArr) {
 424     delete histArr;
 425     histArr=0;
 426   }
 427
 428 }
 429
 430 //______________________________________________
 431 Int_t AliDielectronHF::GetNumberOfBins() const
 432 {
 433   //
 434   // return the number of bins this mixing handler has
 435   //
 436   Int_t size=1;
 437   for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i)
 438     size*=((static_cast<TVectorD*>(fAxes.At(i)))->GetNrows()-1);
 439   return size;
 440 }
 441
 442 //______________________________________________
 443 Bool_t AliDielectronHF::IsPairTypeSelected(Int_t itype)
 444 {
 445   //
 446   // check whether a pair type was selected
 447   //
 448
 449   Bool_t selected = kFALSE;
 450
 451   // fill all
 452   if(fPairType==kAll) selected = kTRUE;
 453
 454   switch(itype) {
 455   case AliDielectron::kEv1PP:
 456   case AliDielectron::kEv1MM:
 457     if(fPairType==kOSandLS)   selected = kTRUE;
 458     break;
 459   case AliDielectron::kEv1PM: selected = kTRUE;
 460     break;
 461   case AliDielectron::kEv1PEv2P:
 462   case AliDielectron::kEv1MEv2P:
 463   case AliDielectron::kEv1PEv2M:
 464   case AliDielectron::kEv1MEv2M:
 465     if(fPairType==kOSandMIX)  selected = kTRUE;
 466     break;
 467   case AliDielectron::kEv2PP:
 468   case AliDielectron::kEv2PM:
 469   case AliDielectron::kEv2MM:
 470     selected = kFALSE;
 471     break;
 472   case AliDielectron::kEv1PMRot:
 473     if(fPairType==kOSandROT)  selected = kTRUE;
 474     break;
 475   }
 476
 477   return selected;
 478
 479 }
 480
 481