PWGLF/FORWARD/analysis2/AliForwardCreateResponseMatrices.cxx

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