]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ANALYSIS/AliCollisionNormalization.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Small corrections to the TENDER code, which will support AliTenderInputHander(current...
[u/mrichter/AliRoot.git] / ANALYSIS / AliCollisionNormalization.cxx
... / ...
CommitLineData
1//-------------------------------------------------------------------------
2// Implementation of Class AliCollisionNormalization
3//
4// This class is used to store the vertex ditributions in the data
5// and in Monte Carlo, needed to compute the real number of
6// collisions a given sample is corresponding to.
7// The strategy matches what described in CERN-THESIS-2009-033 p 119.
8//
9// Author: Michele Floris, CERN
10//-------------------------------------------------------------------------
11
12#include "AliCollisionNormalization.h"
13#include "AliPhysicsSelection.h"
14#include "AliLog.h"
15#include "TFile.h"
16#include "TCanvas.h"
17#include "AliGenPythiaEventHeader.h"
18#include "AliGenDPMjetEventHeader.h"
19#include "AliGenEventHeader.h"
20#include "AliMCEvent.h"
21
22ClassImp(AliCollisionNormalization)
23
24const char * AliCollisionNormalization::fProcLabel[] = {"SD","DD","ND", "Unknown"};
25
26AliCollisionNormalization::AliCollisionNormalization() :
27 TObject(),
28 fNbinsVz(0),
29 fMinVz(0),
30 fMaxVz(0),
31 fZRange(9.99),
32 fIsMC(0),
33 fReferenceXS(0),
34 fVerbose(0),
35 fEnergy(900),
36 fHistVzData (0),
37 fHistProcTypes (0),
38 fHistStatBin0 (0),
39 fHistStat (0)
40{
41
42 // ctor
43
44 fNbinsVz = 30;
45 fMinVz = -15;
46 fMaxVz = +15;
47
48 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
49 fHistVzMCGen[iproc] = 0;
50 fHistVzMCRec[iproc] = 0;
51 fHistVzMCTrg[iproc] = 0;
52 }
53
54
55 BookAllHistos();
56}
57
58AliCollisionNormalization::AliCollisionNormalization(Int_t nbinz, Float_t minz, Float_t maxz):
59 TObject(),
60 fNbinsVz(0),
61 fMinVz(0),
62 fMaxVz(0),
63 fZRange(9.99),
64 fIsMC(0),
65 fReferenceXS(0),
66 fVerbose(0),
67 fEnergy(900),
68 fHistVzData (0),
69 fHistProcTypes (0),
70 fHistStatBin0 (0),
71 fHistStat (0)
72
73{
74
75 // ctor, allows setting binning
76 fNbinsVz = nbinz;
77 fMinVz = minz;
78 fMaxVz = maxz;
79
80 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
81 fHistVzMCGen[iproc] = 0;
82 fHistVzMCRec[iproc] = 0;
83 fHistVzMCTrg[iproc] = 0;
84 }
85
86 BookAllHistos();
87}
88
89AliCollisionNormalization::AliCollisionNormalization(const char * dataFile, const char * dataListName,
90 const char * mcFile, const char * mcListName,
91 const char * eventStatFile) :
92 TObject(),
93 fNbinsVz(0),
94 fMinVz(0),
95 fMaxVz(0),
96 fZRange(9.99),
97 fIsMC(0),
98 fReferenceXS(0),
99 fVerbose(0),
100 fEnergy(900),
101 fHistVzData (0),
102 fHistProcTypes (0),
103 fHistStatBin0 (0),
104 fHistStat (0)
105
106{
107
108 // ctor, loads histograms from file
109 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
110 fHistVzMCGen[iproc] = 0;
111 fHistVzMCRec[iproc] = 0;
112 fHistVzMCTrg[iproc] = 0;
113 }
114
115
116 TFile * fdata = new TFile (dataFile);
117 TFile * fmc = new TFile (mcFile );
118 TFile * fstat = new TFile(eventStatFile);
119
120 if (!fdata->IsOpen() || !fmc->IsOpen() || !fstat->IsOpen()) {
121 AliFatal("Cannot open input file(s)");
122 }
123
124 TList * ldata = (TList*) fdata->Get(dataListName);
125 TList * lmc = (TList*) fmc ->Get(mcListName );
126
127 AliCollisionNormalization * cndata = (AliCollisionNormalization*) ldata->FindObject("AliCollisionNormalization");
128 AliCollisionNormalization * cnmc = (AliCollisionNormalization*) lmc ->FindObject("AliCollisionNormalization");
129
130
131 // Assign or book all histos
132 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
133 fHistVzMCGen[iproc]= cnmc->GetVzMCGen(iproc);
134 fHistVzMCRec[iproc]= cnmc->GetVzMCRec(iproc);
135 fHistVzMCTrg[iproc]= cnmc->GetVzMCTrg(iproc);
136 }
137 fHistVzData = cndata->GetVzData();
138 fHistProcTypes = cnmc->GetHistProcTypes();
139
140 fHistStatBin0 = (TH1F*) fstat->Get("fHistStatistics_Bin0");
141 fHistStat = (TH1F*) fstat->Get("fHistStatistics");
142
143}
144
145
146AliCollisionNormalization::~AliCollisionNormalization(){
147
148 // dtor
149 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
150 if(fHistVzMCGen[iproc]) { delete fHistVzMCGen[iproc] ; fHistVzMCGen[iproc] =0;}
151 if(fHistVzMCRec[iproc]) { delete fHistVzMCRec[iproc] ; fHistVzMCRec[iproc] =0;}
152 if(fHistVzMCTrg[iproc]) { delete fHistVzMCTrg[iproc] ; fHistVzMCTrg[iproc] =0;}
153 }
154
155 if(fHistVzData ) { delete fHistVzData ; fHistVzData =0;}
156 if(fHistStatBin0 ) { delete fHistStatBin0 ; fHistStatBin0 =0;}
157 if(fHistStat ) { delete fHistStat ; fHistStat =0;}
158 if(fHistProcTypes ) { delete fHistProcTypes ; fHistProcTypes =0;}
159
160}
161
162void AliCollisionNormalization::BookAllHistos(){
163 // books all histos
164 // Book histos of vz distributions vs multiplicity
165 // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
166
167 // Do not attach histos to the current directory
168 Bool_t oldStatus = TH1::AddDirectoryStatus();
169 TH1::AddDirectory(kFALSE);
170
171 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
172 fHistVzMCGen [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCGen")+ fProcLabel[iproc] ,"Vz distribution of generated events vs rec multiplicity ");
173 fHistVzMCRec [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fProcLabel[iproc] ,"Vz distribution of reconstructed events vs rec multiplicity");
174 fHistVzMCTrg [iproc] = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fProcLabel[iproc] ,"Vz distribution of triggered events vs rec multiplicity ");
175 }
176 fHistVzData = (TH2F*) BookVzHisto("fHistVzData" ,"Vz distribution of triggered events vs rec multiplicity ");
177 fHistProcTypes = new TH1F ("fHistProcTypes", "Number of events in the different process classes", kNProcs, -0.5 , kNProcs-0.5);
178
179 fHistProcTypes->GetXaxis()->SetBinLabel(kProcSD+1,"SD");
180 fHistProcTypes->GetXaxis()->SetBinLabel(kProcND+1,"ND");
181 fHistProcTypes->GetXaxis()->SetBinLabel(kProcDD+1,"DD");
182 fHistProcTypes->GetXaxis()->SetBinLabel(kProcUnknown+1,"Unknown");
183
184 TH1::AddDirectory(oldStatus);
185
186}
187
188TH1 * AliCollisionNormalization::BookVzHisto(const char * name , const char * title, Bool_t vzOnly) {
189
190 // Book histos of vz distributions vs multiplicity
191 // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
192
193
194 // Do not attach histos to the current directory
195 Bool_t oldStatus = TH1::AddDirectoryStatus();
196 TH1::AddDirectory(kFALSE);
197
198 TH1 * h;
199 Double_t binLimitsVtx[] = {-30,-25,-20,-15,-10,-7,-5.5,-4,-3,-2,-1,0,1,2,3,4,5.5,7,10,15,20,25,30};
200 if (vzOnly) {
201 h = new TH1F(name,title,22,binLimitsVtx);
202 } else {
203 h = new TH2F(name,title,22,binLimitsVtx,100,Double_t(-0.5),Double_t(99.5));
204 }
205
206 h->SetXTitle("V_{z} (cm)");
207 h->SetYTitle("n_{trk}");
208 h->Sumw2();
209
210 TH1::AddDirectory(oldStatus);
211
212 return h;
213
214}
215
216TH2F * AliCollisionNormalization::GetVzMCGen (Int_t procType) {
217
218 // returns MC gen histo. If proc type < 0 sums over all proc types, reweighting the XS
219
220 if(procType>=0) return fHistVzMCGen[procType] ;
221
222 TH2F * sum = (TH2F*) fHistVzMCGen[kProcSD]->Clone();
223 sum->Reset();
224
225 for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
226 sum->Add(fHistVzMCGen[iproc],GetProcessWeight(iproc));
227 }
228
229 return sum;
230}
231TH2F * AliCollisionNormalization::GetVzMCRec (Int_t procType) {
232
233 // returns MC rec histo. If proc type < 0 sums over all proc types, reweighting the XS
234
235 if(procType>=0) return fHistVzMCRec[procType] ;
236
237 TH2F * sum = (TH2F*) fHistVzMCRec[kProcSD]->Clone();
238 sum->Reset();
239
240 for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
241 sum->Add(fHistVzMCRec[iproc],GetProcessWeight(iproc));
242 }
243
244 return sum;
245
246}
247
248
249TH2F * AliCollisionNormalization::GetVzMCTrg (Int_t procType) {
250
251 // returns MC trg histo. If proc type < 0 sums over all proc types, reweighting the XS
252
253 if(procType>=0) return fHistVzMCTrg[procType] ;
254
255 TH2F * sum = (TH2F*) fHistVzMCTrg[kProcSD]->Clone();
256 sum->Reset();
257
258 for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
259 sum->Add(fHistVzMCTrg[iproc],GetProcessWeight(iproc));
260 }
261
262 return sum;
263
264}
265
266Double_t AliCollisionNormalization::ComputeNint() {
267
268 // Compute the number of collisions applying all corrections
269 // TODO: check error propagation
270
271 TH1 * hVzData = fHistVzData->ProjectionX("_px",2,-1,"E"); // Skip zero bin
272 Int_t allEventsWithVertex = (Int_t) fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,
273 2, fHistVzData->GetNbinsY()+1); // include under/overflow!, skip 0 bin
274
275 // Assign histos reweighted with measured XS
276 TH2F * histVzMCGen = GetVzMCGen(-1);
277 TH2F * histVzMCTrg = GetVzMCTrg(-1);
278 TH2F * histVzMCRec = GetVzMCRec(-1);
279
280 // Before we start: print number of input events to the physics
281 // selection: this allows to crosscheck that all runs were
282 // successfully processed (useful if running on grid: you may have a
283 // crash without noticing it).
284 AliInfo(Form("Input Events (No cuts: %d, After Phys. Sel.:%d)",
285 Int_t(fHistStat->GetBinContent(1,1)),
286 Int_t(fHistStat->GetBinContent(fHistStat->GetNbinsX(),1)))); // Fixme: will this change with a different trigger scheme?
287
288 // Get or compute BG. This assumes the CINT1B suite
289 Double_t triggeredEventsWith0MultWithBG = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,1, 1);
290 Double_t bg = 0; // This will include beam gas + accidentals
291 if (fHistStatBin0->GetNbinsY() > 4) { // FIXME: we need a better criterion to decide...
292 AliInfo("Using BG computed by Physics Selection");
293 // WARNING
294 // CHECK IF THE COMPUTATION OF BG OFFSET IS STILL OK, IN CASE OF CHANGES TO THE PHYSICS SELECTION
295 Int_t bgOffset = 0;
296 Int_t nbiny = fHistStatBin0->GetNbinsY();
297 for(Int_t ibiny =1; ibiny <= nbiny; ibiny++){
298 bgOffset++;
299 printf("%d, %s\n", bgOffset, fHistStatBin0->GetYaxis()->GetBinLabel(ibiny) );
300
301 if(!strncmp("All B", fHistStatBin0->GetYaxis()->GetBinLabel(ibiny),5)) break;
302 if((ibiny+1)==nbiny) AliFatal("Cannot compute bg offset");
303 }
304
305 if(fVerbose > 2) AliInfo(Form("BG Offset: %d",bgOffset));
306
307
308 bg = fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(), bgOffset+AliPhysicsSelection::kStatRowBG);
309 bg += fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),bgOffset+AliPhysicsSelection::kStatRowAcc);
310 Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),1);
311 if (cint1B != Int_t(triggeredEventsWith0MultWithBG)) {
312 AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, Int_t(triggeredEventsWith0MultWithBG)));
313 }
314 } else {
315 AliInfo("Computing BG using CINT1A/B/C/E, ignoring intensities");
316 AliError("This will only work for early runs!!! USE BG FROM PHYSICS SELECTION INSTEAD");
317 Int_t icol = fHistStatBin0->GetNbinsX();
318 Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(icol,1);
319 Int_t cint1A = (Int_t) fHistStatBin0->GetBinContent(icol,2);
320 Int_t cint1C = (Int_t) fHistStatBin0->GetBinContent(icol,3);
321 Int_t cint1E = (Int_t) fHistStatBin0->GetBinContent(icol,4);
322 bg = cint1A + cint1C-2*cint1E ; // FIXME: to be changed to take into account ratios of events
323 if (cint1B != triggeredEventsWith0MultWithBG) {
324 AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, (Int_t)triggeredEventsWith0MultWithBG));
325 }
326 }
327
328 Double_t triggeredEventsWith0Mult = triggeredEventsWith0MultWithBG - bg;
329 if(fVerbose > 0) AliInfo(Form("Measured events with vertex: %d",allEventsWithVertex));
330 Double_t bin0 = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,
331 1, 1);
332 if(fVerbose > 0) AliInfo(Form("Zero Bin, Meas: %2.2f, BG: %2.2f, Meas - BG: %2.2f", bin0, bg, triggeredEventsWith0Mult));
333
334 // This pointers are here so that I use the same names used in Jan Fiete's code (allows for faster/easier comparison)
335
336 TH2* eTrig = histVzMCTrg;
337 TH2* eTrigVtx = histVzMCRec;
338 TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1);
339
340 // compute trigger and vertex efficiency
341 TH2 * effVtxTrig = (TH2*) histVzMCRec->Clone("effVtxTrig");
342 effVtxTrig->Reset();
343 effVtxTrig->Divide(histVzMCRec,histVzMCGen,1,1,"B");
344 // Apply correction to data to get corrected events
345 TH2 * correctedEvents = (TH2*) fHistVzData->Clone("correctedEvents");
346 correctedEvents->Divide(effVtxTrig);
347
348 // TH1 * correctedEvents = fHistVzData->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1);
349
350 // loop over vertex bins
351 for (Int_t i = 1; i <= fHistVzData->GetNbinsX(); i++)
352 {
353 Double_t alpha = (Double_t) hVzData->GetBinContent(i) / allEventsWithVertex;
354 Double_t events = alpha * triggeredEventsWith0Mult;
355
356 // if (histVzMCRec->GetBinContent(i,1) == 0)
357 // continue;
358 if (!eTrigVtx_projx->GetBinContent(i) || !eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1))
359 continue;
360
361 Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) *
362 eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
363
364 events *= fZ;
365
366 // multiply with trigger correction
367 Double_t correction = histVzMCGen->GetBinContent(i,1)/histVzMCTrg->GetBinContent(i,1);
368 events *= correction;
369
370 if (fVerbose > 1) Printf(" Bin %d, alpha is %.2f%%, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events,
371 histVzMCTrg->GetBinContent(i,1));
372 correctedEvents->SetBinContent(i, 1, events);
373 }
374
375
376
377 // Integrate correctedEvents over full z range
378 Double_t allEvents = correctedEvents->Integral(0, correctedEvents->GetNbinsX()+1,0, correctedEvents->GetNbinsY()+1);
379 // Integrate correctedEvents over needed z range
380 Double_t allEventsZrange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),
381 0, correctedEvents->GetNbinsY()+1);
382
383 if(fVerbose > 1) AliInfo(Form("Results in |Vz| < %3.3f",fZRange));
384 if(fVerbose > 1) AliInfo(Form(" Events in Bin0: %2.2f, With > 1 track: %2.2f, All corrected: %2.2f",
385 correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),1,1),
386 correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),
387 2,correctedEvents->GetNbinsX()+1),
388 allEventsZrange));
389 if(fVerbose > 1) {
390 Int_t nbin = histVzMCTrg->GetNbinsX();
391 AliInfo(Form("Trigger Efficiency in the zero bin: %3.3f", histVzMCTrg->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1) ));
392 }
393
394
395 AliInfo(Form("Number of collisions in full phase space: %2.2f", allEvents));
396// effVtxTrig->Draw();
397// new TCanvas();
398// correctedEvents->Draw(); // FIXME: debug
399
400 return allEvents;
401}
402
403Long64_t AliCollisionNormalization::Merge(TCollection* list)
404{
405 // Merge a list of AliCollisionNormalization objects with this
406 // (needed for PROOF).
407 // Returns the number of merged objects (including this).
408
409 if (!list)
410 return 0;
411
412 if (list->IsEmpty())
413 return 1;
414
415 TIterator* iter = list->MakeIterator();
416 TObject* obj;
417
418 // collections of all histograms
419 const Int_t nHists = kNProcs*3+5;
420 TList collections[nHists];
421
422 Int_t count = 0;
423 while ((obj = iter->Next())) {
424
425 AliCollisionNormalization* entry = dynamic_cast<AliCollisionNormalization*> (obj);
426 if (entry == 0)
427 continue;
428 Int_t ihist = -1;
429
430 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
431 if (entry->fHistVzMCGen[iproc] ) collections[++ihist].Add(entry->fHistVzMCGen[iproc] );
432 if (entry->fHistVzMCRec[iproc] ) collections[++ihist].Add(entry->fHistVzMCRec[iproc] );
433 if (entry->fHistVzMCTrg[iproc] ) collections[++ihist].Add(entry->fHistVzMCTrg[iproc] );
434 }
435 if (entry->fHistVzData ) collections[++ihist].Add(entry->fHistVzData );
436 if (entry->fHistProcTypes ) collections[++ihist].Add(entry->fHistProcTypes );
437 if (entry->fHistStatBin0 ) collections[++ihist].Add(entry->fHistStatBin0 );
438 if (entry->fHistStat ) collections[++ihist].Add(entry->fHistStat );
439
440 count++;
441 }
442
443 Int_t ihist = -1;
444 for(Int_t iproc = 0; iproc < kNProcs; iproc++){
445 if (fHistVzMCGen[iproc] ) fHistVzMCGen[iproc] ->Merge(&collections[++ihist]);
446 if (fHistVzMCRec[iproc] ) fHistVzMCRec[iproc] ->Merge(&collections[++ihist]);
447 if (fHistVzMCTrg[iproc] ) fHistVzMCTrg[iproc] ->Merge(&collections[++ihist]);
448
449 }
450 if (fHistVzData ) fHistVzData ->Merge(&collections[++ihist]);
451 if (fHistProcTypes ) fHistProcTypes ->Merge(&collections[++ihist]);
452 if (fHistStatBin0 ) fHistStatBin0 ->Merge(&collections[++ihist]);
453 if (fHistStat ) fHistStat ->Merge(&collections[++ihist]);
454
455
456 delete iter;
457
458 return count+1;
459}
460
461void AliCollisionNormalization::FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {
462
463 // Fill MC gen histo and the process types statistics
464
465 Int_t evtype = GetProcessType(mcEvt);
466 // When I fill gen histos, I also fill statistics of process types (used to detemine ratios afterwards).
467 fHistProcTypes->Fill(evtype);
468 fHistVzMCGen[evtype]->Fill(vz,ntrk);
469}
470
471void AliCollisionNormalization::FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt){
472
473 // Fill MC rec histo
474 Int_t evtype = GetProcessType(mcEvt);
475 fHistVzMCRec[evtype]->Fill(vz,ntrk);
476
477}
478void AliCollisionNormalization::FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {
479
480 // Fill MC trg histo
481 Int_t evtype = GetProcessType(mcEvt);
482 fHistVzMCTrg[evtype]->Fill(vz,ntrk);
483}
484
485
486Int_t AliCollisionNormalization::GetProcessType(AliMCEvent * mcEvt) {
487
488 // Determine if the event was generated with pythia or phojet and return the process type
489
490 // Check if mcEvt is fine
491 if (!mcEvt) {
492 AliFatal("NULL mc event");
493 }
494
495 // Determine if it was a pythia or phojet header, and return the correct process type
496 AliGenPythiaEventHeader * headPy = 0;
497 AliGenDPMjetEventHeader * headPho = 0;
498 AliGenEventHeader * htmp = mcEvt->GenEventHeader();
499 if(!htmp) {
500 AliFatal("Cannot Get MC Header!!");
501 }
502 if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
503 headPy = (AliGenPythiaEventHeader*) htmp;
504 } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
505 headPho = (AliGenDPMjetEventHeader*) htmp;
506 } else {
507 AliError("Unknown header");
508 }
509
510 // Determine process type
511 if(headPy) {
512 if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
513 // single difractive
514 return kProcSD;
515 } else if (headPy->ProcessType() == 94) {
516 // double diffractive
517 return kProcDD;
518 }
519 else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {
520 // non difractive
521 return kProcND;
522 }
523 } else if (headPho) {
524 if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
525 // single difractive
526 return kProcSD;
527 } else if (headPho->ProcessType() == 7) {
528 // double diffractive
529 return kProcDD;
530 } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6 && headPho->ProcessType() != 7 ) {
531 // non difractive
532 return kProcND;
533 }
534 }
535
536
537 // no process type found?
538 AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
539 return kProcUnknown;
540}
541
542Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) {
543
544 // Return a weight to be used when combining the MC histos to
545 // compute efficiency, defined as the ratio XS in generator / XS
546 // measured
547
548 Float_t ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND;
549 GetRelativeFractions(fReferenceXS,ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
550
551 static Double_t total = fHistProcTypes->Integral();
552 if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) {
553 AliError("There were unknown processes!!!");
554 }
555 static Double_t SD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total;
556 static Double_t DD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total;
557 static Double_t ND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total;
558
559 if (fVerbose > 2) {
560 AliInfo(Form("Total MC evts: %f",total));
561 AliInfo(Form(" Frac SD %4.4f", SD));
562 AliInfo(Form(" Frac DD %4.4f", DD));
563 AliInfo(Form(" Frac ND %4.4f", ND));
564 }
565
566 switch(proc) {
567 case kProcSD:
568 return ref_SD/SD;
569 break;
570 case kProcDD:
571 return ref_DD/DD;
572 break;
573 case kProcND:
574 return ref_ND/ND;
575 break;
576 default:
577 AliError("Unknown process");
578 }
579
580 return 0;
581
582}
583
584void AliCollisionNormalization::GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND)
585{
586 // Returns fraction of XS (SD, ND, DD) and corresponding error
587 // Stolen from Jan Fiete's drawSystematics macro
588
589 // origin:
590 // -1 = Pythia (test)
591 // 0 = UA5
592 // 1 = Data 1.8 TeV
593 // 2 = Tel-Aviv
594 // 3 = Durham
595 //
596
597 Double_t epsilon = 0.0001;
598 if(TMath::Abs(fEnergy-900)<epsilon) {
599
600 switch (origin)
601 {
602 case -10: // Pythia default at 7 GeV, 50% error
603 AliInfo("PYTHIA x-sections");
604 ref_SD = 0.192637; error_SD = ref_SD * 0.5;
605 ref_DD = 0.129877; error_DD = ref_DD * 0.5;
606 ref_ND = 0.677486; error_ND = 0;
607 break;
608
609 case -1: // Pythia default at 900 GeV, as test
610 AliInfo("PYTHIA x-sections");
611 ref_SD = 0.223788;
612 ref_DD = 0.123315;
613 ref_ND = 0.652897;
614 break;
615
616 case 0: // UA5
617 AliInfo("UA5 x-sections a la first paper");
618 ref_SD = 0.153; error_SD = 0.05;
619 ref_DD = 0.080; error_DD = 0.05;
620 ref_ND = 0.767; error_ND = 0;
621 break;
622
623 case 10: // UA5
624 AliInfo("UA5 x-sections hadron level definition for Pythia");
625 // Fractions in Pythia with UA5 cuts selection for SD
626 // ND: 0.688662
627 // SD: 0.188588 --> this should be 15.3
628 // DD: 0.122750
629 ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189;
630 ref_DD = 0.095; error_DD = 0.06;
631 ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
632 break;
633
634 case 11: // UA5
635 AliInfo("UA5 x-sections hadron level definition for Phojet");
636 // Fractions in Phojet with UA5 cuts selection for SD
637 // ND: 0.783573
638 // SD: 0.151601 --> this should be 15.3
639 // DD: 0.064827
640 ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152;
641 ref_DD = 0.095; error_DD = 0.06;
642 ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
643 break;
644 case 2: // tel-aviv model
645 AliInfo("Tel-aviv model x-sections");
646 ref_SD = 0.171;
647 ref_DD = 0.094;
648 ref_ND = 1 - ref_SD - ref_DD;
649 break;
650
651 case 3: // durham model
652 AliInfo("Durham model x-sections");
653 ref_SD = 0.190;
654 ref_DD = 0.125;
655 ref_ND = 1 - ref_SD - ref_DD;
656 break;
657 default:
658 AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy));
659 }
660 }
661 else if(TMath::Abs(fEnergy-1800)<epsilon) {
662 switch (origin)
663 {
664 case 20: // E710, 1.8 TeV
665 AliInfo("E710 x-sections hadron level definition for Pythia");
666 // ND: 0.705709
667 // SD: 0.166590 --> this should be 15.9
668 // DD: 0.127701
669 ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167;
670 ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43;
671 ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
672 break;
673
674 case 21: // E710, 1.8 TeV
675 AliInfo("E710 x-sections hadron level definition for Phojet");
676 // ND: 0.817462
677 // SD: 0.125506 --> this should be 15.9
678 // DD: 0.057032
679 ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126;
680 ref_DD = 0.075 * 1.43; error_DD = 0.02 * 1.43;
681 ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
682 break;
683
684 case 1: // data 1.8 TeV
685 AliInfo("??? x-sections");
686 ref_SD = 0.152;
687 ref_DD = 0.092;
688 ref_ND = 1 - ref_SD - ref_DD;
689 break;
690 default:
691 AliFatal(Form("Unknown origin %d, Energy %f", origin, fEnergy));
692 }
693 }
694 else {
695 AliFatal(Form("Unknown energy %f", fEnergy));
696 }
697
698}
699
700