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