]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/background/AliAnalysisTaskBGvsTime.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Important update of the Physics selection: the BG cut based on cluster vs tracklets...
[u/mrichter/AliRoot.git] / PWG1 / background / AliAnalysisTaskBGvsTime.cxx
CommitLineData
dbec3946 1// AliAnalysisTaskBGvsTime
2
3// This task computes various histos as a function of time, using a
4// time stamp made from orbit, bx number and period. It is used for
5// studies of BG vs time and of luminosity.
6
7// The histograms are booked dinamically when they are called the
8// first time. This simplifies the handling of different trigger
9// classes (you don't need to know them in advance). The merging is
10// then complicated, but it is handled via the external class
11// AliHistoListWrapper.
12
13// Different sets of cuts can be used, and the binning of vs time
14// histos can be defined dynamically (see setters)
15
16// Can also process MC, in which case the efficiency x acceptance for
17// tracks used in the rate (luminosity) histos is also computed, and
18// put in a histo in which each bin corresponds to a given rate histo
19// (eta < 0.8 and pt > 0.5 or pt > 1.0 at the time I'm writing this).
20
d5f1fb98 21// The rates are corrected for the dead time, computed using trigger
22// scalers. WARNING: the dead time computation will NOT WORK on CAF
23// (it assumes consecutive events).
24
25// Author: Michele Floris, CERN
26
dbec3946 27#include "AliAnalysisTaskBGvsTime.h"
28#include "AliESDInputHandler.h"
29#include "TString.h"
30#include "AliVParticle.h"
31#include "AliESDInputHandlerRP.h"
32#include "TTree.h"
33#include "TClonesArray.h"
34#include "TBranch.h"
35//#include "AliITSRecPoint.h"
36#include "AliMultiplicity.h"
37#include "AliLog.h"
38#include "AliAnalysisManager.h"
dbec3946 39#include <iostream>
40#include "TFile.h"
41#include "TCanvas.h"
42#include "AliHistoListWrapper.h"
43#include "AliTriggerAnalysis.h"
44#include "TMath.h"
45#include "AliPhysicsSelection.h"
46#include "AliBackgroundSelection.h"
47#include "AliESDtrackCuts.h"
48#include "AliStack.h"
49#include "AliMCEvent.h"
50#include "TDatime.h"
51
52using namespace std;
53
54ClassImp(AliAnalysisTaskBGvsTime)
55
56
57
58AliAnalysisTaskBGvsTime::AliAnalysisTaskBGvsTime()
59 : AliAnalysisTaskSE("TaskBGvsTime"),
60 fESD(0),fListHisto(0),fListWrapper(0),fStartTime(0),fEndTime(0),
61 fNMultBins(0), fMultBins(0),fFirstTimeStamp(0), fLastTimeStamp(0), fBinWidth(10), fNoCuts(0), fPhysicsSelection(0), fUsePhysicsSelection(0),
62 fUseZeroBin(0), fIsMC(0), fSkipV0(0), fSkipZeroBin(0), fUseBunchInt(0), fHistoTimeStampVsUTC(0), fHistoTimeStampDiffUTC(0)
63
64{
65 // constructor
66
67 DefineOutput(1, AliHistoListWrapper::Class());
68 DefineOutput(2, AliPhysicsSelection::Class());
69 // DefineOutput(2, TH1I::Class());
70
71}
72AliAnalysisTaskBGvsTime::AliAnalysisTaskBGvsTime(const char * name)
73 : AliAnalysisTaskSE(name),fESD (0),fListHisto(0),fListWrapper(0),fStartTime(0),fEndTime(0),
74 fNMultBins(0),fMultBins(0),fFirstTimeStamp(0), fLastTimeStamp(0), fBinWidth(10), fNoCuts(0), fPhysicsSelection(0), fUsePhysicsSelection(0),
75 fUseZeroBin(0), fIsMC(0), fSkipV0(0), fSkipZeroBin(0), fUseBunchInt(0), fHistoTimeStampVsUTC(0), fHistoTimeStampDiffUTC(0)
76{
77 //
78 // Standard constructur which should be used
79 //
80
81 DefineOutput(1, AliHistoListWrapper::Class());
82 DefineOutput(2, AliPhysicsSelection::Class());
83 // DefineOutput(2, TH1I::Class());
84
85}
86
87AliAnalysisTaskBGvsTime::AliAnalysisTaskBGvsTime(const AliAnalysisTaskBGvsTime& obj) :
88 AliAnalysisTaskSE(obj) ,fESD (0),fListHisto(0),fListWrapper(0),fStartTime(0),fEndTime(0),
89 fNMultBins(0),fMultBins(0),fFirstTimeStamp(0), fLastTimeStamp(0), fBinWidth(10), fNoCuts(0), fPhysicsSelection(0), fUsePhysicsSelection(0),
90 fUseZeroBin(0), fIsMC(0), fSkipV0(0), fSkipZeroBin(0), fUseBunchInt(0), fHistoTimeStampVsUTC(0), fHistoTimeStampDiffUTC(0)
91{
92 //copy ctor
93 fESD = obj.fESD ;
94 fListHisto= obj.fListHisto;
95 fListWrapper= obj.fListWrapper;
96 fStartTime= obj.fStartTime;
97 fEndTime= obj.fEndTime;
98
99 fNMultBins= obj.fNMultBins;
100 fMultBins= obj.fMultBins;
101 fFirstTimeStamp= obj.fFirstTimeStamp;
102 fLastTimeStamp= obj.fLastTimeStamp;
103 fBinWidth= obj.fBinWidth;
104 fNoCuts= obj.fNoCuts;
105 fPhysicsSelection= obj.fPhysicsSelection;
106 fUsePhysicsSelection= obj.fUsePhysicsSelection;
107 fUseZeroBin= obj.fUseZeroBin;
108 fIsMC= obj.fIsMC;
109 fSkipV0= obj.fSkipV0;
110 fSkipZeroBin= obj.fSkipZeroBin;
111 fUseBunchInt= obj.fUseBunchInt;
112 fHistoTimeStampVsUTC = obj.fHistoTimeStampVsUTC;
113 fHistoTimeStampDiffUTC = obj.fHistoTimeStampDiffUTC;
114
115}
116
117AliAnalysisTaskBGvsTime::~AliAnalysisTaskBGvsTime(){
118
119 // destructor
120 if(fListWrapper) {
121 delete fListWrapper;
122 fListWrapper = 0;
123 }
124 if(fMultBins) {
125 delete fMultBins;
126 fMultBins = 0;
127 }
128 if(fPhysicsSelection) {
129 delete fPhysicsSelection;
130 fPhysicsSelection = 0;
131 }
132 // Histos should not be destroyed: fListWrapper is owner!
133}
134void AliAnalysisTaskBGvsTime::UserCreateOutputObjects()
135{
136 // Called once
137 fListWrapper = new AliHistoListWrapper("histoListWrapper","histoListWrapper");
138 fListHisto = fListWrapper->GetList();
139 Float_t lenght = fEndTime - fStartTime;
140 Int_t nBins = TMath::FloorNint(lenght/fBinWidth)+1;
141 Int_t end = nBins*fBinWidth;
142 // cover 20 days in UTC. I'll set offset by hand below.
143 Int_t lenghtUTC = 20*3600*24;
144 Int_t nBinsUTC = TMath::FloorNint(Double_t(lenghtUTC)/fBinWidth)+1;
145 Int_t endUTC = nBinsUTC*fBinWidth;
146 if (!fIsMC) {
147 Bool_t oldStatus = TH1::AddDirectoryStatus();
148 TH1::AddDirectory(kFALSE);
149 fHistoTimeStampVsUTC = new TH2F ("fHistoTimeStampVsUTC", "fHistoTimeStampVsUTC", nBins, -0.5, end-0.5, nBinsUTC, -0.5, endUTC-0.5);
150 fListHisto->Add(fHistoTimeStampVsUTC);
151 fHistoTimeStampDiffUTC = new TH1F("fHistoTimeStampDiffUTC", "fHistoTimeStampDiffUTC", nBinsUTC, -0.5, endUTC-0.5);
152 fListHisto->Add(fHistoTimeStampDiffUTC);
153 TH1::AddDirectory(oldStatus);
154 }
155
156 fPhysicsSelection = new AliPhysicsSelection();
d5f1fb98 157 fPhysicsSelection->SetUseBXNumbers();
dbec3946 158 if(fIsMC) fPhysicsSelection->SetAnalyzeMC();
159 fPhysicsSelection->SetBin0Callback(this->GetName());
160
161 fPhysicsSelection->SetSkipTriggerClassSelection();//
162
1ea7a921 163 // AliBackgroundSelection * bg = new AliBackgroundSelection();
164 // bg->SetDeltaPhiCut(10);
165 // fPhysicsSelection->AddBackgroundIdentification(bg);
dbec3946 166
167 if (fSkipV0) {
168 fPhysicsSelection->SetSkipV0();
169 }
170
171}
172
173
174void AliAnalysisTaskBGvsTime::UserExec(Option_t *)
175{
176
177 /* PostData(0) is taken care of by AliAnalysisTaskSE */
178 PostData(1,fListWrapper);
179 PostData(2,fPhysicsSelection);
180
d5f1fb98 181
dbec3946 182 const float etaCut = 0.8;
183
184
185 fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
2787875e 186 if(!fESD) {
187 AliFatal("Cannot get ESD");
188 }
dbec3946 189 if (strcmp(fESD->ClassName(),"AliESDEvent")) {
190 AliFatal("Not processing ESDs");
191 }
192
193 if (fUseBunchInt) {
194 fPhysicsSelection->SetComputeBG();
195 }
196
d5f1fb98 197 // Get V0 flags and trigger flags
198 static AliTriggerAnalysis * triggerAnalysis = new AliTriggerAnalysis();
199 Bool_t v0A = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0A);
200 Bool_t v0C = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0C);
201 Bool_t v0ABG = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0ABG);
202 Bool_t v0CBG = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0CBG);
203 Bool_t mb1Offline = triggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kMB1);
204 Bool_t isInV0 = !v0ABG && !v0CBG && ((v0A && !v0C) || (v0C && !v0A)); // try to select beam gas in CINT1A/C events: require one v0 interaction (but not the other) and not BG hit in v0. This should select events produced in between the 2 v0s and boosted forward.
205
206
dbec3946 207 // If it is MC: fill histo of generated events for efficiency calculations
208 if (fIsMC) {
209 Bool_t atLeastPt1MC = kFALSE;
210 Bool_t atLeastPt05MC = kFALSE;
211
212 if (!fMCEvent) {
213 AliError("No MC info found");
214 } else {
215
216 //loop on the MC event
217 for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
218 AliMCParticle *mcPart = (AliMCParticle*)fMCEvent->GetTrack(ipart);
219
220 // We don't care about neutrals and non-physical primaries
221 if(mcPart->Charge() == 0) continue;
222 if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue;
223
224 // Kinematic cuts:
225 if (TMath::Abs(mcPart->Eta()) < etaCut) {
d5f1fb98 226 if (mcPart->Pt() > 0.5) atLeastPt05MC = kTRUE;
227 if (mcPart->Pt() > 1.0) atLeastPt1MC = kTRUE;
dbec3946 228 }
229 if (atLeastPt1MC && atLeastPt05MC) break; // no need to look for other tracks
230 }
d5f1fb98 231 if (atLeastPt1MC ) {
232 GetEfficiencyHisto(kEffStepGen)->Fill(kEffPt1);
233 if(mb1Offline) GetEfficiencyHisto(kEffStepTrig)->Fill(kEffPt1);
234 }
235 if (atLeastPt05MC) {
236 GetEfficiencyHisto(kEffStepGen)->Fill(kEffPt05 );
237 if(mb1Offline) GetEfficiencyHisto(kEffStepTrig)->Fill(kEffPt05);
238 }
dbec3946 239 }
240 }
241
242
243
244
245 // CUTS
246 const AliMultiplicity* mult = fESD->GetMultiplicity();
247 if (!mult){
248 AliFatal("No multiplicity object"); // TODO: Should this be fatal?
249 }
250 // get number of SPD clusters
251 Int_t spdClusters = 0;
252 for(Int_t ilayer = 0; ilayer < 2; ilayer++){
253 spdClusters += mult->GetNumberOfITSClusters(ilayer);
254 }
255
256 // Bool_t isZeroBin = kFALSE;
257 Bool_t isZeroBin = IsEventInBinZero();
258 if(fUseZeroBin && !isZeroBin) return;
259 if(fSkipZeroBin && isZeroBin ) return;
260
261 Bool_t physelDecision = fPhysicsSelection->IsCollisionCandidate(fESD);
262
263 if (fUsePhysicsSelection) {
264 if(!physelDecision) return;
265 }
266 else if (!fNoCuts) {if (spdClusters < 2) return;}// At least 2 clusters
267 // else {AliInfo("No Cuts");}
268
269
270
271 // get time stamp
272 Long64_t timeStampBX = 0;
273 timeStampBX = fESD->GetBunchCrossNumber();
d5f1fb98 274 timeStampBX += (Long64_t) 3564 * (fESD->GetOrbitNumber() + fESD->GetPeriodNumber() * 16777216);
dbec3946 275 timeStampBX = (Long64_t) (25e-9 * timeStampBX);
276 if (fFirstTimeStamp == 0) {
277 fFirstTimeStamp = timeStampBX;
278 fLastTimeStamp = timeStampBX;
279 }
280 if (timeStampBX < fFirstTimeStamp) {
281 AliError("Time stamp not monothonic!");
282 fFirstTimeStamp = timeStampBX;
283 }
284 if (timeStampBX > fLastTimeStamp) {
285 fLastTimeStamp = timeStampBX;
286 }
287 timeStampBX -= fStartTime;
288
289 Long64_t timeStamp = timeStampBX;
290
291 Long64_t timeStampGDC = fESD->GetTimeStamp();
d5f1fb98 292
dbec3946 293 static TDatime timeOffsetDT(2009,12,5,0,0,0);
294 static Long64_t timeOffset = timeOffsetDT.Convert();
d5f1fb98 295
296
297
298 // Get trigger scalers for dead time calculation (only data)
299 // Only CINT1B (at least for the time being)
300 AliESDHeader* esdheader = (AliESDHeader*)fESD->GetHeader();
301 static ULong64_t L0 = 0;
302 static ULong64_t L2 = 0;
303 if (fESD->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")&&!fIsMC) {
304 AliTriggerScalersRecordESD* scalrecord = (AliTriggerScalersRecordESD*)esdheader->GetTriggerScalersRecord();
305 const AliTriggerScalersESD* scalers = scalrecord->GetTriggerScalersForClass(2); //2 is the cint1b class index in the trigger mask
306 L0 = scalers->GetLOCB(); //L0 before any vetos
307 L2 = scalers->GetL2CA(); //L2 after vetos
308 }
dbec3946 309
310
311 // loop over trigger classes in the event
312 TObjArray * tokens = 0;
313 if(fIsMC) {
314 // in case of montecarlo I override the trigger class to CINT1B for latter compatibility
315 tokens = new TObjArray;
316 tokens->SetOwner();
317 tokens->Add(new TObjString("CINT1B-ABCE-NOPF-ALL"));
318 }
319 else {
320 TString trgClasses = fESD->GetFiredTriggerClasses();
321 tokens = trgClasses.Tokenize(" ");
322 }
323 TIter iter(tokens);
324
325 while(TObjString * tok = (TObjString*) iter.Next()){
326 // clean up trigger name
327 TString trg = tok->GetString();
328 trg.Strip(TString::kTrailing, ' ');
329 trg.Strip(TString::kLeading, ' ');
330 // print selected events in !CIN1B trigs:
331 // if(!trg.Contains("CINT1B"))
332// Printf("File: %s, IEV: %d, TRG: %s, Orbit: 0x%x, Period: %d, BC: %d\n",
333// ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(),
334// trg.Data(),
335// fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
336
337
338
339 // Fill histos
340 GetVsTimeHistoAll(trg.Data())->Fill(timeStamp);
341 GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"ALL")->Fill(timeStamp);
342 GetVsTimeHisto(trg.Data(),-1, fESD->GetBunchCrossNumber(),"ALL")->Fill(timeStamp);
343 GetVsTimeHisto(trg.Data(),spdClusters,-1, "ALL")->Fill(timeStamp);
344 GetVsTimeHisto(trg.Data(),-1 ,-1, "ALL")->Fill(timeStamp);
345 if (isInV0) {
346 GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"inV0")->Fill(timeStamp);
347 GetVsTimeHisto(trg.Data(),-1 ,fESD->GetBunchCrossNumber(),"inV0")->Fill(timeStamp);
348 GetVsTimeHisto(trg.Data(),-1 ,-1, "inV0")->Fill(timeStamp);
349 }
350
351 // In order to compute mean multiplicity of spd clusters in a time
352 // bin, we integrate the multiplicity in that bin, and then we
353 // divide for the number of events in that bin in terminate.
354
355 // Is the error computed correctly if we fill with a weight ? Looping to be sure...
356 for (Int_t iclus = 0; iclus < spdClusters; iclus ++) {
357 GetVsTimeHisto((TString("hMultSPDvsTime_")+trg).Data())->Fill(timeStamp);
358 GetVsTimeHisto(trg.Data(),-1, fESD->GetBunchCrossNumber(),"ALL","MultSPDvsTime")->Fill(timeStamp);
359 GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"ALL","MultSPDvsTime")->Fill(timeStamp);
360 GetVsTimeHisto(trg.Data(),spdClusters,-1, "ALL","MultSPDvsTime")->Fill(timeStamp);
361 GetVsTimeHisto(trg.Data(),-1,-1, "ALL","MultSPDvsTime")->Fill(timeStamp);
362 if (isInV0) {
363 GetVsTimeHisto(trg.Data(),-1, fESD->GetBunchCrossNumber(),"inV0","MultSPDvsTime")->Fill(timeStamp);
364 GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"inV0","MultSPDvsTime")->Fill(timeStamp);
365 GetVsTimeHisto(trg.Data(),-1,-1, "inV0","MultSPDvsTime")->Fill(timeStamp);
366 }
367
368 }
369
370 // Keep cluster distribution for reference
371 GetDistributionHisto(trg.Data(),kDistSPDMult)->Fill(spdClusters);
372 if(isInV0) GetDistributionHisto(trg.Data(),kDistSPDMult,"_inV0")->Fill(spdClusters);
373
d5f1fb98 374 // Distribution of hits per its layer:
375 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
376 GetDistributionHisto(trg.Data(),kDistClsITSLayer)->Fill(ilayer,mult->GetNumberOfITSClusters(ilayer));// fill weighting with the number of CLS
377 }
378
379
380 // Dead time vs time stamp & time offset
381 // WARNING THIS WON'T WORK ON CAF, NOR GRID: REDO IT WITH A NEW MERGEABLE OBJECT
382 if (trg=="CINT1B-ABCE-NOPF-ALL"&&!fIsMC) {
383
384 // Fill time difference histos (only for CINt1B)
385 if (!fIsMC) {
386 fHistoTimeStampVsUTC->Fill(timeStampBX,timeStampGDC-timeOffset);
387 fHistoTimeStampDiffUTC->Fill(timeStampGDC-timeOffset-timeStampBX);
388 }
389
390 static ULong64_t oldL0 = 0; // L0 counts at the beginning of this bin
391 static ULong64_t oldL2 = 0; // L2 counts at the beginning of this bin
392 static ULong64_t prevL0 = 0; // L0 counts in the previous event
393 static ULong64_t prevL2 = 0; // L2 counts in the previous event
394
395 static ULong64_t oldTime = 0; // time stamp at the beginning of this bin
b56fd745 396 // static ULong64_t previousTime; // timestamp in the previous event
d5f1fb98 397
398 static Int_t prevbin = -1; // bin in the previous event
399
400 Int_t bin = GetDeadTimeHisto(trg.Data())->FindBin(timeStamp);
401
402 if (prevbin == -1) { // first event
403 prevbin = bin;
404 oldL0 = L0;
405 oldL2 = L2;
406 oldTime = timeStamp;
407
408 } else if (prevbin != bin) {
409 // New bin: let's fill the previous one
410 Double_t dL0 = Double_t(prevL0 - oldL0);
d5f1fb98 411 Double_t dL2 = Double_t(prevL2 - oldL2);
d5f1fb98 412
413 // Double_t deadtime = Double_t(1 - dL2/dL0); // interested in relative fraction of dead time
414 Double_t deadtime = Double_t(dL2/dL0); // interested in relative fraction of dead time
b56fd745 415 Double_t edeadtime = TMath::Sqrt(dL2*(dL0-dL2)/dL0/dL0/dL0); // Binomial error
d5f1fb98 416
b56fd745 417// cout << "DEADTIME " << endl;
418// cout << L0 << " " << dL0 << " " << oldL0 << " " << prevL0 << endl;
419// cout << L2 << " " << dL2 << " " << oldL2 << " " << prevL2 << endl;
420// cout << deadtime << endl;
d5f1fb98 421
422
423 GetDeadTimeHisto(trg.Data())->SetBinContent(prevbin, deadtime );
424 GetDeadTimeHisto(trg.Data())->SetBinError (prevbin, edeadtime);
425
426 oldL0 = L0;
427 oldL2 = L2;
428 oldTime = timeStamp;
429
430 }
431 prevbin = bin;
432 prevL0 = L0;
433 prevL2 = L2;
434
435 }
dbec3946 436
437 // TPC track multiplicity vs time
438
439 // Selection by andrea.dainese@pd.infn.it
440 // 15.03.2010
441
442
443 // Primary vertex
444 Bool_t badVertex = kFALSE;
445 const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
446 if(vertex->GetNContributors()<1) {
447 // SPD vertex
448 vertex = fESD->GetPrimaryVertexSPD();
449 if(vertex->GetNContributors()<1) {
450 badVertex = kTRUE;
451 }
452 }
453
454 // Fill vertex distribution
455 if ( ((vertex->IsFromVertexerZ() && vertex->GetDispersion()<=0.02) || !vertex->IsFromVertexerZ()) && !badVertex ) {
456 GetDistributionHisto(trg.Data(),kDistVertex)->Fill(vertex->GetZ());
457 if (vertex->IsFromVertexerZ()) GetDistributionHisto(trg.Data(),kDistVertexZ) ->Fill(vertex->GetZ());
458 else GetDistributionHisto(trg.Data(),kDistVertex3D)->Fill(vertex->GetZ());
459
460 if(isInV0) GetDistributionHisto(trg.Data(),kDistVertex,"_inV0")->Fill(vertex->GetZ());
461 }
462
463 // apply a cut |zVertex| < CUT, if needed
464
465 // Track cuts (except d0 cut)
466 //------- TPC track selection --------
467 // Selection by andrea.dainese@pd.infn.it
468 // 15.03.2010
b56fd745 469
470 Bool_t selectPrimaries=kTRUE;
471 static AliESDtrackCuts* esdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selectPrimaries);
dbec3946 472
473 // loop on tracks
474 Int_t ntracks = fESD->GetNumberOfTracks();
475 if(badVertex) ntracks = -1; // skip loop if the vertex is bad
476
477 // flags for luminosity histos
478 Bool_t atLeastPt1 = kFALSE;
479 Bool_t atLeastPt05 = kFALSE;
480
481 for (Int_t iTrack = 0; iTrack<ntracks; iTrack++) {
482 AliESDtrack * track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTrack));
483 // for each track
484
485 // track quality cuts
486 if(!esdtrackCutsITSTPC->AcceptTrack(track)) continue;
487
488 // bring it to the primary vertex and compute impact parameters
489 if(!track->RelateToVertex(vertex,fESD->GetMagneticField(),kVeryBig)) continue; // this is already done in AliReconstruction...
490
491 // track-to-vertex cut (see below)
492 if(!SelectOnImpPar(track)) continue;
493
494 // Fill histos (TPC Multiplicity vs TIME)
495 GetVsTimeHisto(trg.Data(),-1, fESD->GetBunchCrossNumber(),"ALL","MultTPCvsTime")->Fill(timeStamp);
496 GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"ALL","MultTPCvsTime")->Fill(timeStamp);
497 GetVsTimeHisto(trg.Data(),spdClusters,-1, "ALL","MultTPCvsTime")->Fill(timeStamp);
498 GetVsTimeHisto(trg.Data(),-1,-1, "ALL","MultTPCvsTime")->Fill(timeStamp);
499 if ((v0A || v0C) && !(v0ABG || v0CBG)) {
500 GetVsTimeHisto(trg.Data(),-1, fESD->GetBunchCrossNumber(),"inV0","MultTPCvsTime")->Fill(timeStamp);
501 GetVsTimeHisto(trg.Data(),spdClusters,fESD->GetBunchCrossNumber(),"inV0","MultTPCvsTime")->Fill(timeStamp);
502 GetVsTimeHisto(trg.Data(),-1,-1, "inV0","MultTPCvsTime")->Fill(timeStamp);
503 }
504
dbec3946 505
506 // has the event at least one track satisfying the required conditions? (used for rate/luminosity)
507 if (TMath::Abs(track->Eta()) < etaCut) {
d5f1fb98 508 // Fill histo (pt distribution with standard cuts)
509 GetDistributionHisto(trg.Data(),kDistPt)->Fill(track->Pt());
dbec3946 510 if (track->Pt() > 0.5) {
511 atLeastPt05 = kTRUE;
512 }
513 if (track->Pt() > 1.0) {
514 atLeastPt1 = kTRUE;
515 }
516 }
517 }
518
b56fd745 519 // TEMPORARY LOOP: FILL DISTRIBUTION OF PT FOR 2 CLASSES OF EVENTS
520 for (Int_t iTrack = 0; iTrack<ntracks; iTrack++) {
521 AliESDtrack * track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTrack));
522 // for each track
523
524 // track quality cuts
525 if(!esdtrackCutsITSTPC->AcceptTrack(track)) continue;
526
527 // bring it to the primary vertex and compute impact parameters
528 if(!track->RelateToVertex(vertex,fESD->GetMagneticField(),kVeryBig)) continue; // this is already done in AliReconstruction...
529
530 // track-to-vertex cut (see below)
531 if(!SelectOnImpPar(track)) continue;
532 if (TMath::Abs(track->Eta()) < etaCut){
533 if (atLeastPt05) GetDistributionHisto(trg.Data(),kDistPt, "_atLeast05")->Fill(track->Pt());
534 if (atLeastPt1) GetDistributionHisto(trg.Data(),kDistPt, "_atLeast1" )->Fill(track->Pt());
535 }
536 }
537 // END OF TEMPORARY LOOP
538
dbec3946 539 // Fill histos for luminosity: rate of events with at least one
540 // track in the pseudo rapidity region |eta| < 0.8 and pt > 0.5 or
541 // 1 GeV
542 if (atLeastPt05) {
543 GetVsTimeHisto(GetVsTimeHistoForLuminosityName(trg.Data(), 0.5))->Fill(timeStamp);
d5f1fb98 544 if(fIsMC) GetEfficiencyHisto(kEffStepRec)->Fill(kEffPt05);
dbec3946 545 }
546 if (atLeastPt1) {
547
548 if(GetVsTimeHisto(GetVsTimeHistoForLuminosityName(trg.Data(), 1.0))->Fill(timeStamp) < 0) {
549 AliWarning(Form("Timestamp out of range %lld", timeStamp));
550 };
d5f1fb98 551 if(fIsMC) GetEfficiencyHisto(kEffStepRec)->Fill(kEffPt1);
552 }
553
554
555 // TPC TRACKS: loose selection in order to keep some BG track
556 ntracks = fESD->GetNumberOfTracks(); // Ignore vertex selecion
557 for(Int_t itrack = 0; itrack < ntracks; itrack++){
558 AliESDtrack * track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(itrack));
559 // for each track
560 // Same as AliESDtrackCuts::GetStandardTPCOnlyTrackCuts() , but not using DCA cut
561
562 AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
563
564 esdTrackCuts->SetMinNClustersTPC(50);
565 esdTrackCuts->SetMaxChi2PerClusterTPC(4);
566 esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
567
568// esdTrackCuts->SetMaxDCAToVertexZ(3.2);
569// esdTrackCuts->SetMaxDCAToVertexXY(2.4);
570// esdTrackCuts->SetDCAToVertex2D(kTRUE);
571
572 if (!esdTrackCuts->AcceptTrack(track)) continue;
573 // Fill pt and DCA distribution
574 // pt
575 GetDistributionHisto(trg.Data(),kDistPtLoose)->Fill(track->Pt());
576 if(isInV0) GetDistributionHisto(trg.Data(),kDistPtLoose, "_inV0")->Fill(track->Pt());
577 // dca (on the xy plane)
578 Float_t xy_dca, z_dca;
579 track->GetImpactParameters(xy_dca,z_dca);
580 GetDistributionHisto(trg.Data(),kDistDCATPC)->Fill(xy_dca);
581 if(isInV0) GetDistributionHisto(trg.Data(),kDistDCATPC, "_inV0")->Fill(track->Pt());
582
dbec3946 583 }
584
d5f1fb98 585
dbec3946 586 }
587
588}
589
590void AliAnalysisTaskBGvsTime::Terminate(Option_t *){
591
592 // normalize and rescale histos
593 AliInfo("Normalizing multiplicity histo") ;
594
595 fListWrapper = dynamic_cast<AliHistoListWrapper*> (GetOutputData(1));
596 if (!fListWrapper){
597 AliError("Cannot get list wrapper");
598 }
599 fListHisto = fListWrapper->GetList();
600 fHistoTimeStampDiffUTC = (TH1F*) fListHisto->FindObject("fHistoTimeStampDiffUTC");
601 fHistoTimeStampVsUTC = (TH2F*) fListHisto->FindObject("fHistoTimeStampVsUTC");
602
603 // Divide rate histos for the number of events
604 TIterator * iter = fListHisto->MakeIterator();
605 TH1 * h = 0;
606 while ((h = (TH1*) iter->Next())) {
607 if((TString(h->GetName())).Contains("hMultSPDvsTime_") ||
608 (TString(h->GetName())).Contains("hMultTPCvsTime_")
609 ){
610
611 AliInfo(Form("Normalizing %s",h->GetName()));
612 // continue;
613 TString histoname = h->GetName();
614 histoname.ReplaceAll("MultSPDvsTime","");
615 histoname.ReplaceAll("MultTPCvsTime","");
616
617 TH1* hev = (TH1*)fListHisto->FindObject(histoname.Data());
618 if(!hev) {
619 AliError(Form(" -> Cannot find events histo %s",histoname.Data()));
620 continue;
621 }
622 AliInfo (Form(" with histo %s",hev->GetName()));
623 // Errors on ev num should be ignored in the division
624 Int_t nbin = h->GetNbinsX();
625 for(Int_t ibin = 1; ibin <=nbin; ibin++) {
626 if(hev->GetBinContent(ibin) != 0) {
627 h->SetBinContent(ibin, h->GetBinContent(ibin)/ hev->GetBinContent(ibin));
628 h->SetBinError (ibin, h->GetBinError(ibin) / hev->GetBinContent(ibin));
629 }
630 }
631 }
632 if((TString(h->GetName())).Contains("hRate")){
633 h->Scale(1.,"width"); // divide for bin width to obtain a rate
634 }
635 }
636
637 // Compute Efficiency:
638 if(fIsMC) {
639 TH1F* hgen =(TH1F*) fListHisto->FindObject("hEffGen");
d5f1fb98 640 TH1F* htrg =(TH1F*) fListHisto->FindObject("hEffTrig");
dbec3946 641 TH1F* hrec =(TH1F*) fListHisto->FindObject("hEffRec");
d5f1fb98 642 if (!hgen || !hrec || !htrg) {
dbec3946 643 AliError("Cannot find eff histos");
644 }
645 else {
646 TH1F* heff = (TH1F*) fListHisto->FindObject("hEffRatio");
d5f1fb98 647 TH1F* hefftrg = (TH1F*) fListHisto->FindObject("hEffRatioTrg");
dbec3946 648 if (heff) {
649 AliWarning("hEffRatio already in output list?");
650 }
651 else {
652 heff = (TH1F*) hgen->Clone("hEffRatio");
d5f1fb98 653 }
dbec3946 654 heff->Reset();
655 heff->Divide(hrec,hgen,1,1,"B");
656 fListHisto->Add(heff);
d5f1fb98 657
658 if (hefftrg) {
659 AliWarning("hEffRatioTrg already in output list?");
660 }
661 else {
662 hefftrg = (TH1F*) hgen->Clone("hEffRatioTrg");
663 }
664 hefftrg->Reset();
665 hefftrg->Divide(htrg,hgen,1,1,"B");
666 fListHisto->Add(hefftrg);
667
dbec3946 668 }
669 }
670
671
672 AliInfo(Form("Time interval: %lld -- %lld",fFirstTimeStamp,fLastTimeStamp));
673
674 AliInfo("Saving physics selection histos");
675 fPhysicsSelection = dynamic_cast<AliPhysicsSelection*> (GetOutputData(2));
676
677 TFile* fout = new TFile("event_stat.root", "RECREATE");
678
679 if (fPhysicsSelection)
680 {
681 fPhysicsSelection->Print();
682 fPhysicsSelection->SaveHistograms("physics_selection");
683 }
684
685 fout->Write();
686 fout->Close();
687
688}
689
690TH1F * AliAnalysisTaskBGvsTime::BookVsTimeHisto(const char * name, const char * title){
691
692 // Book istograms vs time
693
694 AliInfo(Form("Booking histo %s",name));
695
696 Bool_t oldStatus = TH1::AddDirectoryStatus();
697 TH1::AddDirectory(kFALSE);
698
699 // static Int_t nBins = fEndTime - fStartTime + 1;
700 // Compute bin width and range.
701 // if time interval is not a multiple of bin width drop the last bin
702
703 Float_t lenght = fEndTime - fStartTime;
704 Int_t nBins = TMath::FloorNint(lenght/fBinWidth)+1;
705 Float_t * bins = new Float_t[nBins+1];
706 for(Int_t ibin = 0; ibin <= nBins; ibin++){
707 Float_t edge = ibin*fBinWidth - 0.5;
708 bins[ibin] = edge < (lenght - 0.5) ? edge : (lenght-0.5);
709 }
710
711 // Int_t end = nBins*fBinWidth;
712
713 // TH1F * h = new TH1F(name,title, nBins, - 0.5, fEndTime - fStartTime + 0.5);
714 TH1F * h = new TH1F(name,title, nBins, bins);
715 fListHisto->Add(h);
716 h->SetXTitle("orbit");
717 h->Sumw2();
718
719 TH1::AddDirectory(oldStatus);
720
d5f1fb98 721
63424920 722 delete [] bins;
dbec3946 723 return h;
724}
725
726TH1F * AliAnalysisTaskBGvsTime::GetVsTimeHisto(const char * name) {
727
728 // Returns vs time histo. If not existing, creates it.
729
730 TH1F * h = (TH1F*) fListHisto->FindObject(name);
731 if(!h) h = BookVsTimeHisto(name,name);
732 return h;
733
734}
735
d5f1fb98 736TH1F * AliAnalysisTaskBGvsTime::GetEfficiencyHisto(Int_t step) {
dbec3946 737
738 // Return efficiency histo. If not existing, creates it.
739 // 1 bin per category
740
741 // Probability of reconstructing one event with at least one
742 // particle in a given category
743 // This is only created for MC.
744
745 // If the first argument is true, returns the histo at generation
746 // level, otherwise returns the histo at rec level. Fill this histo
747 // with elements of the enum kEffPt1, kEffPt05, ...
748
d5f1fb98 749 const char * name = 0;
750 switch(step) {
751 case kEffStepGen:
752 name = "hEffGen";
753 break;
754 case kEffStepRec:
755 name = "hEffRec";
756 break;
757 case kEffStepTrig:
758 name = "hEffTrig";
759 break;
760 }
761
dbec3946 762
763 TH1F * h = (TH1F*) fListHisto->FindObject(name);
764 if(!h) {
765 Bool_t oldStatus = TH1::AddDirectoryStatus();
766 TH1::AddDirectory(kFALSE);
767 h = new TH1F (name,name, kNEff, -0.5, kNEff-0.5);
768 h->GetXaxis()->SetBinLabel(h->FindBin(kEffPt1) , "pt > 1.0 GeV");
769 h->GetXaxis()->SetBinLabel(h->FindBin(kEffPt05), "pt > 0.5 GeV");
770 h->Sumw2();
771 fListHisto->Add(h);
772 TH1::AddDirectory(oldStatus);
773 }
774
775 return h;
776
777}
778
d5f1fb98 779const char * AliAnalysisTaskBGvsTime::GetVsTimeHistoForLuminosityName(const char * triggerClass, Float_t ptmin) {
dbec3946 780
781 // Compose name of rate histos
782
783 static TString name;
784 name = "hRate_";
d5f1fb98 785 name += triggerClass;
dbec3946 786
787 name += Form ("_etamin_0.8_ptmin_%2.2f",ptmin);
788
789 return name.Data();
790}
791
792
d5f1fb98 793const char * AliAnalysisTaskBGvsTime::GetVsTimeHistoName(const char * triggerClass, Int_t nclusters, Int_t bx, const char * v0flag, const char *prefix){
dbec3946 794
795 // compose name of vs time histos w/ different cuts
796
797 static TString name;
798 name = "h";
799 name = name +prefix+"_";
d5f1fb98 800 name += triggerClass;
dbec3946 801
802 if (nclusters >= 0) {
803 if(fMultBins){
804 // Add multiplicity label
805 Int_t selected_bin = -1;
806 for(Int_t ibin = 0; ibin < (fNMultBins-1); ibin++){
807 if (nclusters < fMultBins[ibin+1]) {
808 selected_bin = ibin;
809 break;
810 }
811 }
812 if(selected_bin < 0) {
813 AliError(Form("Cannot get mult bin - %d", nclusters));
814 selected_bin = (fNMultBins-1);
815 }
816 name += Form("_SPDCls_%d-%d",fMultBins[selected_bin],fMultBins[selected_bin+1]);
817 // AliInfo(Form("Mult bin: %d, %d, %s", nclusters, selected_bin, name.Data()));
818 }
819 } else {
820 name += "_SPDCls_ALL";
821 }
822 if(bx < 0) name = name +"_BXALL_"+v0flag;
823 else name = name +"_BX"+long(bx)+"_"+v0flag;
824 return name.Data();
825}
d5f1fb98 826const char * AliAnalysisTaskBGvsTime::GetVsTimeHistoNameAll(const char * triggerClass) {
dbec3946 827
828 // Compose default name of vstime histos in a given trigger class
829
830 static TString name;
831 name = "h_";
d5f1fb98 832 name += triggerClass;
dbec3946 833 return name.Data();
834
835}
836
d5f1fb98 837TH1F * AliAnalysisTaskBGvsTime::GetDistributionHisto(const char * triggerClass, Int_t dist, const char * suffix) {
dbec3946 838
839 // Returns distributions histos. If not existing, creates it.
840
841 // Possible distributions:
842 // - SPD cluster
843 // - TPC tracks multiplicity
844 // - TPC tracks pt
845 // - Vz
d5f1fb98 846 // - TPC tracks DCA
847 // - cluster per ITS layer
dbec3946 848
849 TString name ;
850 TString title ;
851 TString xtitle;
852 Int_t nbin =0;
ca3f2e06 853 Float_t min =0;
854 Float_t max =0;
dbec3946 855
856
857 if (dist == kDistSPDMult) {
858 name = "hMultSPD_";
859 title = "SPD Cluster Multiplicity (";
860 xtitle = "SPD Clusters";
861 nbin = 50;
862 min = 0;
863 max = 100;
864 } else if (dist == kDistTPCMult) {
865 name = "hMultTPC_";
866 title = "TPC Tracks Multiplicity (";
867 xtitle = "TPC Tracks";
868 nbin = 25;
869 min = 0;
870 max = 50;
d5f1fb98 871 } else if (dist == kDistPtLoose) {
872 name = "hPtLoose_";
873 title = "p_{T} distribution - TPC, loose cuts (";
874 xtitle = "p_{T} (GeV)";
875 nbin = 50;
876 min = 0;
877 max = 10;
dbec3946 878 } else if (dist == kDistPt) {
879 name = "hPt_";
d5f1fb98 880 title = "p_{T} distribution - TPC, standard cuts, |#eta|<0.8 (";
dbec3946 881 xtitle = "p_{T} (GeV)";
d5f1fb98 882 nbin = 100;
dbec3946 883 min = 0;
884 max = 10;
885 } else if (dist == kDistVertex || dist == kDistVertexZ || dist == kDistVertex3D ) {
886 nbin = 120;
887 min = -30;
888 max = 30;
889 name = "hVz_";
890 title = "V_{z} distribution";
891 if (dist == kDistVertexZ) {
892 title += " - Vertexer Z (";
893 name+= "Z_";
894 }
895 else if (dist == kDistVertex3D) {
896 title += " - Vertexer 3D (";
897 name+= "3D_";
898 }
899 else title += " (";
900 xtitle = "V_{z} (cm)";
901 } else if (dist == kDistDCATPC) {
902 nbin = 50;
903 min = -0.5;
904 max = +0.5;
905 name = "hDCA_";
d5f1fb98 906 title = "TPC DCA distribution - loose cuts (";
907 xtitle = "DCA";
908 } else if (dist == kDistClsITSLayer) {
909 nbin = 6;
910 min = -0.5;
911 max = 5.5;
912 name = "hClsITS_";
913 title = "Cluster per ITS layer (";
914 xtitle = "ITS Layer";
915 }else {
dbec3946 916 AliError(Form("Distribution type not supported: %d",dist));
917 return 0;
918 }
919
920
d5f1fb98 921 name += triggerClass;
dbec3946 922 if (suffix) name += suffix;
923 title = title+name+")";
924
925 TH1F* h = (TH1F*) fListHisto->FindObject(name.Data());
926
927 if(!h) h = BookDistributionHisto(name.Data(), title.Data(), xtitle.Data(), nbin, min, max);
928 return h;
929
930}
931
ca3f2e06 932TH1F * AliAnalysisTaskBGvsTime::BookDistributionHisto(const char * name, const char * title, const char * xtitle, Int_t nbin, Float_t min, Float_t max) {
dbec3946 933
934 // Book distributions histos
935
936 AliInfo(Form("Booking histo %s",name));
937
938 Bool_t oldStatus = TH1::AddDirectoryStatus();
939 TH1::AddDirectory(kFALSE);
940
941 TH1F * h = new TH1F (name,title,nbin,min,max);
942 h->Sumw2();
943 h->SetXTitle(xtitle);
944 // h->SetYTitle("N");
945 fListHisto->Add(h);
946 TH1::AddDirectory(oldStatus);
947
948 return h;
949
950}
951
d5f1fb98 952TH1F * AliAnalysisTaskBGvsTime::GetDeadTimeHisto(const char * triggerClass) {
953 // returns histo of dead time vs timestamp for a given trigger
954 // class. If the histo does not exist, it books it.
955 TString name = "hDeadTime_" ;
956 TString title = "Dead Time vs TimeStamp (";
957
958 name += triggerClass;
959 title = title+triggerClass+")";
960
961
962 TH1F* h = (TH1F*) fListHisto->FindObject(name.Data());
963
964 if(!h) {
965 h = BookVsTimeHisto(name.Data(), title.Data());
966 h->SetYTitle("deadtime");
967 }
968 return h;
969
970
971}
972
973
974// TH2F * AliAnalysisTaskBGvsTime::BookDeadTimeHisto(const char * name, const char * title) {
975// // Book dead time vs time stamp histos
976
977// AliInfo(Form("Booking histo %s",name));
978
979// Bool_t oldStatus = TH1::AddDirectoryStatus();
980// TH1::AddDirectory(kFALSE);
981
982// Float_t lenght = fEndTime - fStartTime;
983// Int_t nBins = TMath::FloorNint(lenght/fBinWidth)+1;
984// Double_t * bins = new Double_t[nBins+1];
985// for(Int_t ibin = 0; ibin <= nBins; ibin++){
986// Float_t edge = ibin*fBinWidth - 0.5;
987// bins[ibin] = edge < (lenght - 0.5) ? edge : (lenght-0.5);
988// }
989
990// TH2F * h = new TH2F (name,title,nBins,bins,200,0.,2.);
991// h->Sumw2();
992// h->SetXTitle("Time (s)");
993// h->SetYTitle("Dead Time");
994// fListHisto->Add(h);
995// TH1::AddDirectory(oldStatus);
996
997// delete bins;
998
999// return h;
1000
1001// }
1002
1003
dbec3946 1004
1005Bool_t AliAnalysisTaskBGvsTime::IsEventInBinZero() {
1006
1007 // Returns true if an event is to be assigned to the zero bin
1008
1009 Bool_t isZeroBin = kTRUE;
1010 const AliESDEvent* esd=fESD;
1011 const AliMultiplicity* mult = esd->GetMultiplicity();
1012 if (!mult){
1013 Printf("AliAnalysisTaskBGvsTime::IsBinZero: Can't get mult object");
1014 return kFALSE;
1015 }
1016 Int_t ntracklet = mult->GetNumberOfTracklets();
1017 const AliESDVertex * vtxESD = esd->GetPrimaryVertexSPD();
1018 if(vtxESD) {
1019 // If there is a vertex from vertexer z with delta phi > 0.02 we
1020 // don't consider it rec (we keep the event in bin0). If quality
1021 // is good eneough we check the number of tracklets
1022 // if the vertex is more than 15 cm away, this is autamatically bin0
1023 if( TMath::Abs(vtxESD->GetZ()) <= 15 ) {
1024 if (vtxESD->IsFromVertexerZ()) {
1025 if (vtxESD->GetDispersion()<=0.02 ) {
1026 if(ntracklet>0) isZeroBin = kFALSE;
1027 }
1028 } else if(ntracklet>0) isZeroBin = kFALSE; // if the event is not from Vz we chek the n of tracklets
1029 }
1030 }
1031 return isZeroBin;
1032
1033}
1034
1035Bool_t AliAnalysisTaskBGvsTime::SelectOnImpPar(AliESDtrack* t) {
1036 // from andrea dainese
1037 // cut on transverse impact parameter
1038 Float_t d0z0[2],covd0z0[3];
1039 t->GetImpactParameters(d0z0,covd0z0);
1040 Float_t sigma= 0.0050+0.0060/TMath::Power(t->Pt(),0.9);
1041 Float_t d0max = 7.*sigma;
1042 if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
1043 return kFALSE;
1044}
1045
1046
1047
1048