1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //_________________________________________________________________________
17 // Count events with different selections
19 // It produces a histogram with the number of events with 9 bins:
20 // 0: all events (that passed the physics selection if it was on)
21 // 1: same but cross check that event pointer did exist
22 // 2: passes vertex cut
23 // 3: passes track number cut, tracks for eta < 0.8
29 // 9: not pileup from SPD
34 // 14: 10 && 2 && 3 && 5
38 // Author: Gustavo Conesa Balbastre (LPSC)
40 //_________________________________________________________________________
46 #include <TProfile2D.h>
48 #include <TClonesArray.h>
49 #include <TGeoGlobalMagField.h>
50 #include "AliAODHeader.h"
51 //#include "AliTriggerAnalysis.h"
52 #include "AliESDEvent.h"
53 #include "AliAODEvent.h"
54 #include "AliESDtrackCuts.h"
55 #include "AliAnalysisManager.h"
56 #include "AliInputEventHandler.h"
58 #include "AliAnalysisTaskCounter.h"
59 ClassImp(AliAnalysisTaskCounter)
61 //______________________________________________________________
62 AliAnalysisTaskCounter::AliAnalysisTaskCounter(const char *name)
63 : AliAnalysisTaskSE(name),
64 fAcceptFastCluster(kTRUE),
66 fTrackMultEtaCut(0.8),
68 fOutputContainer(0x0),
69 fESDtrackCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
70 //fTriggerAnalysis (new AliTriggerAnalysis),
71 fCurrFileName(0), fCheckMCCrossSection(kFALSE),
73 fhXVertex(0), fhYVertex(0), fhZVertex(0),
74 fhXGoodVertex(0),fhYGoodVertex(0),fhZGoodVertex(0),
75 fhCentrality(0), fhEventPlaneAngle(0),
76 fh1Xsec(0), fh1Trials(0)
79 DefineOutput(1, TList::Class());
82 //______________________________________________
83 AliAnalysisTaskCounter::AliAnalysisTaskCounter()
84 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskCounter"),
85 fAcceptFastCluster(kTRUE),
87 fTrackMultEtaCut(0.8),
89 fOutputContainer(0x0),
90 fESDtrackCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
91 //fTriggerAnalysis (new AliTriggerAnalysis),
92 fCurrFileName(0), fCheckMCCrossSection(kFALSE),
94 fhXVertex(0), fhYVertex(0), fhZVertex(0),
95 fhXGoodVertex(0),fhYGoodVertex(0),fhZGoodVertex(0),
96 fhCentrality(0), fhEventPlaneAngle(0),
97 fh1Xsec(0), fh1Trials(0)
100 DefineOutput(1, TList::Class());
103 //__________________________________________________
104 AliAnalysisTaskCounter::~AliAnalysisTaskCounter()
108 if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
112 fOutputContainer->Delete() ;
113 delete fOutputContainer ;
116 if(fESDtrackCuts) delete fESDtrackCuts;
117 //if(fTriggerAnalysis) delete fTriggerAnalysis;
122 //____________________________________________________
123 void AliAnalysisTaskCounter::UserCreateOutputObjects()
127 fOutputContainer = new TList();
129 if(fCheckMCCrossSection)
131 fh1Xsec = new TH1F("hXsec","xsec from pyxsec.root",1,0,1);
132 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
133 fOutputContainer->Add(fh1Xsec);
135 fh1Trials = new TH1F("hTrials","trials root file",1,0,1);
136 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
137 fOutputContainer->Add(fh1Trials);
140 fhZVertex = new TH1F("hZVertex", " Z vertex distribution" , 200 , -50 , 50 ) ;
141 fhZVertex->SetXTitle("v_{z} (cm)");
142 fOutputContainer->Add(fhZVertex);
144 fhZGoodVertex = new TH1F("hZGoodVertex", " Good Z vertex distribution" , 200 , -50 , 50 ) ;
145 fhZGoodVertex->SetXTitle("v_{z} (cm)");
146 fOutputContainer->Add(fhZGoodVertex);
148 fhXVertex = new TH1F("hXVertex", " X vertex distribution" , 200 , -2 , 2 ) ;
149 fhXVertex->SetXTitle("v_{x} (cm)");
150 fOutputContainer->Add(fhXVertex);
152 fhXGoodVertex = new TH1F("hXGoodVertex", " Good X vertex distribution" , 200 , -2 , 2 ) ;
153 fhXGoodVertex->SetXTitle("v_{x} (cm)");
154 fOutputContainer->Add(fhXGoodVertex);
156 fhYVertex = new TH1F("hYVertex", " Y vertex distribution" , 200 , -2 , 2 ) ;
157 fhYVertex->SetXTitle("v_{y} (cm)");
158 fOutputContainer->Add(fhYVertex);
160 fhYGoodVertex = new TH1F("hYGoodVertex", " Good Y vertex distribution" , 200 , -2 , 2 ) ;
161 fhYGoodVertex->SetXTitle("v_{y} (cm)");
162 fOutputContainer->Add(fhYGoodVertex);
164 fhCentrality = new TH1F("hCentrality","Number of events in centrality bin, |vz|<10 cm, method <V0M> ",100,0.,100.) ;
165 fhCentrality->SetXTitle("Centrality bin");
166 fOutputContainer->Add(fhCentrality) ;
168 fhEventPlaneAngle=new TH1F("hEventPlaneAngle","Number of events in event plane, |vz|<10 cm, method <V0> ",100,0.,TMath::Pi()) ;
169 fhEventPlaneAngle->SetXTitle("EP angle (rad)");
170 fOutputContainer->Add(fhEventPlaneAngle) ;
172 fhNEvents = new TH1I("hNEvents", "Number of analyzed events", 21, 0, 21) ;
173 //fhNEvents->SetXTitle("Selection");
174 fhNEvents->SetYTitle("# events");
175 fhNEvents->GetXaxis()->SetBinLabel(1 ,"1 = PS");
176 fhNEvents->GetXaxis()->SetBinLabel(2 ,"2 = 1 & ESD");
177 fhNEvents->GetXaxis()->SetBinLabel(3 ,"3 = 2 & |Z|<10");
178 fhNEvents->GetXaxis()->SetBinLabel(4 ,"4 = 2 & |eta|<0.8");
179 fhNEvents->GetXaxis()->SetBinLabel(5 ,"5 = 3 & 4");
180 fhNEvents->GetXaxis()->SetBinLabel(6 ,"6 = 2 & V0AND");
181 fhNEvents->GetXaxis()->SetBinLabel(7 ,"7 = 3 & 6");
182 fhNEvents->GetXaxis()->SetBinLabel(8 ,"8 = 4 & 6");
183 fhNEvents->GetXaxis()->SetBinLabel(9 ,"9 = 5 & 6");
184 fhNEvents->GetXaxis()->SetBinLabel(10,"10 = 2 & not pileup");
185 fhNEvents->GetXaxis()->SetBinLabel(11,"11 = 2 & good vertex");
186 fhNEvents->GetXaxis()->SetBinLabel(12,"12 = 3 & 11");
187 fhNEvents->GetXaxis()->SetBinLabel(13,"13 = 4 & 11");
188 fhNEvents->GetXaxis()->SetBinLabel(14,"14 = 6 & 11");
189 fhNEvents->GetXaxis()->SetBinLabel(15,"15 = 9 & 11");
190 fhNEvents->GetXaxis()->SetBinLabel(16,"16 = 10 & 11");
191 fhNEvents->GetXaxis()->SetBinLabel(17,"17 = 6 & 10");
192 fhNEvents->GetXaxis()->SetBinLabel(17,"17 = 1 & |Z|<50");
193 fhNEvents->GetXaxis()->SetBinLabel(18,"18 = Reject EMCAL 1");
194 fhNEvents->GetXaxis()->SetBinLabel(19,"19 = 18 & 2");
195 fhNEvents->GetXaxis()->SetBinLabel(20,"20 = Reject EMCAL 2");
196 fhNEvents->GetXaxis()->SetBinLabel(21,"21 = 20 & 2");
198 fOutputContainer->Add(fhNEvents);
200 fOutputContainer->SetOwner(kTRUE);
202 PostData(1,fOutputContainer);
207 //_______________________________________________
208 void AliAnalysisTaskCounter::UserExec(Option_t *)
211 // Called for each event
213 //printf("___ Event __ %d __\n",(Int_t)Entry());
217 fhNEvents->Fill(0.5);
219 AliVEvent * event = InputEvent();
220 AliESDEvent * esdevent = dynamic_cast<AliESDEvent*> (event);
221 AliAODEvent * aodevent = dynamic_cast<AliAODEvent*> (event);
223 // Init mag field for tracks in case of ESDs, needed, not clear why
224 if (!TGeoGlobalMagField::Instance()->GetField() && esdevent) esdevent->InitMagneticField();
226 TString triggerclasses = event->GetFiredTriggerClasses();
228 //printf("Trigger class fired: %s \n",event->GetFiredTriggerClasses().Data());
230 if (triggerclasses.Contains("FAST") && !triggerclasses.Contains("ALL") && !fAcceptFastCluster)
232 //printf("Do not count events from fast cluster, trigger name %s\n",triggerclasses.Data());
236 fhNEvents->Fill(1.5);
239 Bool_t bSelectVZ = kFALSE;
240 Bool_t bV0AND = kFALSE;
241 Bool_t bPileup = kFALSE;
242 Bool_t bGoodV = kFALSE;
243 Bool_t bSelectTrack = kFALSE;
246 //---------------------------------
247 //Get the primary vertex, cut on Z
248 //---------------------------------
250 event->GetPrimaryVertex()->GetXYZ(v) ;
251 fhXVertex->Fill(v[0]);
252 fhYVertex->Fill(v[1]);
253 fhZVertex->Fill(v[2]);
255 if(TMath::Abs(v[2]) < fZVertexCut)
258 fhNEvents->Fill(2.5);
260 //else printf("Vertex out %f \n",v[2]);
262 //--------------------------------------------------
263 //Count tracks, cut on number of tracks in eta < 0.8
264 //--------------------------------------------------
265 Int_t nTracks = event->GetNumberOfTracks() ;
266 for (Int_t itrack = 0; itrack < nTracks; itrack++)
267 {////////////// track loop
268 AliVTrack * track = (AliVTrack*)event->GetTrack(itrack) ; // retrieve track from esd
271 if(esdevent && !fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
274 if(aodevent && !((AliAODTrack*)track)->IsHybridGlobalConstrainedGlobal()) continue ;
276 //Do not count tracks out of acceptance cut
277 if(TMath::Abs(track->Eta())< fTrackMultEtaCut) trackMult++;
280 //printf("AliAnalysisTaskCounter::UserExec() - Track Mult %d \n",trackMult);
282 //--------------------------------------------------
283 // At least one track
284 //--------------------------------------------------
287 bSelectTrack = kTRUE;
288 fhNEvents->Fill(3.5);
289 if(bSelectVZ) fhNEvents->Fill(4.5);
292 //---------------------------------
294 //---------------------------------
296 //if(esdevent) bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esdevent, AliTriggerAnalysis::kV0AND);
297 AliVVZERO* v0 = fInputEvent->GetVZEROData();
298 bV0AND = ((v0->GetV0ADecision()==1) && (v0->GetV0CDecision()==1));
302 fhNEvents->Fill(5.5);
303 if (bSelectVZ) fhNEvents->Fill(6.5);
304 if (bSelectTrack) fhNEvents->Fill(7.5);
305 if (bSelectVZ && bSelectTrack) fhNEvents->Fill(8.5);
308 //---------------------------------
310 //---------------------------------
311 bPileup = event->IsPileupFromSPD(3, 0.8, 3., 2., 5.); //Default values, if not it does not compile
312 //bPileup = event->IsPileupFromSPD();
316 fhNEvents->Fill(9.5);
317 if(bV0AND) fhNEvents->Fill(16.5);
320 //---------------------------------
322 //---------------------------------
323 bGoodV = CheckForPrimaryVertex();
325 //Remove events with vertex (0,0,0), bad vertex reconstruction
326 if(TMath::Abs(v[0]) < 1.e-6 &&
327 TMath::Abs(v[1]) < 1.e-6 &&
328 TMath::Abs(v[2]) < 1.e-6) bGoodV = kFALSE;
332 fhXGoodVertex->Fill(v[0]);
333 fhYGoodVertex->Fill(v[1]);
334 fhZGoodVertex->Fill(v[2]);
336 fhNEvents->Fill(10.5);
337 if(bSelectVZ) fhNEvents->Fill(11.5);
338 if(bSelectTrack) fhNEvents->Fill(12.5);
339 if(bV0AND) fhNEvents->Fill(13.5);
340 if(bSelectVZ && bSelectTrack && bV0AND)
341 fhNEvents->Fill(14.5);
342 if(!bPileup) fhNEvents->Fill(15.5);
344 if(TMath::Abs(v[2]) < 10.)
346 if(InputEvent()->GetCentrality())
347 fhCentrality->Fill(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
349 if(InputEvent()->GetEventplane())
351 Float_t ep = InputEvent()->GetEventplane()->GetEventplane("V0", InputEvent());
353 ep+=TMath::Pi()/2.; // put same range as for <Q> method, [0,pi]
355 fhEventPlaneAngle->Fill(ep);
361 //printf("AliAnalysisTaskCounter::UserExec() : z vertex %d, good vertex %d, v0and %d, pile up %d, track mult %d\n ", bSelectVZ, bGoodV, bV0AND, bPileup, trackMult);
363 // Events that could be rejected in EMCAL
364 // LHC11a, SM4 and some SM3 events cut with this
365 Bool_t bEMCALRejected = kFALSE;
366 for (Int_t i = 0; i < InputEvent()->GetNumberOfCaloClusters(); i++)
368 AliVCluster *clus = InputEvent()->GetCaloCluster(i);
370 if ((clus->E() > 500 && clus->GetNCells() > 200 ) || clus->GetNCells() > 200)
373 //printf("Counter: Reject event with cluster: E %f, ncells %d\n",clus->E(),clus->GetNCells());
375 fhNEvents->Fill(17.5);
376 if(bSelectVZ) fhNEvents->Fill(18.5);
377 bEMCALRejected = kTRUE;
383 //LHC11a, 3 last runs, cut with this
386 // Count number of cells in SM3 with energy larger than 0.1, cut on this number
389 for(Int_t icell = 0; icell < event->GetEMCALCells()->GetNumberOfCells(); icell++)
391 if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
392 if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;
396 if(triggerclasses.Contains("EMC")) ncellcut = 35;
398 if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 )
400 //printf("Counter: reject event with ncells in SM3: ncells %d\n",ncells);
402 fhNEvents->Fill(19.5);
403 if(bSelectVZ) fhNEvents->Fill(20.5);
408 PostData(1,fOutputContainer);
413 //____________________________________________________
414 Bool_t AliAnalysisTaskCounter::CheckForPrimaryVertex()
416 //Check if the vertex was well reconstructed, copy of conversion group
418 AliESDEvent * esdevent = dynamic_cast<AliESDEvent*> (InputEvent());
419 AliAODEvent * aodevent = dynamic_cast<AliAODEvent*> (InputEvent());
423 if(esdevent->GetPrimaryVertex()->GetNContributors() > 0)
428 if(esdevent->GetPrimaryVertex()->GetNContributors() < 1)
431 if(esdevent->GetPrimaryVertexSPD()->GetNContributors() > 0)
436 if(esdevent->GetPrimaryVertexSPD()->GetNContributors() < 1)
444 if (aodevent->GetPrimaryVertex() != NULL)
446 if(aodevent->GetPrimaryVertex()->GetNContributors() > 0)
452 if(aodevent->GetPrimaryVertexSPD() != NULL)
454 if(aodevent->GetPrimaryVertexSPD()->GetNContributors() > 0)
460 AliWarning(Form("Number of contributors from bad vertex type:: %s",aodevent->GetPrimaryVertex()->GetName()));
472 //_____________________________________________________
473 void AliAnalysisTaskCounter::FinishTaskOutput()
475 // Put in the output some event summary histograms
477 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
478 AliInputEventHandler *inputH = dynamic_cast<AliInputEventHandler*>(am->GetInputEventHandler());
480 TH2F *histStat = dynamic_cast<TH2F*>(inputH->GetStatistics());
481 TH2F *histBin0 = dynamic_cast<TH2F*>(inputH->GetStatistics("BIN0"));
484 fOutputContainer->Add(histStat);
485 else if(DebugLevel() > 1)
486 printf("AliAnalysisTaskCounter::FinishTaskOutput() - Stat histogram not available check, \n if ESDs, that AliPhysicsSelection was on, \n if AODs, if EventStat_temp.root exists \n");
489 fOutputContainer->Add(histBin0);
494 //_____________________________________
495 Bool_t AliAnalysisTaskCounter::Notify()
498 // Implemented Notify() to read the cross sections
499 // and number of trials from pyxsec.root
502 if(!fCheckMCCrossSection) return kTRUE;
504 // Fetch the aod also from the input in,
505 // have todo it in notify
507 Float_t xsection = 0;
511 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
512 if(!tree) return kFALSE;
514 TFile *curfile = tree->GetCurrentFile();
516 if(!curfile) return kFALSE;
518 if(fCurrFileName == curfile->GetName()) return kFALSE;
520 fCurrFileName = TString(curfile->GetName());
522 if(!fh1Xsec||!fh1Trials)
524 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
528 Bool_t ok = PythiaInfoFromFile(fCurrFileName,xsection,trials);
530 if(!ok) return kFALSE;
532 fh1Xsec->Fill("<#sigma>",xsection);
534 // construct a poor man average trials
535 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
537 if(trials >= nEntries && nEntries > 0.) fAvgTrials = trials/nEntries;
539 fh1Trials->Fill("#sum{ntrials}",trials);
541 printf("AliAnalysisTaskCounter::Notify() - xs %f, trial %f, avg trials %f\n",xsection,trials, fAvgTrials);
543 if(fDebug) Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
548 //_____________________________________________________________________________________________________
549 Bool_t AliAnalysisTaskCounter::PythiaInfoFromFile(TString file,Float_t & xsec,Float_t & trials)
552 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
553 // This is to called in Notify and should provide the path to the AOD/ESD file
558 if(file.Contains("root_archive.zip#"))
560 Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
561 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
562 Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
563 file.Replace(pos+1,pos2-pos1,"");
567 // not an archive take the basename....
568 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
571 //Printf("%s",file.Data());
573 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
576 // next trial fetch the histgram file
577 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
580 // not a severe condition but inciate that we have no information
585 // find the tlist we want to be independtent of the name so use the Tkey
586 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
593 TList *list = dynamic_cast<TList*>(key->ReadObj());
600 xsec = ((TProfile*)list->FindObject("h1Xsec")) ->GetBinContent(1);
601 trials = ((TH1F*) list->FindObject("h1Trials"))->GetBinContent(1);
604 } // no tree pyxsec.root
607 TTree *xtree = (TTree*)fxsec->Get("Xsection");
615 Double_t xsection = 0;
616 xtree->SetBranchAddress("xsection",&xsection);
617 xtree->SetBranchAddress("ntrials",&ntrials);