2 // **************************************
3 // Task used for the correction of determiantion of reconstructed jet spectra
4 // Compares input (gen) and output (rec) jets
5 // *******************************************
8 /**************************************************************************
9 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
11 * Author: The ALICE Off-line Project. *
12 * Contributors are mentioned in the code where appropriate. *
14 * Permission to use, copy, modify and distribute this software and its *
15 * documentation strictly for non-commercial purposes is hereby granted *
16 * without fee, provided that the above copyright notice appears in all *
17 * copies and that both the copyright notice and this permission notice *
18 * appear in the supporting documentation. The authors make no claims *
19 * about the suitability of this software for any purpose. It is *
20 * provided "as is" without express or implied warranty. *
21 **************************************************************************/
28 #include <TInterpreter.h>
37 #include <TLorentzVector.h>
39 #include <TClonesArray.h>
40 #include "TDatabasePDG.h"
42 #include "AliAnalysisTaskJetServices.h"
43 #include "AliCentrality.h"
44 #include "AliAnalysisDataContainer.h"
45 #include "AliAnalysisDataSlot.h"
46 #include "AliAnalysisManager.h"
47 #include "AliJetFinder.h"
48 #include "AliJetHeader.h"
49 #include "AliJetReader.h"
50 #include "AliJetReaderHeader.h"
51 #include "AliUA1JetHeaderV1.h"
52 #include "AliESDEvent.h"
53 #include "AliAODEvent.h"
54 #include "AliAODHandler.h"
55 #include "AliAODTrack.h"
56 #include "AliAODVZERO.h"
57 #include "AliAODJet.h"
58 #include "AliAODMCParticle.h"
59 #include "AliMCEventHandler.h"
60 #include "AliMCEvent.h"
62 #include "AliGenPythiaEventHeader.h"
63 #include "AliJetKineReaderHeader.h"
64 #include "AliGenCocktailEventHeader.h"
65 #include "AliInputEventHandler.h"
66 #include "AliPhysicsSelection.h"
67 #include "AliTriggerAnalysis.h"
69 #include "AliAnalysisHelperJetTasks.h"
71 ClassImp(AliAnalysisTaskJetServices)
73 AliAODVZERO* AliAnalysisTaskJetServices::fgAODVZERO = NULL;
74 AliAODHeader* AliAnalysisTaskJetServices::fgAODHeader = NULL;
75 TClonesArray* AliAnalysisTaskJetServices::fgAODVertices = NULL;
77 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices():
80 fUsePhysicsSelection(kFALSE),
82 fFilterAODCollisions(kFALSE),
83 fPhysicsSelectionFlag(AliVEvent::kMB),
88 fCollisionType(kPbPb),
98 fMaxCosmicAngle(0.01),
100 fTrackRecEtaWindow(0.9),
113 fh1PtHardTrials(0x0),
114 fh1SelectionInfoESD(0x0),
115 fh1EventCutInfoESD(0),
119 fh2ReducedTrigger(0),
120 fh2CentralityTriggerESD(0),
121 fh2CentralityTrigger(0),
122 fh2TriggerCount(0x0),
123 fh2ESDTriggerCount(0x0),
125 fh2ESDTriggerVtx(0x0),
126 fh2ESDTriggerRun(0x0),
128 fh1NCosmicsPerEvent(0x0),
142 fh2RPCentrality(0x0),
143 fh2RPACentrality(0x0),
144 fh2RPCCentrality(0x0),
145 fTriggerAnalysis(0x0),
148 fRunRange[0] = fRunRange[1] = 0;
151 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
152 AliAnalysisTaskSE(name),
153 fUseAODInput(kFALSE),
154 fUsePhysicsSelection(kFALSE),
156 fFilterAODCollisions(kFALSE),
157 fPhysicsSelectionFlag(AliVEvent::kMB),
158 fSelectionInfoESD(0),
162 fCollisionType(kPbPb),
172 fMaxCosmicAngle(0.01),
174 fTrackRecEtaWindow(0.9),
187 fh1PtHardTrials(0x0),
188 fh1SelectionInfoESD(0x0),
189 fh1EventCutInfoESD(0),
193 fh2ReducedTrigger(0),
194 fh2CentralityTriggerESD(0),
195 fh2CentralityTrigger(0),
196 fh2TriggerCount(0x0),
197 fh2ESDTriggerCount(0x0),
199 fh2ESDTriggerVtx(0x0),
200 fh2ESDTriggerRun(0x0),
202 fh1NCosmicsPerEvent(0x0),
216 fh2RPCentrality(0x0),
217 fh2RPACentrality(0x0),
218 fh2RPCCentrality(0x0),
219 fTriggerAnalysis(0x0),
222 fRunRange[0] = fRunRange[1] = 0;
223 DefineOutput(1,TList::Class());
228 Bool_t AliAnalysisTaskJetServices::Notify()
231 // Implemented Notify() to read the cross sections
232 // and number of trials from pyxsec.root
235 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
236 Float_t xsection = 0;
241 TFile *curfile = tree->GetCurrentFile();
243 Error("Notify","No current file");
246 if(!fh1Xsec||!fh1Trials){
247 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
250 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
251 fh1Xsec->Fill("<#sigma>",xsection);
252 // construct a poor man average trials
253 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
254 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
255 fh1Trials->Fill("#sum{avg ntrials}",ftrials);
260 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
264 // Create the output container
267 fRandomizer = new TRandom3(0);
268 if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
271 if(!fHistList)fHistList = new TList();
272 fHistList->SetOwner();
273 PostData(1, fHistList); // post data in any case once
275 Bool_t oldStatus = TH1::AddDirectoryStatus();
276 TH1::AddDirectory(kFALSE);
277 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
278 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
279 fHistList->Add(fh1Xsec);
281 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
282 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
283 fHistList->Add(fh1Trials);
285 fh1AvgTrials = new TH1F("fh1AvgTrials","avg trials root file",1,0,1);
286 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
287 fHistList->Add(fh1AvgTrials);
289 const Int_t nBinPt = 125;
290 Double_t binLimitsPt[nBinPt+1];
291 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
293 binLimitsPt[iPt] = 0.0;
296 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
301 fh1CentralityESD = new TH1F("fh1CentralityESD","cent",103,-1,102);
302 fHistList->Add(fh1CentralityESD);
304 fh1Centrality = new TH1F("fh1Centrality","cent",103,-1,102);
305 fHistList->Add(fh1Centrality);
307 fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
308 fHistList->Add(fh1RP);
310 fh2ReducedTrigger = new TH2F("fh2ReducedTrigger","red trigger;red trigger;centrality",1<<fNTrigger,-0.5,(1<<fNTrigger)-0.5,100,-0.5,99.5);
311 fHistList->Add(fh2ReducedTrigger);
313 fh2CentralityTriggerESD = new TH2F("fh2CentralityTriggerESD",";cent;trigger no",103,-1,102,fNTrigger,-0.5,fNTrigger-0.5);
314 fHistList->Add(fh2CentralityTriggerESD);
316 fh2CentralityTrigger = new TH2F("fh2CentralityTrigger",";cent;trigger no",103,-1,102,fNTrigger,-0.5,fNTrigger-0.5);
317 fHistList->Add(fh2CentralityTrigger);
319 for(int i = 0;i<fNTrigger;++i){
320 fh2CentralityTriggerESD->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
321 fh2CentralityTrigger->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
325 fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
326 fHistList->Add(fh2TriggerCount);
328 fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
329 fHistList->Add(fh2ESDTriggerCount);
330 const Int_t nBins = 6*kConstraints;
331 fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
332 fHistList->Add(fh2TriggerVtx);
334 fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
335 fHistList->Add(fh2ESDTriggerVtx);
338 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
339 fHistList->Add(fh1PtHard);
340 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
341 fHistList->Add(fh1PtHardTrials);
342 fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
343 AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
344 fHistList->Add(fh1SelectionInfoESD);
346 fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
347 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
348 fHistList->Add(fh1EventCutInfoESD);
350 // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
351 // 3 triggers BB BE/EB EE
353 fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Eventclass vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
354 fHistList->Add(fh2ESDTriggerRun);
356 fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
357 fHistList->Add(fh2VtxXY);
358 // =========== Switch on Sumw2 for all histos ===========
359 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
360 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
365 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
369 fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
370 fHistList->Add(fh1NCosmicsPerEvent),
373 fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle from tracks" , 20, 0.,100., 180, 0, TMath::Pi());
374 fHistList->Add(fh2RPCentrality);
377 fh2RPACentrality = new TH2F("fh2RPACentrality" ,"Reaction Plane Angle from vzero A" , 20, 0.,100., 180, 0, TMath::Pi());
378 fHistList->Add(fh2RPACentrality);
380 fh2RPCCentrality = new TH2F("fh2RPCCentrality" ,"Reaction Plane Angle from vzero C" , 20, 0.,100., 180, 0, TMath::Pi());
381 fHistList->Add(fh2RPCCentrality);
383 fh2RPAC = new TH2F("fh2RPAC" ,"Reaction Plane Angle vzero a vs c" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
384 fHistList->Add( fh2RPAC);
386 fh2RPAT = new TH2F("fh2RPAT" ,"Reaction Plane Angle vzero a vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
387 fHistList->Add( fh2RPAT);
389 fh2RPCT = new TH2F("fh2RPCT" ,"Reaction Plane Angle vzero c vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
390 fHistList->Add( fh2RPCT);
392 fh2XYA = new TH2F("fh2XYA" ,"XY vzeroa ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
393 fHistList->Add(fh2XYA);
394 fh2XYC = new TH2F("fh2XYC" ,"XY vzeroc ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
395 fHistList->Add(fh2XYC);
398 fp1RPXA = new TProfile("fp1RPXA","mean vzeroa x vs run number;run;x",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
399 fHistList->Add(fp1RPXA);
401 fp1RPYA = new TProfile("fp1RPYA","mean vzeroa y vs run number;run;y",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
402 fHistList->Add(fp1RPYA);
405 fp1RPXC = new TProfile("fp1RPXC","mean vzeroc x vs run number;run;x",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
406 fHistList->Add(fp1RPXC);
408 fp1RPYC = new TProfile("fp1RPYC","mean vzeroa y vs run number;run;y",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
409 fHistList->Add(fp1RPYC);
412 TH1::AddDirectory(oldStatus);
414 // Add an AOD branch for replication
415 if(fNonStdFile.Length()){
416 if (fDebug > 1) AliInfo("Replicating header");
417 fgAODHeader = new AliAODHeader;
418 AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
419 if (fDebug > 1) AliInfo("Replicating vzeros");
420 fgAODVZERO = new AliAODVZERO;
421 AddAODBranch("AliAODVZERO",&fgAODVZERO,fNonStdFile.Data());
422 if (fDebug > 1) AliInfo("Replicating primary vertices");
423 fgAODVertices = new TClonesArray("AliAODVertex",3);
424 fgAODVertices->SetName("vertices");
425 AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
429 void AliAnalysisTaskJetServices::Init()
434 if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
438 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
442 // Execute analysis for current event
445 AliAODEvent *aod = 0;
446 AliESDEvent *esd = 0;
450 AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
451 AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
452 AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,0); // set slection to false
453 fSelectionInfoESD = 0; // reset
454 fEventCutInfoESD = 0; // reset
455 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
458 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
463 aod = dynamic_cast<AliAODEvent*>(InputEvent());
465 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
471 // assume that the AOD is in the general output...
474 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
477 esd = dynamic_cast<AliESDEvent*>(InputEvent());
481 Printf("Vertices %d",aod->GetNumberOfVertices());
482 Printf("tracks %d",aod->GetNumberOfTracks());
483 Printf("jets %d",aod->GetNJets());
485 fSelectionInfoESD |= kNoEventCut;
486 fEventCutInfoESD |= kNoEventCut;
488 Bool_t esdVtxValid = false;
489 Bool_t esdVtxIn = false;
490 Bool_t aodVtxValid = false;
491 Bool_t aodVtxIn = false;
493 // Apply additional constraints
494 Bool_t esdEventSelected = IsEventSelected(esd);
495 Bool_t esdEventPileUp = IsEventPileUp(esd);
496 Bool_t esdEventCosmic = IsEventCosmic(esd);
498 Bool_t aodEventSelected = IsEventSelected(aod);
500 Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
503 fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
504 fh1EventCutInfoESD->Fill(fEventCutInfoESD);
506 if(esdEventSelected) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kVertexIn;
507 if(esdEventPileUp) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsPileUp;
508 if(esdEventCosmic) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsCosmic;
509 if(physicsSelection) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kPhysicsSelection;
514 // here we have all selection information, fill histogram
515 for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
516 if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
518 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD);
520 if(esd&&aod&&fDebug){
521 if(esdEventSelected&&!aodEventSelected){
522 Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
523 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
524 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
525 Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
527 Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
532 // loop over all possible triggers for esd
535 if(fCollisionType==kPbPb){
536 if(aod)cent = aod->GetHeader()->GetCentrality();
537 if(fDebug)Printf("%s:%d %3.3f",(char*)__FILE__,__LINE__,cent);
538 if(cent<0)cent = 101;
547 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
548 esdVtxValid = IsVertexValid(vtxESD);
549 esdVtxIn = IsVertexIn(vtxESD);
550 if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
551 if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
552 if(cent<=80&&esdVtxIn){
553 aodH->SetFillAOD(kTRUE);
554 aodH->SetFillExtension(kTRUE);
559 Float_t zvtx = vtxESD->GetZ();
560 Int_t iCl = GetEventClass(esd);
561 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
562 Bool_t cand = physicsSelection;
564 if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
565 fh2ESDTriggerCount->Fill(0.,kAllTriggered);
566 fh2ESDTriggerCount->Fill(iCl,kAllTriggered);
568 fh2ESDTriggerCount->Fill(0.,kSelectedALICE);
569 fh2ESDTriggerCount->Fill(iCl,kSelectedALICE);
570 fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
572 // if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
574 fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
575 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
576 fh2ESDTriggerVtx->Fill(iCl,zvtx);
578 fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
579 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
580 fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
583 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
584 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
585 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
589 if(cand&&esdVtxIn&&iCl<5){
590 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
591 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
592 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
593 fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
594 fh2ESDTriggerCount->Fill(iCl,kSelected);
595 fh2ESDTriggerCount->Fill(0.,kSelected);
596 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
597 if(esd->GetCentrality()){
598 Float_t tmpCent = 100;
599 tmpCent = esd->GetCentrality()->GetCentralityPercentile("V0M");
600 if(tmpCent<0)tmpCent = 101;
601 fh1CentralityESD->Fill(tmpCent);
603 for(int it = 0;it<fNTrigger;++it){
604 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
605 fh2CentralityTriggerESD->Fill(tmpCent,it);
609 fh2ReducedTrigger->Fill(ir,cent);
617 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
618 aodVtxValid = IsVertexValid(vtxAOD);
619 aodVtxIn = IsVertexIn(vtxAOD);
620 Float_t zvtx = vtxAOD->GetZ();
621 Int_t iCl = GetEventClass(aod);
622 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
623 Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
624 if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
625 fh2TriggerCount->Fill(0.,kAllTriggered);
626 fh2TriggerCount->Fill(iCl,kAllTriggered);
628 fh2TriggerCount->Fill(0.,kSelectedALICE);
629 fh2TriggerCount->Fill(iCl,kSelectedALICE);
630 fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
633 fh2TriggerCount->Fill(0.,kTriggeredVertex);
634 fh2TriggerCount->Fill(iCl,kTriggeredVertex);
635 fh2TriggerVtx->Fill(iCl,zvtx);
637 fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
638 fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
639 fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
642 fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
643 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
644 fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
648 if(cand&&aodVtxIn&&iCl<5){
649 fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
650 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
651 fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
652 fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
653 fh2TriggerCount->Fill(iCl,kSelected);
654 fh2TriggerCount->Fill(0.,kSelected);
655 fh1Centrality->Fill(cent);
656 for(int it = 0;it<fNTrigger;++it){
657 if(fInputHandler->IsEventSelected()&fTriggerBit[it])fh2CentralityTrigger->Fill(cent,it);
659 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
660 if(aodH&&cand&&fFilterAODCollisions&&!esd){
661 if(fCentrality<=80&&aodVtxIn){
662 aodH->SetFillAOD(kTRUE);
663 aodH->SetFillExtension(kTRUE);
668 GetListOfTracks(&recTracks);
669 CalculateReactionPlaneAngleVZERO(aod);
670 fRPAngle = aod->GetHeader()->GetEventplane();
671 fh1RP->Fill(fRPAngle);
672 fh2RPCentrality->Fill(fCentrality,fRPAngle);
673 fh2RPACentrality->Fill(fCentrality,fPsiVZEROA);
674 fh2RPCCentrality->Fill(fCentrality,fPsiVZEROC);
675 fh2RPAC->Fill(fPsiVZEROA,fPsiVZEROC);
676 fh2RPAT->Fill(fPsiVZEROA,fRPAngle);
677 fh2RPCT->Fill(fPsiVZEROC,fRPAngle);
678 if(fRPMethod==kRPTracks)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fRPAngle); // set slection to false
679 else if(fRPMethod==kRPVZEROA)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
680 else if(fRPMethod==kRPVZEROC)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
682 if(fUseAODInput&&fCentrality<=80){
683 if(fFilterAODCollisions&&aod){
684 aodH->SetFillAOD(kTRUE);
685 aodH->SetFillExtension(kTRUE);
691 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
695 Double_t nTrials = 1; // Trials for MC trigger
697 fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
698 AliMCEvent* mcEvent = MCEvent();
699 // AliStack *pStack = 0;
701 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
703 nTrials = pythiaGenHeader->Trials();
704 ptHard = pythiaGenHeader->GetPtHard();
705 int iProcessType = pythiaGenHeader->ProcessType();
707 // 12 f+barf -> f+barf
712 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
713 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
714 fh1PtHard->Fill(ptHard);
715 fh1PtHardTrials->Fill(ptHard,nTrials);
717 }// if pythia gen header
723 if(fNonStdFile.Length()&&aod){
725 *fgAODHeader = *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
726 Double_t q[kRPMethods] = {fRPAngle,fPsiVZEROA,fPsiVZEROC};
727 fgAODHeader->SetQTheta(q,kRPMethods);
730 *fgAODVZERO = *(dynamic_cast<AliAODVZERO*>(aod->GetVZEROData()));
733 fgAODVertices->Delete();
734 TClonesArray &vertices = *fgAODVertices;
735 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
736 // we only use some basic information,
742 vtxAOD->GetXYZ(pos); // position
743 vtxAOD->GetCovMatrix(covVtx); //covariance matrix
745 AliAODVertex * vtx = new(vertices[jVertices++])
746 AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
747 vtx->SetName(vtxAOD->GetName());
748 vtx->SetTitle(vtxAOD->GetTitle());
749 TString vtitle = vtxAOD->GetTitle();
750 vtx->SetNContributors(vtxAOD->GetNContributors());
752 // Add SPD "main" vertex
753 const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
754 vtxS->GetXYZ(pos); // position
755 vtxS->GetCovMatrix(covVtx); //covariance matrix
756 AliAODVertex * mVSPD = new(vertices[jVertices++])
757 AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
758 mVSPD->SetName(vtxS->GetName());
759 mVSPD->SetTitle(vtxS->GetTitle());
760 mVSPD->SetNContributors(vtxS->GetNContributors());
762 // Add tpc only vertex
764 const AliESDVertex *vtxT = esd->GetPrimaryVertexTPC();
765 vtxT->GetXYZ(pos); // position
766 vtxT->GetCovMatrix(covVtx); //covariance matrix
767 AliAODVertex * mVTPC = new(vertices[jVertices++])
768 AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
769 mVTPC->SetName(vtxT->GetName());
770 mVTPC->SetTitle(vtxT->GetTitle());
771 mVTPC->SetNContributors(vtxT->GetNContributors());
776 PostData(1, fHistList);
779 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
780 if(!esd)return kFALSE;
781 const AliESDVertex *vtx = esd->GetPrimaryVertex();
782 return IsVertexIn(vtx); // vertex in calls vertex valid
785 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
787 fgAODVertices->Delete();
788 delete fgAODVertices;
791 if(fgAODHeader)delete fgAODHeader;
792 if(fgAODVZERO)delete fgAODVZERO;
797 delete [] fTriggerBit;
798 delete [] fTriggerName;
803 void AliAnalysisTaskJetServices::SetV0Centroids(TProfile *xa,TProfile *ya,TProfile *xc, TProfile *yc){
806 if(fp1CalibRPXA)delete fp1CalibRPXA;
807 fp1CalibRPXA = (TProfile*)xa->Clone(Form("%sCalib",xa->GetName()));
810 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
813 if(fp1CalibRPYA)delete fp1CalibRPYA;
814 fp1CalibRPYA = (TProfile*)ya->Clone(Form("%sCalib",ya->GetName()));
817 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
820 if(fp1CalibRPXC)delete fp1CalibRPXC;
821 fp1CalibRPXC = (TProfile*)xc->Clone(Form("%sCalib",xc->GetName()));
824 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
827 if(fp1CalibRPYC)delete fp1CalibRPYC;
828 fp1CalibRPYC = (TProfile*)yc->Clone(Form("%sCalib",yc->GetName()));
831 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
835 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
836 if(!aod)return kFALSE;
837 const AliAODVertex *vtx = aod->GetPrimaryVertex();
838 return IsVertexIn(vtx); // VertexIn calls VertexValid
841 Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
843 // We can check the type of the vertex by its name and title
844 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
846 // Tracks PrimaryVertex VertexerTracksNoConstraint
847 // Tracks PrimaryVertex VertexerTracksWithConstraint
848 // TPC TPCVertex VertexerTracksNoConstraint
849 // TPC TPCVertex VertexerTracksWithConstraint
850 // SPD SPDVertex vertexer: 3D
851 // SPD SPDVertex vertexer: Z
854 Printf("%s:%d No ESD vertex found",(char*)__FILE__,__LINE__);
857 Int_t nCont = vtx->GetNContributors();
859 fEventCutInfoESD |= kContributorsCut1;
861 fEventCutInfoESD |= kContributorsCut2;
863 fEventCutInfoESD |= kContributorsCut3;
868 if(nCont<3)return kFALSE;
870 // do not want tpc only primary vertex
871 TString vtxName(vtx->GetName());
872 if(vtxName.Contains("TPCVertex")){
873 fEventCutInfoESD |= kVertexTPC;
876 if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
877 if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
880 TString vtxTitle(vtx->GetTitle());
881 if(vtxTitle.Contains("vertexer: Z")){
882 if(vtx->GetDispersion()>0.02)return kFALSE;
884 fEventCutInfoESD |= kSPDDispersionCut;
889 Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
891 // We can check the type of the vertex by its name and title
892 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
894 // Tracks PrimaryVertex VertexerTracksNoConstraint
895 // TPC TPCVertex VertexerTracksNoConstraint
896 // SPD SPDVertex vertexer: 3D
897 // SPD SPDVertex vertexer: Z
900 Printf("%s:%d No AOD vertex found",(char*)__FILE__,__LINE__);
904 TString vtxName(vtx->GetName());
906 Printf(" n contrib %d",vtx->GetNContributors());
907 Printf("vtxName: %s",vtxName.Data());
911 // if(vtx->GetNContributors()<3)return kFALSE;
912 // do not want tpc only primary vertex
914 if(vtxName.Contains("TPCVertex"))return kFALSE;
916 // no dispersion yet...
918 TString vtxTitle(vtx->GetTitle());
919 if(vtxTitle.Contains("vertexer: Z")){
920 if(vtx->GetDispersion()>0.02)return kFALSE;
927 Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
929 if(!IsVertexValid(vtx))return kFALSE;
931 Float_t zvtx = vtx->GetZ();
932 Float_t yvtx = vtx->GetY();
933 Float_t xvtx = vtx->GetX();
941 if(TMath::Abs(zvtx)>fVtxZCut){
944 fEventCutInfoESD |= kVertexZCut;
945 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
946 if(r2>(fVtxRCut*fVtxRCut)){
949 fEventCutInfoESD |= kVertexRCut;
954 Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
956 if(!IsVertexValid(vtx))return kFALSE;
958 Float_t zvtx = vtx->GetZ();
959 Float_t yvtx = vtx->GetY();
960 Float_t xvtx = vtx->GetX();
966 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
968 Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
972 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
973 if(!esd)return kFALSE;
974 return esd->IsPileupFromSPD();
977 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
978 if(!esd)return kFALSE;
979 // add track cuts for which we look for cosmics...
981 Bool_t isCosmic = kFALSE;
982 Int_t nTracks = esd->GetNumberOfTracks();
983 Int_t nCosmicCandidates = 0;
985 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
986 AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
987 if (!track1) continue;
988 UInt_t status1 = track1->GetStatus();
989 //If track is ITS stand alone track, skip the track
990 if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
991 if(track1->Pt()<fPtMinCosmic) continue;
992 //Start 2nd track loop to look for correlations
993 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
994 AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
995 if(!track2) continue;
996 UInt_t status2 = track2->GetStatus();
997 //If track is ITS stand alone track, skip the track
998 if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
999 if(track2->Pt()<fPtMinCosmic) continue;
1000 //Check if back-to-back
1001 Double_t mom1[3],mom2[3];
1002 track1->GetPxPyPz(mom1);
1003 track2->GetPxPyPz(mom2);
1004 TVector3 momv1(mom1[0],mom1[1],mom1[2]);
1005 TVector3 momv2(mom2[0],mom2[1],mom2[2]);
1006 Float_t theta = (float)(momv1.Phi()-momv2.Phi());
1007 if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
1009 Float_t deltaPhi = track1->Phi()-track2->Phi();
1010 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
1012 Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
1013 if(rIsol<fRIsolMinCosmic) continue;
1015 if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
1016 nCosmicCandidates+=1;
1023 fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
1029 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
1032 if(fCollisionType==kPbPb){
1033 if(esd->GetCentrality()){
1034 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
1036 if(cent<0)cent = 101;
1037 if(cent>80||cent<0)return 5;
1038 if(cent>50)return 4;
1039 if(cent>30)return 3;
1040 if(cent>10)return 2;
1047 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
1049 if(fCollisionType==kPbPb){
1050 Float_t cent = aod->GetHeader()->GetCentrality();
1051 if(cent>80||cent<0)return 5;
1052 if(cent>50)return 4;
1053 if(cent>30)return 3;
1054 if(cent>10)return 2;
1062 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
1064 // Terminate analysis
1068 Bool_t AliAnalysisTaskJetServices::CalculateReactionPlaneAngleVZERO(AliAODEvent *aod){
1070 // const Double_t arr_eta[11]={-3.7, -3.2, -2.7, -2.2, -1.7,0, 2.8, 3.4, 3.9, 4.5,5.1};
1071 if(!aod)return kFALSE;
1072 static bool bFirst = true;
1073 static Double_t v0phi[64] = {0,};
1077 for(int iArm = 0; iArm<2; iArm++){
1078 for(int iRing = 0; iRing<4 ; iRing++){
1079 for(int iSec = 0; iSec<8 ; iSec++){
1080 v0phi[is] = 22.5 + 45. * iSec;
1081 v0phi[is] *= TMath::Pi()/180;
1082 // cout<< is <<" "<< v0phi[is]<<endl;
1091 const AliAODVZERO *aodVZERO = aod->GetVZEROData();
1092 Double_t numYZNA = 0,numXZNA = 0,sumZNA = 0;
1093 Double_t meanXA = 0,meanYA = 0;
1095 Double_t numYZNC = 0,numXZNC = 0,sumZNC = 0;
1096 Double_t meanXC = 0,meanYC = 0;
1100 static Int_t iOldRun = -1;
1101 static Int_t iFoundBin = -1;
1103 if(aod->GetRunNumber()!=iOldRun&&(fp1CalibRPYA)){
1104 // search only or the bin in case of new runs
1106 Int_t ib = fp1CalibRPYA->FindBin(aod->GetRunNumber());
1107 Float_t err = fp1CalibRPYA->GetBinError(ib);
1108 if(err>0){// value can be zero...
1114 while(iFoundBin<0&&(ibLo>0||ibUp<=fp1CalibRPYA->GetNbinsX())){
1115 err = fp1CalibRPYA->GetBinError(ibLo);
1120 err = fp1CalibRPYA->GetBinError(ibUp);
1121 if(err>0)iFoundBin = ibUp;
1127 iOldRun = aod->GetRunNumber();
1130 if(fDebug)Printf("%s:%d iFoundBin %d",(char*)__FILE__,__LINE__,iFoundBin);
1132 if(iFoundBin>0&&(fp1CalibRPYA)){
1133 meanXA = fp1CalibRPXA->GetBinContent(iFoundBin);
1134 meanYA = fp1CalibRPYA->GetBinContent(iFoundBin);
1135 meanXC = fp1CalibRPXC->GetBinContent(iFoundBin);
1136 meanYC = fp1CalibRPYC->GetBinContent(iFoundBin);
1139 if(fDebug)Printf("%s:%d iFoundBin %1.3E %1.3E %1.3E %1.3E",(char*)__FILE__,__LINE__,meanXA,meanYA,meanXC,meanYC);
1141 for (int i=0; i<64; i++) {
1142 Double_t mult = aodVZERO->GetMultiplicity(i);
1143 Double_t phi = v0phi[i];
1145 if (i<32) { //C-side
1146 Double_t wZNC= mult;
1147 numYZNC += sin(2.*phi)*wZNC;
1148 numXZNC += cos(2.*phi)*wZNC;
1151 else if(i>31){ //A-side
1153 numYZNA += sin(2.*phi)*wZNA;
1154 numXZNA += cos(2.*phi)*wZNA;
1168 XC = numXZNC/sumZNC;
1169 YC = numYZNC/sumZNC;
1170 fPsiVZEROC = 0.5*TMath::ATan2(YC-meanYC, XC-meanXC);
1171 if(fPsiVZEROC>TMath::Pi()){fPsiVZEROC-=TMath::Pi();}
1172 if(fPsiVZEROC<0){fPsiVZEROC+=TMath::Pi();}
1175 fh2XYC->Fill(XC-meanXC,YC-meanYC); // control
1176 fp1RPXC->Fill(aod->GetRunNumber(),XC);
1177 fp1RPYC->Fill(aod->GetRunNumber(),YC);
1183 XA = numXZNA/sumZNA;
1184 YA = numYZNA/sumZNA;
1185 fPsiVZEROA = 0.5*TMath::ATan2(YA-meanYA, XA-meanXA);
1186 if(fPsiVZEROA>TMath::Pi()){fPsiVZEROA-=TMath::Pi();}
1187 if(fPsiVZEROA<0){fPsiVZEROA+=TMath::Pi();}
1189 fh2XYA->Fill(XA-meanXA,YA-meanYA); // control
1190 fp1RPXA->Fill(aod->GetRunNumber(),XA);
1191 fp1RPYA->Fill(aod->GetRunNumber(),YA);
1198 Int_t AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
1200 AliAODEvent *aod = 0;
1201 if(fUseAODInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1202 else aod = AODEvent();
1206 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1207 AliAODTrack *tr = aod->GetTrack(it);
1208 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1209 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1210 if(tr->Pt()<fMinTrackPt)continue;
1219 void AliAnalysisTaskJetServices::SetNTrigger(Int_t n){
1222 delete [] fTriggerName;
1223 fTriggerName = new TString [fNTrigger];
1224 delete [] fTriggerBit;fTriggerBit = 0;
1225 fTriggerBit = new UInt_t [fNTrigger];
1232 void AliAnalysisTaskJetServices::SetTrigger(Int_t i,UInt_t it,const char* c){
1237 Printf("%p",&fTriggerName[i]);
1239 fTriggerBit[i] = it;
1240 fTriggerName[i] = c; // placement new
1241 // new(&fTriggerName[i]) TString(c); // placement new
1242 Printf("%s",fTriggerName[i].Data());