PWGLF/FORWARD/analysis2/AliForwardCreateResponseMatrices.cxx

   1 /**
   2  * @file   AliForwardCreateResponseMatrices.cxx
   3  * @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
   4  * @date   Thu Feb  7 01:01:24 2013
   5  *
   6  * @brief
   7  *
   8  *
   9  * @ingroup pwglf_forward_multdist
  10  */
  11
  12 #include <TH1D.h>
  13 #include "AliForwardCreateResponseMatrices.h"
  14 #include "AliForwardMultiplicityDistribution.h"
  15 #include "AliAODForwardMult.h"
  16 #include "AliAODCentralMult.h"
  17 #include "AliAODEvent.h"
  18 #include "AliFMDMCEventInspector.h"
  19 #include "AliAODMCHeader.h"
  20 #include "AliAODMCParticle.h"
  21
  22
  23 ClassImp(AliForwardCreateResponseMatrices)
  24 #if 0
  25 ; // This is for Emacs - do not delete
  26 #endif
  27 //_____________________________________________________________________
  28 AliForwardCreateResponseMatrices::AliForwardCreateResponseMatrices()
  29   : AliBaseAODTask(),
  30     fBins(),
  31     fIsSelected(false)
  32 {
  33   //
  34   // Default Constructor
  35   //
  36 }
  37
  38 //_____________________________________________________________________
  39 AliForwardCreateResponseMatrices::AliForwardCreateResponseMatrices(const char* name)
  40   : AliBaseAODTask(name,"AliForwardCreateResponseMatrices"),
  41     fBins(),
  42     fIsSelected(false)
  43 {
  44   //
  45   // Constructor
  46   //
  47 }
  48
  49 //_____________________________________________________________________
  50 Bool_t
  51 AliForwardCreateResponseMatrices::CheckEvent(const AliAODForwardMult& fwd)
  52 {
  53   fIsSelected = AliBaseAODTask::CheckEvent(fwd);
  54   // We always return true, so that we can process MC truth;
  55   return true;
  56 }
  57
  58 //_____________________________________________________________________
  59 Bool_t AliForwardCreateResponseMatrices::Book()
  60 {
  61   //
  62   // Create Output Objects
  63   //
  64   TH1D* dndetaSumForward = new TH1D("dndetaSumForward","dndetaSumForward",
  65                                     200,-4,6);
  66   TH1D* dndetaSumCentral = new TH1D("dndetaSumCentral","dndetaSumCentral",
  67                                     200,-4,6);
  68   TH1D* dndetaSumMC = new TH1D("dndetaSumMC","dndetaSumMC", 200,-4,6);
  69
  70   fSums->Add(dndetaSumForward);
  71   fSums->Add(dndetaSumCentral);
  72   fSums->Add(dndetaSumMC);
  73
  74   TH1D* dndetaEventForward = new TH1D();
  75   TH1D* dndetaEventCentral = new TH1D();
  76   TH1D* dndetaEventMC = new TH1D();
  77
  78   fSums->Add(dndetaEventForward);
  79   fSums->Add(dndetaEventCentral);
  80   fSums->Add(dndetaEventMC);
  81
  82   //Loop over all individual eta bins, and define their hisograms
  83   TIter next(&fBins);
  84   Bin * bin = 0;
  85   while ((bin = static_cast<Bin*>(next()))) {
  86     bin->CreateOutputObjectss(fSums, 400);
  87   }
  88
  89   return true;
  90
  91 }
  92
  93
  94 //_____________________________________________________________________
  95 Bool_t AliForwardCreateResponseMatrices::Event(AliAODEvent& aod)
  96 {
  97   //
  98   // User Exec
  99   //
 100   // Here, we used to get ForwardMC!
 101   AliAODForwardMult* AODforward = GetForward(aod);
 102   if (!AODforward) return false;
 103
 104   // Here, we used to get CentralClusters!
 105   AliAODCentralMult* AODcentral = GetCentral(aod);
 106   if (!AODcentral) return false;
 107
 108   TH2D& forward = AODforward->GetHistogram();
 109   TH2D& central = AODcentral->GetHistogram();
 110
 111   //check if event is NSD according to either MC truth or analysis (ESD)
 112   Bool_t isMCNSD  = kFALSE;
 113   Bool_t isESDNSD = kFALSE;
 114   if(AODforward->IsTriggerBits(AliAODForwardMult::kMCNSD))   isMCNSD=kTRUE;
 115   if(AODforward->IsTriggerBits(AliAODForwardMult::kNSD))     isESDNSD=kTRUE;
 116
 117   //primary dndeta dist
 118   TH2D*    primHist   = GetPrimary(aod);
 119
 120   TH1D* dndetaSumForward    = (TH1D*)fSums->FindObject("dndetaSumForward");
 121   TH1D* dndetaSumCentral    = (TH1D*)fSums->FindObject("dndetaSumCentral");
 122   TH1D* dndetaSumMC         = (TH1D*)fSums->FindObject("dndetaSumMC");
 123   TH1D* dndetaEventForward  = forward.ProjectionX("dndetaForward",1,
 124                                                   forward.GetNbinsY(),"");
 125   TH1D* dndetaEventCentral  = central.ProjectionX("dndetaCentral",1,
 126                                                   central.GetNbinsY(),"");
 127   TH1D* dndetaEventMC       = primHist->ProjectionX("dndetaMC",1,
 128                                                     primHist->GetNbinsY(),"");
 129
 130   // underflow eta bin of forward/central carry information on whether
 131   // there is acceptance. 1= acceptance, 0= no acceptance
 132   TH1D* normEventForward = 0;
 133   TH1D* normEventCentral = 0;
 134   normEventForward       = forward.ProjectionX("dndetaEventForwardNorm",0,0,"");
 135   normEventCentral       = central.ProjectionX("dndetaEventCentralNorm",0,0,"");
 136
 137   dndetaSumForward->Add(dndetaEventForward);
 138   dndetaSumCentral->Add(dndetaEventCentral);
 139   dndetaSumMC->Add(dndetaEventMC);
 140
 141
 142   Double_t VtxZ = AODforward->GetIpZ();
 143
 144   // loop over all eta bins, and create response matrices,
 145   // trigger-vertex bias histograms etc
 146   TIter next(&fBins);
 147   Bin * bin = 0;
 148   while ((bin = static_cast<Bin*>(next()))) {
 149     bin->Process(dndetaEventForward, dndetaEventCentral,
 150                  normEventForward,   normEventCentral,
 151                  dndetaEventMC, VtxZ, fIsSelected,
 152                  isMCNSD, isESDNSD,  aod);
 153   }
 154
 155   return true;
 156 }
 157
 158
 159 //=====================================================================
 160 AliForwardCreateResponseMatrices::Bin::Bin(Double_t etaLow, Double_t etaHigh)
 161   : TNamed("", ""),
 162     fEtaLow(etaLow),          // low eta limit
 163     fEtaHigh(etaHigh),        // high eta limit
 164     fHist(),                  // multiplicity histogram
 165     fHistMC(),                // multiplicity histogram MC truth primaries
 166     fAcceptance(),            // histogram showing the 'holes' in acceptance. BinContent of 1 shows a hole, and BinContent of 10 shows data coverage
 167     fVtxZvsNdataBins(),       // VtxZ vs. number of data acceptance bins (normalised to the eta range)
 168     fResponseMatrix(),        // Response matrix (MC truth vs. analysed multiplicity)
 169     fResponseMatrixPlus05(),  // Response matrix with analysed multiplicity scaled up by 5%
 170     fResponseMatrixPlus075(), // Response matrix  with analysed multiplicity scaled up by 7.5%
 171     fResponseMatrixPlus10(),  // Response matrix with analysed multiplicity scaled up by 10%
 172     fResponseMatrixMinus05(), // Response matrix with analysed multiplicity scaled down by 5%
 173     fResponseMatrixMinus075(),// Response matrix with analysed multiplicity scaled down by 7.55%
 174     fResponseMatrixMinus10(), // Response matrix with analysed multiplicity scaled down by 10%
 175     fResponseMatrixMinusSys(),// Response matrix with analysed multiplicity scaled up by event mult uncertainty
 176     fResponseMatrixPlusSys(), // Response matrix with analysed multiplicity scaled down by event mult uncertainty
 177     fESDNSD(),                // number of events found as NSD by the analysis vs. multiplicity
 178     fMCNSD(),                 // number of events found as NSD by the MC truth vs. multiplicity
 179     fMCESDNSD(),              // number of events found as NSD by both analysis and MC truth vs. multiplicity
 180     fTriggerBias()             // histogram for trigger vertex bias correction
 181 {
 182   //
 183   // Constructor
 184   //
 185   const char* name = AliForwardMultiplicityDistribution::FormBinName(etaLow,etaHigh);
 186
 187   SetName(name);
 188   SetTitle(Form("%+4.1f < #eta < %+4.1f", fEtaLow, fEtaHigh));
 189
 190 }
 191 //_____________________________________________________________________
 192 AliForwardCreateResponseMatrices::Bin::Bin()
 193   : TNamed(),
 194     fEtaLow(),          // low eta limit
 195     fEtaHigh(),        // high eta limit
 196     fHist(),                  // multiplicity histogram
 197     fHistMC(),                // multiplicity histogram MC truth primaries
 198     fAcceptance(),            // histogram showing the 'holes' in acceptance. BinContent of 1 shows a hole, and BinContent of 10 shows data coverage
 199     fVtxZvsNdataBins(),       // VtxZ vs. number of data acceptance bins (normalised to the eta range)
 200     fResponseMatrix(),        // Response matrix (MC truth vs. analysed multiplicity)
 201     fResponseMatrixPlus05(),  // Response matrix with analysed multiplicity scaled up by 5%
 202     fResponseMatrixPlus075(), // Response matrix  with analysed multiplicity scaled up by 7.5%
 203     fResponseMatrixPlus10(),  // Response matrix with analysed multiplicity scaled up by 10%
 204     fResponseMatrixMinus05(), // Response matrix with analysed multiplicity scaled down by 5%
 205     fResponseMatrixMinus075(),// Response matrix with analysed multiplicity scaled down by 7.55%
 206     fResponseMatrixMinus10(), // Response matrix with analysed multiplicity scaled down by 10%
 207     fResponseMatrixMinusSys(),// Response matrix with analysed multiplicity scaled up by event mult uncertainty
 208     fResponseMatrixPlusSys(), // Response matrix with analysed multiplicity scaled down by event mult uncertainty
 209     fESDNSD(),                // number of events found as NSD by the analysis vs. multiplicity
 210     fMCNSD(),                 // number of events found as NSD by the MC truth vs. multiplicity
 211     fMCESDNSD(),              // number of events found as NSD by both analysis and MC truth vs. multiplicity
 212     fTriggerBias()             // histogram for trigger vertex bias correction
 213 {
 214   //
 215   // Default Constructor
 216   //
 217 }
 218 //_____________________________________________________________________
 219 void AliForwardCreateResponseMatrices::Bin::CreateOutputObjectss(TList* cont,  Int_t max)
 220 {
 221   //
 222   // Define eta bin output histos
 223   //
 224   TList* out = new TList;
 225   out->SetName(GetName());
 226   cont->Add(out);
 227
 228   fHist                    = new TH1D("mult", GetTitle(), max, -0.5, max-.5);
 229   fHistMC                  = new TH1D("multTruth", GetTitle(), max, -0.5, max-.5);
 230   fVtxZvsNdataBins         = new TH2D("VtxZvsNdataBins", "VtxZ vs dataAcceptance/etaRange;z-vtz;dataAcceptance/etaRange", 20, -10,10, 130,0,1.3);
 231   fAcceptance              = new TH2D("Acceptance","Acceptance;#eta;z-vtx",200,-4, 6 , 20,-10,10);
 232   fResponseMatrix          = new TH2D("responseMatrix","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
 233   fResponseMatrixPlus05    = new TH2D("responseMatrixPlus05","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
 234   fResponseMatrixPlus075   = new TH2D("responseMatrixPlus075","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
 235   fResponseMatrixPlus10    = new TH2D("responseMatrixPlus10","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
 236   fResponseMatrixMinus05   = new TH2D("responseMatrixMinus05","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
 237   fResponseMatrixMinus075  = new TH2D("responseMatrixMinus075","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
 238   fResponseMatrixMinus10   = new TH2D("responseMatrixMinus10","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
 239   fResponseMatrixMinusSys  = new TH2D("responseMatrixMinusSys","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
 240   fResponseMatrixPlusSys   = new TH2D("responseMatrixPlusSys","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
 241   fTriggerBias             = new TH1D("triggerBias","triggerBias;MC_{truth};Correction}", max, -0.5, max-.5);
 242   fMCNSD                   = new TH1D("fMCNSD","fMCNSD", 300,-0.5,299.5);
 243   fESDNSD                  = new TH1D("fESDNSD","fESDNSD", 300,-0.5,299.5);
 244   fMCESDNSD                = new TH1D("fMCESDNSD","fMCESDNSD", 300,-0.5,299.5);
 245
 246   out->Add(fMCNSD);
 247   out->Add(fESDNSD);
 248   out->Add(fMCESDNSD);
 249   out->Add(fHist);
 250   out->Add(fHistMC);
 251   out->Add(fVtxZvsNdataBins);
 252   out->Add(fAcceptance);
 253   out->Add(fResponseMatrix);
 254   out->Add(fResponseMatrixPlus05);
 255   out->Add(fResponseMatrixPlus075);
 256   out->Add(fResponseMatrixPlus10);
 257   out->Add(fResponseMatrixMinus05);
 258   out->Add(fResponseMatrixMinus075);
 259   out->Add(fResponseMatrixMinus10);
 260   out->Add(fResponseMatrixPlusSys);
 261   out->Add(fResponseMatrixMinusSys);
 262   out->Add(fTriggerBias);
 263
 264 }
 265
 266
 267 //_____________________________________________________________________
 268 void
 269 AliForwardCreateResponseMatrices::Bin::Process(TH1D* dndetaForward,
 270                                                TH1D* dndetaCentral,
 271                                                TH1D* normForward,
 272                                                TH1D* normCentral,
 273                                                TH1D* mc,
 274                                                Double_t VtxZ,
 275                                                Bool_t selectedTrigger,
 276                                                Bool_t isMCNSD,
 277                                                Bool_t isESDNSD,
 278                                                const AliAODEvent& aodevent)
 279 {
 280   //
 281   // Process a single eta bin
 282   //
 283
 284   // Diagnostics on event acceptance
 285   Int_t    first = dndetaForward->GetXaxis()->FindBin(fEtaLow);
 286   Int_t    last  = dndetaForward->GetXaxis()->FindBin(fEtaHigh-.0001);
 287   Double_t acceptanceBins=0;
 288   for(Int_t n= first;n<=last;n++){
 289     if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0){
 290       acceptanceBins++;
 291     }
 292     fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ), 1);
 293     if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0)
 294       fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ),10);
 295   }
 296   fVtxZvsNdataBins->Fill(VtxZ, (Double_t)acceptanceBins/(last-first+1));
 297
 298
 299
 300   Double_t c        = 0;
 301   Double_t e2       = 0;
 302   Double_t cPrimary = 0;
 303   Double_t e2Primary= 0;
 304
 305   for (Int_t i = first; i <= last; i++){
 306     Double_t cForward  = 0;
 307     Double_t cCentral  = 0;
 308     Double_t e2Forward = 0;
 309     Double_t e2Central = 0;
 310     Double_t cMC  = 0;
 311     Double_t e2MC = 0;
 312     cMC= mc->GetBinContent(i);
 313     e2MC= mc->GetBinError(i) * mc->GetBinError(i);
 314     cPrimary+=cMC;
 315     e2Primary+=e2MC;
 316     if (normForward->GetBinContent(i) > 0) {
 317       cForward  = dndetaForward->GetBinContent(i);
 318       e2Forward = dndetaForward->GetBinError(i) * dndetaForward->GetBinError(i);
 319     }
 320     if (normCentral->GetBinContent(i) > 0) {
 321       cCentral  = dndetaCentral->GetBinContent(i);
 322       e2Central = dndetaCentral->GetBinError(i) * dndetaCentral->GetBinError(i);
 323     }
 324     Double_t cc  = 0;
 325     Double_t ee2 = 0;
 326     if (cCentral > 0 && cForward > 0) {
 327       cc  = 0.5 * (cForward  + cCentral);
 328       ee2 = 0.5 * (e2Forward + e2Central);
 329     }
 330     else if (cCentral > 0) {
 331       cc  = cCentral;
 332       ee2 = e2Central;
 333     }
 334     else if (cForward > 0) {
 335       cc  = cForward;
 336       ee2 = e2Forward;
 337     }
 338     c  += cc;
 339     e2 += ee2;
 340   }
 341
 342   // retreive MC particles from event
 343   TClonesArray* mcArray = (TClonesArray*)aodevent.FindListObject(AliAODMCParticle::StdBranchName());
 344   if(!mcArray){
 345     AliWarning("No MC array found in AOD. Try making it again.");
 346     return;
 347   }
 348   AliAODMCHeader* header =
 349     dynamic_cast<AliAODMCHeader*>(aodevent.
 350                                   FindListObject(AliAODMCHeader::StdBranchName()));
 351   if (!header) {
 352     AliWarning("No header file found.");
 353     return;
 354   }
 355
 356
 357   Int_t ntracks = mcArray->GetEntriesFast();
 358   // Track loop: find MC truth multiplicity in selected eta bin
 359   Double_t mult = 0;
 360   for (Int_t it = 0; it < ntracks; it++) {
 361     AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
 362     if (!particle) {
 363       AliError(Form("Could not receive track %d", it));
 364       continue;
 365     }
 366     if (!particle->IsPhysicalPrimary()) continue;
 367     if (particle->Charge() == 0) continue;
 368     if(particle->Eta()>fEtaLow&&particle->Eta()<fEtaHigh-0.0001)
 369       mult++;
 370   }
 371   //fill fMCNSD with multiplicity of MC truth NSD events with MC truth |vtxz|<4
 372   if(header->GetVtxZ()>-4&&header->GetVtxZ()<4){
 373     if(isMCNSD){
 374       fMCNSD->Fill(mult);
 375     }
 376   }
 377   //fill fESDNSD with multiplicity from events with a reconstructed NSD trigger and reconstructed |vtxz|<4
 378   if(VtxZ>-4&&VtxZ<4){
 379     if(isESDNSD){
 380       fESDNSD->Fill(mult);
 381     }
 382   }
 383   // fullfilling both requirements of fMCNSD and fESDNSD
 384   if(header->GetVtxZ()>-4&&header->GetVtxZ()<4&&VtxZ>-4&&VtxZ<4&&isMCNSD&&isESDNSD){
 385     fMCESDNSD->Fill(mult);
 386   }
 387
 388
 389
 390 //Systematic errors from here
 391
 392   Double_t fmd=0;
 393   Double_t spd=0;
 394   Double_t overlap=0;
 395
 396   // number of eta bins in fmd, spd and overlap respectively
 397   for(Int_t i = first;i<=last;i++){
 398     if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)<1)
 399       fmd++;
 400     if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)>0)
 401       overlap++;
 402     if(normCentral->GetBinContent(i)>0&&normForward->GetBinContent(i)<1)
 403       spd++;
 404   }
 405
 406   Double_t sysErrorSquared = 0;
 407
 408   // estimate of systematic uncertainty on the event multiplicity - hardcoded :(. estimates taken from Hans Hjersing Dalsgaard or Casper Nygaard phd theses.
 409   Double_t fmdSysError= 0.08;
 410   Double_t spdSysError= 0.04;
 411   Double_t total = 0;
 412   total= fmd + spd + overlap;
 413   if(total){
 414     // Combined systematc event uncertainty, by weighting with the number of eta-bins of fmd-only, spd-only and the overlap.
 415     sysErrorSquared= (fmd*TMath::Power(fmdSysError,2)+ spd*TMath::Power(spdSysError,2)+
 416                       0.5*overlap*TMath::Power(fmdSysError,2)+ 0.5*overlap*TMath::Power(spdSysError,2))/total;
 417   }
 418
 419
 420   if(selectedTrigger){
 421     fHist->Fill(c);
 422     fHistMC->Fill(cPrimary);
 423     fResponseMatrix->Fill(cPrimary, c);
 424     fResponseMatrixPlusSys->Fill(cPrimary,(1+TMath::Sqrt(sysErrorSquared))*c);
 425     fResponseMatrixMinusSys->Fill(cPrimary,(1-TMath::Sqrt(sysErrorSquared))*c);
 426     fResponseMatrixPlus05->Fill(cPrimary, 1.05*c);
 427     fResponseMatrixPlus075->Fill(cPrimary, 1.075*c);
 428     fResponseMatrixPlus10->Fill(cPrimary, 1.1*c);
 429     fResponseMatrixMinus05->Fill(cPrimary, 0.95*c);
 430     fResponseMatrixMinus075->Fill(cPrimary, 0.925*c);
 431     fResponseMatrixMinus10->Fill(cPrimary, 0.9*c);
 432
 433   }
 434
 435 }
 436
 437
 438
 439
 440 //_____________________________________________________________________
 441 //
 442 //
 443 // EOF