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(1)
  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   fRefObj.SetOwner(kTRUE);
  66 }
  67
  68 //______________________________________________
  69 AliDielectronHF::AliDielectronHF(const char* name, const char* title) :
  70   TNamed(name, title),
  71   fArrPairType(0x0),
  72   fPairType(kOSonly),
  73   fSignalsMC(0x0),
  74   fAxes(kMaxCuts),
  75   fHasMC(kFALSE),
  76   fStepGenerated(kFALSE),
  77   fRefObj(1)
  78 {
  79   //
  80   // Named Constructor
  81   //
  82   for (Int_t i=0; i<kMaxCuts; ++i){
  83     fVarCuts[i]=0;
  84     fVarCutType[i]=0;
  85     fBinType[i]=kStdBin;
  86   }
  87   fAxes.SetOwner(kTRUE);
  88   fRefObj.SetOwner(kTRUE);
  89 }
  90
  91 //______________________________________________
  92 AliDielectronHF::~AliDielectronHF()
  93 {
  94   //
  95   // Default Destructor
  96   //
  97   fAxes.Delete();
  98 }
  99
 100 //________________________________________________________________
 101 void AliDielectronHF::AddRefHist(TH1 *obj, UInt_t vars[4])
 102 {
 103   //
 104   // store reference object and its varaibles
 105   //
 106
 107   //  UInt_t val[2]={AliDielectronVarManager::kM,AliDielectronVarManager::kPt};
 108   AliDielectronHistos::StoreVariables(obj,vars);
 109   AliDielectronHistos::AdaptNameTitle(obj,"Pair");
 110   obj->SetName(Form("HF_%s",obj->GetName()));
 111   fRefObj.AddLast(obj);
 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   // fill only if we have asked for these steps
 186   if(!fStepGenerated) return;
 187
 188   AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
 189   AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
 190
 191   AliDielectronMC* dieMC = AliDielectronMC::Instance();
 192
 193   Int_t mLabel1 = dieMC->GetMothersLabel(label1);    // should work for both ESD and AOD
 194   Int_t mLabel2 = dieMC->GetMothersLabel(label2);
 195
 196   // check the same mother option
 197   AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
 198   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
 199   if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
 200
 201   // fill the leg variables
 202   Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
 203   Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
 204   AliDielectronVarManager::Fill(part1,valuesLeg1);
 205   AliDielectronVarManager::Fill(part2,valuesLeg2);
 206
 207   // fill the pair and event variables
 208   Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
 209   AliDielectronVarManager::Fill(dieMC->GetMCEvent(), valuesPair);
 210   AliDielectronVarManager::FillVarMCParticle2(part1,part2,valuesPair);
 211
 212   if(part1->Charge()*part2->Charge()<0) {
 213     //valuesPair[AliDielectronVarManager::kPairType]=1;
 214 //     for(Int_t i=0; i<fAxes.GetEntriesFast(); i++) {
 215 //       printf("[D]: %s %f %f %f \n",
 216 //           AliDielectronVarManager::GetValueName(fVarCuts[i]),
 217 //           valuesLeg1[fVarCuts[i]], valuesLeg2[fVarCuts[i]], valuesPair[fVarCuts[i]]);
 218 //     }
 219     Fill(nSignal+fSignalsMC->GetEntries(), valuesPair,  valuesLeg1, valuesLeg2);
 220   }
 221   // only OS at the moment
 222   // else if(part1->Charge()>0)
 223   //   valuesPair[AliDielectronVarManager::kPairType]=0;
 224   // else
 225   //   valuesPair[AliDielectronVarManager::kPairType]=2; // if one of the two particles is neutral, the pair will go here
 226
 227   return;
 228 }
 229 //______________________________________________
 230 void AliDielectronHF::Fill(Int_t pairIndex, const AliDielectronPair *particle)
 231 {
 232   //
 233   // fill histograms for event, pair and daughter cuts and pair types
 234   //
 235
 236   // only OS pairs in case of MC
 237   if(fHasMC && pairIndex!=AliDielectron::kEv1PM) return;
 238
 239   // only selected pair types in case of data
 240   if(!IsPairTypeSelected(pairIndex)) return;
 241
 242   // get event and pair variables
 243   Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
 244   AliDielectronVarManager::Fill(particle,valuesPair);
 245
 246   // get leg variables
 247   Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]={0};
 248   AliDielectronVarManager::Fill(particle->GetFirstDaughter(),valuesLeg1);
 249   Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]={0};
 250   AliDielectronVarManager::Fill(particle->GetSecondDaughter(),valuesLeg2);
 251
 252   // fill
 253   if(!fHasMC) { Fill(pairIndex, valuesPair,  valuesLeg1, valuesLeg2); }
 254   if(fHasMC && fSignalsMC) {
 255     for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 256       if(AliDielectronMC::Instance()->IsMCTruth(particle, (AliDielectronSignalMC*)fSignalsMC->At(i)))
 257         Fill(i, valuesPair,  valuesLeg1, valuesLeg2);
 258     }
 259   }
 260
 261   return;
 262
 263 }
 264
 265 //______________________________________________
 266 void AliDielectronHF::Fill(Int_t Index, Double_t * const valuesPair, Double_t * const valuesLeg1, Double_t * const valuesLeg2)
 267 {
 268   //
 269   // main fill function using index and values as input
 270   //
 271
 272   TObjArray *histArr = static_cast<TObjArray*>(fArrPairType.At(Index));
 273   if(!histArr) return;
 274
 275   Int_t size  = GetNumberOfBins();
 276   // loop over all histograms
 277   for(Int_t ihist=0; ihist<size; ihist++) {
 278
 279     Int_t sizeAdd   = 1;
 280     Bool_t selected = kTRUE;
 281
 282     // loop over all cut variables
 283     Int_t nvars = fAxes.GetEntriesFast();
 284     for(Int_t ivar=0; ivar<nvars; ivar++) {
 285
 286       // get bin limits
 287       TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
 288       Int_t nbins    = bins->GetNrows()-1;
 289
 290       // bin limits for current ivar bin
 291       Int_t ibin   = (ihist/sizeAdd)%nbins;
 292       Double_t lowEdge = (*bins)[ibin];
 293       Double_t upEdge  = (*bins)[ibin+1];
 294       switch(fBinType[ivar]) {
 295       case kStdBin:     upEdge=(*bins)[ibin+1];     break;
 296       case kBinToMax:   upEdge=(*bins)[nbins];      break;
 297       case kBinFromMin: lowEdge=(*bins)[0];         break;
 298       case kSymBin:     upEdge=(*bins)[nbins-ibin];
 299         if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
 300         break;
 301       }
 302
 303       // leg variable
 304       if(fVarCutType[ivar]) {
 305         if( (valuesLeg1[fVarCuts[ivar]] < lowEdge || valuesLeg1[fVarCuts[ivar]] >= upEdge) ||
 306             (valuesLeg2[fVarCuts[ivar]] < lowEdge || valuesLeg2[fVarCuts[ivar]] >= upEdge) ) {
 307           selected=kFALSE;
 308           break;
 309         }
 310       }
 311       else { // pair and event variables
 312         if( (valuesPair[fVarCuts[ivar]] < lowEdge || valuesPair[fVarCuts[ivar]] >= upEdge) ) {
 313           selected=kFALSE;
 314           break;
 315         }
 316       }
 317
 318       sizeAdd*=nbins;
 319     } //end of var cut loop
 320
 321     // do not fill the histogram
 322     if(!selected) continue;
 323
 324     // fill the object with Pair and event values (TODO: this needs to be changed)
 325     TObjArray *tmp = (TObjArray*) histArr->At(ihist);
 326     TString title = tmp->GetName();
 327     AliDebug(10,title.Data());
 328     for(Int_t i=0; i<tmp->GetEntriesFast(); i++) {
 329       AliDielectronHistos::FillValues(tmp->At(i), valuesPair);
 330     }
 331     //    AliDebug(10,Form("Fill var %d %s value %f in %s \n",fVar,AliDielectronVarManager::GetValueName(fVar),valuesPair[fVar],tmp->GetName()));
 332   } //end of hist loop
 333
 334 }
 335
 336 //______________________________________________
 337 void AliDielectronHF::Init()
 338 {
 339   //
 340   // initialise event buffers
 341   //
 342
 343   // has MC signals
 344   fHasMC=AliDielectronMC::Instance()->HasMC();
 345   Int_t steps = 0;
 346   if(fHasMC) steps=fSignalsMC->GetEntries();
 347   if(fStepGenerated) steps*=2;
 348
 349   // init pair type array
 350   fArrPairType.SetName(Form("%s_HF",GetName()));
 351   if(fHasMC) fArrPairType.Expand(steps);
 352   else fArrPairType.Expand(AliDielectron::kEv1PMRot+1);
 353
 354   Int_t size  = GetNumberOfBins();
 355   AliDebug(10,Form("Creating a histo array with size %d \n",size));
 356
 357   Int_t sizeAdd  = 1;
 358
 359   // fill object array with the histograms
 360   TObjArray *histArr = new TObjArray();
 361   histArr->Expand(size);
 362
 363   //  printf("fRefObj %p \n",fRefObj);
 364   for(Int_t ihist=0; ihist<size; ihist++) {
 365     histArr->AddAt(fRefObj.Clone(""), ihist);
 366     //histArr->AddAt(fRefObj.Clone(Form("h%04d",ihist)), ihist);
 367   }
 368
 369   // loop over all cut variables
 370   Int_t nvars = fAxes.GetEntriesFast();
 371   for(Int_t ivar=0; ivar<nvars; ivar++) {
 372
 373     // get bin limits
 374     TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
 375     Int_t nbins    = bins->GetNrows()-1;
 376
 377
 378     // loop over all histograms an set unique titles
 379     for(Int_t ihist=0; ihist<size; ihist++) {
 380
 381       // get the lower limit for current ivar bin
 382       Int_t ibin   = (ihist/sizeAdd)%nbins;
 383       Double_t lowEdge = (*bins)[ibin];
 384       Double_t upEdge  = (*bins)[ibin+1];
 385       switch(fBinType[ivar]) {
 386       case kStdBin:     upEdge=(*bins)[ibin+1];     break;
 387       case kBinToMax:   upEdge=(*bins)[nbins];      break;
 388       case kBinFromMin: lowEdge=(*bins)[0];         break;
 389       case kSymBin:     upEdge=(*bins)[nbins-ibin];
 390         if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins]; // to avoid low>up
 391         break;
 392       }
 393
 394       TObjArray *tmp= (TObjArray*) histArr->At(ihist);
 395       TString title = tmp->GetName();
 396       if(!ivar)             title ="";
 397       if( ivar)             title+=":";
 398       if(fVarCutType[ivar]) title+="Leg";
 399       title+=AliDielectronVarManager::GetValueName(fVarCuts[ivar]);
 400       title+=Form("#%.2f#%.2f",lowEdge,upEdge);
 401       tmp->SetName(title.Data());
 402       AliDebug(10,title.Data());
 403     }
 404     sizeAdd*=nbins;
 405   }
 406
 407   // copy array to the selected pair types/ MC sources
 408   if(fHasMC) {
 409     for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
 410         TString title = Form("(Signal: %s)",fSignalsMC->At(i)->GetTitle());
 411         fArrPairType[i]=(TObjArray*)histArr->Clone(title.Data());
 412         if(fStepGenerated)  {
 413           title+=" MC truth";
 414           fArrPairType[i+fSignalsMC->GetEntries()]=(TObjArray*)histArr->Clone(title.Data());
 415         }
 416       }
 417    }
 418    else {
 419     for(Int_t i=0; i<AliDielectron::kEv1PMRot+1; i++) {
 420       if(IsPairTypeSelected(i)) fArrPairType[i]=(TObjArray*)histArr->Clone(Form("%s",AliDielectron::PairClassName(i)));
 421       else fArrPairType[i]=0x0;
 422     }
 423   }
 424
 425   // clean up
 426   if(histArr) {
 427     delete histArr;
 428     histArr=0;
 429   }
 430
 431 }
 432
 433 //______________________________________________
 434 Int_t AliDielectronHF::GetNumberOfBins() const
 435 {
 436   //
 437   // return the number of bins this mixing handler has
 438   //
 439   Int_t size=1;
 440   for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i)
 441     size*=((static_cast<TVectorD*>(fAxes.At(i)))->GetNrows()-1);
 442   return size;
 443 }
 444
 445 //______________________________________________
 446 Bool_t AliDielectronHF::IsPairTypeSelected(Int_t itype)
 447 {
 448   //
 449   // check whether a pair type was selected
 450   //
 451
 452   Bool_t selected = kFALSE;
 453
 454   // fill all
 455   if(fPairType==kAll) selected = kTRUE;
 456
 457   switch(itype) {
 458   case AliDielectron::kEv1PP:
 459   case AliDielectron::kEv1MM:
 460     if(fPairType==kOSandLS)   selected = kTRUE;
 461     break;
 462   case AliDielectron::kEv1PM: selected = kTRUE;
 463     break;
 464   case AliDielectron::kEv1PEv2P:
 465   case AliDielectron::kEv1MEv2P:
 466   case AliDielectron::kEv1PEv2M:
 467   case AliDielectron::kEv1MEv2M:
 468     if(fPairType==kOSandMIX)  selected = kTRUE;
 469     break;
 470   case AliDielectron::kEv2PP:
 471   case AliDielectron::kEv2PM:
 472   case AliDielectron::kEv2MM:
 473     selected = kFALSE;
 474     break;
 475   case AliDielectron::kEv1PMRot:
 476     if(fPairType==kOSandROT)  selected = kTRUE;
 477     break;
 478   }
 479
 480   return selected;
 481
 482 }
 483
 484