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),
97 fMaxCosmicAngle(0.01),
99 fTrackRecEtaWindow(0.9),
109 fh1PtHardTrials(0x0),
110 fh1SelectionInfoESD(0x0),
111 fh1EventCutInfoESD(0),
115 fh2TriggerCount(0x0),
116 fh2ESDTriggerCount(0x0),
118 fh2ESDTriggerVtx(0x0),
119 fh2ESDTriggerRun(0x0),
121 fh1NCosmicsPerEvent(0x0),
135 fh2RPCentrality(0x0),
136 fh2RPACentrality(0x0),
137 fh2RPCCentrality(0x0),
138 fTriggerAnalysis(0x0),
141 fRunRange[0] = fRunRange[1] = 0;
144 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
145 AliAnalysisTaskSE(name),
146 fUseAODInput(kFALSE),
147 fUsePhysicsSelection(kFALSE),
149 fFilterAODCollisions(kFALSE),
150 fPhysicsSelectionFlag(AliVEvent::kMB),
151 fSelectionInfoESD(0),
155 fCollisionType(kPbPb),
164 fMaxCosmicAngle(0.01),
166 fTrackRecEtaWindow(0.9),
176 fh1PtHardTrials(0x0),
177 fh1SelectionInfoESD(0x0),
178 fh1EventCutInfoESD(0),
182 fh2TriggerCount(0x0),
183 fh2ESDTriggerCount(0x0),
185 fh2ESDTriggerVtx(0x0),
186 fh2ESDTriggerRun(0x0),
188 fh1NCosmicsPerEvent(0x0),
202 fh2RPCentrality(0x0),
203 fh2RPACentrality(0x0),
204 fh2RPCCentrality(0x0),
205 fTriggerAnalysis(0x0),
208 fRunRange[0] = fRunRange[1] = 0;
209 DefineOutput(1,TList::Class());
214 Bool_t AliAnalysisTaskJetServices::Notify()
217 // Implemented Notify() to read the cross sections
218 // and number of trials from pyxsec.root
221 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
222 Float_t xsection = 0;
227 TFile *curfile = tree->GetCurrentFile();
229 Error("Notify","No current file");
232 if(!fh1Xsec||!fh1Trials){
233 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
236 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
237 fh1Xsec->Fill("<#sigma>",xsection);
238 // construct a poor man average trials
239 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
240 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
245 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
249 // Create the output container
252 fRandomizer = new TRandom3(0);
253 if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
256 if(!fHistList)fHistList = new TList();
257 fHistList->SetOwner();
258 PostData(1, fHistList); // post data in any case once
260 Bool_t oldStatus = TH1::AddDirectoryStatus();
261 TH1::AddDirectory(kFALSE);
262 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
263 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
264 fHistList->Add(fh1Xsec);
266 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
267 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
268 fHistList->Add(fh1Trials);
270 const Int_t nBinPt = 100;
271 Double_t binLimitsPt[nBinPt+1];
272 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
274 binLimitsPt[iPt] = 0.0;
277 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
282 fh1CentralityESD = new TH1F("fh1CentralityESD","cent",103,-1,102);
283 fHistList->Add(fh1CentralityESD);
285 fh1Centrality = new TH1F("fh1Centrality","cent",103,-1,102);
286 fHistList->Add(fh1Centrality);
288 fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
289 fHistList->Add(fh1RP);
291 fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
292 fHistList->Add(fh2TriggerCount);
294 fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
295 fHistList->Add(fh2ESDTriggerCount);
296 const Int_t nBins = 6*kConstraints;
297 fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
298 fHistList->Add(fh2TriggerVtx);
300 fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
301 fHistList->Add(fh2ESDTriggerVtx);
304 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
305 fHistList->Add(fh1PtHard);
306 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
307 fHistList->Add(fh1PtHardTrials);
308 fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
309 AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
310 fHistList->Add(fh1SelectionInfoESD);
312 fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
313 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
314 fHistList->Add(fh1EventCutInfoESD);
316 // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
317 // 3 triggers BB BE/EB EE
319 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);
320 fHistList->Add(fh2ESDTriggerRun);
322 fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
323 fHistList->Add(fh2VtxXY);
324 // =========== Switch on Sumw2 for all histos ===========
325 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
326 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
331 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
335 fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
336 fHistList->Add(fh1NCosmicsPerEvent),
339 fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle from tracks" , 20, 0.,100., 180, 0, TMath::Pi());
340 fHistList->Add(fh2RPCentrality);
343 fh2RPACentrality = new TH2F("fh2RPACentrality" ,"Reaction Plane Angle from vzero A" , 20, 0.,100., 180, 0, TMath::Pi());
344 fHistList->Add(fh2RPACentrality);
346 fh2RPCCentrality = new TH2F("fh2RPCCentrality" ,"Reaction Plane Angle from vzero C" , 20, 0.,100., 180, 0, TMath::Pi());
347 fHistList->Add(fh2RPCCentrality);
349 fh2RPAC = new TH2F("fh2RPAC" ,"Reaction Plane Angle vzero a vs c" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
350 fHistList->Add( fh2RPAC);
352 fh2RPAT = new TH2F("fh2RPAT" ,"Reaction Plane Angle vzero a vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
353 fHistList->Add( fh2RPAT);
355 fh2RPCT = new TH2F("fh2RPCT" ,"Reaction Plane Angle vzero c vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
356 fHistList->Add( fh2RPCT);
358 fh2XYA = new TH2F("fh2XYA" ,"XY vzeroa ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
359 fHistList->Add(fh2XYA);
360 fh2XYC = new TH2F("fh2XYC" ,"XY vzeroc ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
361 fHistList->Add(fh2XYC);
364 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);
365 fHistList->Add(fp1RPXA);
367 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);
368 fHistList->Add(fp1RPYA);
371 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);
372 fHistList->Add(fp1RPXC);
374 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);
375 fHistList->Add(fp1RPYC);
378 TH1::AddDirectory(oldStatus);
380 // Add an AOD branch for replication
381 if(fNonStdFile.Length()){
382 if (fDebug > 1) AliInfo("Replicating header");
383 fgAODHeader = new AliAODHeader;
384 AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
385 if (fDebug > 1) AliInfo("Replicating vzeros");
386 fgAODVZERO = new AliAODVZERO;
387 AddAODBranch("AliAODVZERO",&fgAODVZERO,fNonStdFile.Data());
388 if (fDebug > 1) AliInfo("Replicating primary vertices");
389 fgAODVertices = new TClonesArray("AliAODVertex",3);
390 fgAODVertices->SetName("vertices");
391 AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
395 void AliAnalysisTaskJetServices::Init()
400 if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
404 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
408 // Execute analysis for current event
411 AliAODEvent *aod = 0;
412 AliESDEvent *esd = 0;
416 AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
417 AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
418 AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,0); // set slection to false
419 fSelectionInfoESD = 0; // reset
420 fEventCutInfoESD = 0; // reset
421 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
424 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
429 aod = dynamic_cast<AliAODEvent*>(InputEvent());
431 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
437 // assume that the AOD is in the general output...
440 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
443 esd = dynamic_cast<AliESDEvent*>(InputEvent());
447 Printf("Vertices %d",aod->GetNumberOfVertices());
448 Printf("tracks %d",aod->GetNumberOfTracks());
449 Printf("jets %d",aod->GetNJets());
451 fSelectionInfoESD |= kNoEventCut;
452 fEventCutInfoESD |= kNoEventCut;
454 Bool_t esdVtxValid = false;
455 Bool_t esdVtxIn = false;
456 Bool_t aodVtxValid = false;
457 Bool_t aodVtxIn = false;
461 if(!fTriggerAnalysis){
462 fTriggerAnalysis = new AliTriggerAnalysis;
463 fTriggerAnalysis->SetAnalyzeMC(fMC);
464 fTriggerAnalysis->SetSPDGFOThreshhold(1);
466 // fTriggerAnalysis->FillTriggerClasses(esd);
467 Bool_t v0A = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0A);
468 Bool_t v0C = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0C);
469 Bool_t v0ABG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0ABG);
470 Bool_t v0CBG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0CBG);
471 Bool_t spdFO = fTriggerAnalysis->SPDFiredChips(esd, 0);;
472 if(v0A)fSelectionInfoESD |= AliAnalysisHelperJetTasks::kV0A;
473 if(v0C)fSelectionInfoESD |= AliAnalysisHelperJetTasks::kV0C;
474 if(!(v0ABG||v0CBG))fSelectionInfoESD |= AliAnalysisHelperJetTasks::kNoV0BG;
475 if(spdFO)fSelectionInfoESD |= AliAnalysisHelperJetTasks::kSPDFO;
478 // Apply additional constraints
479 Bool_t esdEventSelected = IsEventSelected(esd);
480 Bool_t esdEventPileUp = IsEventPileUp(esd);
481 Bool_t esdEventCosmic = IsEventCosmic(esd);
483 Bool_t aodEventSelected = IsEventSelected(aod);
485 Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
488 fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
489 fh1EventCutInfoESD->Fill(fEventCutInfoESD);
491 if(esdEventSelected) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kVertexIn;
492 if(esdEventPileUp) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsPileUp;
493 if(esdEventCosmic) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsCosmic;
494 if(physicsSelection) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kPhysicsSelection;
497 // here we have all selection information, fill histogram
498 for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
499 if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
501 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD);
503 if(esd&&aod&&fDebug){
504 if(esdEventSelected&&!aodEventSelected){
505 Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
506 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
507 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
508 Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
510 Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
515 // loop over all possible triggers for esd
518 if(fCollisionType==kPbPb){
519 if(aod)cent = aod->GetHeader()->GetCentrality();
520 if(fDebug)Printf("%s:%d %3.3f",(char*)__FILE__,__LINE__,cent);
521 if(cent<0)cent = 101;
530 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
531 esdVtxValid = IsVertexValid(vtxESD);
532 esdVtxIn = IsVertexIn(vtxESD);
533 if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
534 if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
535 if(cent<=80&&esdVtxIn){
536 aodH->SetFillAOD(kTRUE);
537 aodH->SetFillExtension(kTRUE);
542 Float_t zvtx = vtxESD->GetZ();
543 Int_t iCl = GetEventClass(esd);
544 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
545 Bool_t cand = physicsSelection;
547 if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
548 fh2ESDTriggerCount->Fill(0.,kAllTriggered);
549 fh2ESDTriggerCount->Fill(iCl,kAllTriggered);
551 fh2ESDTriggerCount->Fill(0.,kSelectedALICE);
552 fh2ESDTriggerCount->Fill(iCl,kSelectedALICE);
553 fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
555 // if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
557 fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
558 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
559 fh2ESDTriggerVtx->Fill(iCl,zvtx);
561 fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
562 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
563 fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
566 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
567 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
568 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
572 if(cand&&esdVtxIn&&iCl<5){
573 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
574 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
575 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
576 fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
577 fh2ESDTriggerCount->Fill(iCl,kSelected);
578 fh2ESDTriggerCount->Fill(0.,kSelected);
579 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
580 if(esd->GetCentrality()){
581 Float_t tmpCent = 100;
582 tmpCent = esd->GetCentrality()->GetCentralityPercentile("V0M");
583 if(tmpCent<0)tmpCent = 101;
584 fh1CentralityESD->Fill(tmpCent);
592 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
593 aodVtxValid = IsVertexValid(vtxAOD);
594 aodVtxIn = IsVertexIn(vtxAOD);
595 Float_t zvtx = vtxAOD->GetZ();
596 Int_t iCl = GetEventClass(aod);
597 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
598 Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
599 if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
600 fh2TriggerCount->Fill(0.,kAllTriggered);
601 fh2TriggerCount->Fill(iCl,kAllTriggered);
603 fh2TriggerCount->Fill(0.,kSelectedALICE);
604 fh2TriggerCount->Fill(iCl,kSelectedALICE);
605 fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
608 fh2TriggerCount->Fill(0.,kTriggeredVertex);
609 fh2TriggerCount->Fill(iCl,kTriggeredVertex);
610 fh2TriggerVtx->Fill(iCl,zvtx);
612 fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
613 fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
614 fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
617 fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
618 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
619 fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
623 if(cand&&aodVtxIn&&iCl<5){
624 fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
625 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
626 fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
627 fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
628 fh2TriggerCount->Fill(iCl,kSelected);
629 fh2TriggerCount->Fill(0.,kSelected);
630 fh1Centrality->Fill(cent);
631 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
632 if(aodH&&cand&&fFilterAODCollisions&&!esd){
633 if(fCentrality<=80&&aodVtxIn){
634 aodH->SetFillAOD(kTRUE);
635 aodH->SetFillExtension(kTRUE);
640 GetListOfTracks(&recTracks);
641 CalculateReactionPlaneAngleVZERO(aod);
642 fRPAngle = aod->GetHeader()->GetEventplane();
643 fh1RP->Fill(fRPAngle);
644 fh2RPCentrality->Fill(fCentrality,fRPAngle);
645 fh2RPACentrality->Fill(fCentrality,fPsiVZEROA);
646 fh2RPCCentrality->Fill(fCentrality,fPsiVZEROC);
647 fh2RPAC->Fill(fPsiVZEROA,fPsiVZEROC);
648 fh2RPAT->Fill(fPsiVZEROA,fRPAngle);
649 fh2RPCT->Fill(fPsiVZEROC,fRPAngle);
650 if(fRPMethod==kRPTracks)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fRPAngle); // set slection to false
651 else if(fRPMethod==kRPVZEROA)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
652 else if(fRPMethod==kRPVZEROC)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
654 if(fUseAODInput&&fCentrality<=80){
655 if(fFilterAODCollisions&&aod){
656 aodH->SetFillAOD(kTRUE);
657 aodH->SetFillExtension(kTRUE);
663 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
667 Double_t nTrials = 1; // Trials for MC trigger
669 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
670 AliMCEvent* mcEvent = MCEvent();
671 // AliStack *pStack = 0;
673 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
675 nTrials = pythiaGenHeader->Trials();
676 ptHard = pythiaGenHeader->GetPtHard();
677 int iProcessType = pythiaGenHeader->ProcessType();
679 // 12 f+barf -> f+barf
684 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
685 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
686 fh1PtHard->Fill(ptHard);
687 fh1PtHardTrials->Fill(ptHard,nTrials);
689 }// if pythia gen header
695 if(fNonStdFile.Length()&&aod){
697 *fgAODHeader = *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
698 Double_t q[kRPMethods] = {fRPAngle,fPsiVZEROA,fPsiVZEROC};
699 fgAODHeader->SetQTheta(q,kRPMethods);
702 *fgAODVZERO = *(dynamic_cast<AliAODVZERO*>(aod->GetVZEROData()));
705 fgAODVertices->Delete();
706 TClonesArray &vertices = *fgAODVertices;
707 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
708 // we only use some basic information,
714 vtxAOD->GetXYZ(pos); // position
715 vtxAOD->GetCovMatrix(covVtx); //covariance matrix
717 AliAODVertex * vtx = new(vertices[jVertices++])
718 AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
719 vtx->SetName(vtxAOD->GetName());
720 vtx->SetTitle(vtxAOD->GetTitle());
721 TString vtitle = vtxAOD->GetTitle();
722 vtx->SetNContributors(vtxAOD->GetNContributors());
724 // Add SPD "main" vertex
725 const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
726 vtxS->GetXYZ(pos); // position
727 vtxS->GetCovMatrix(covVtx); //covariance matrix
728 AliAODVertex * mVSPD = new(vertices[jVertices++])
729 AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
730 mVSPD->SetName(vtxS->GetName());
731 mVSPD->SetTitle(vtxS->GetTitle());
732 mVSPD->SetNContributors(vtxS->GetNContributors());
734 // Add tpc only vertex
736 const AliESDVertex *vtxT = esd->GetPrimaryVertexTPC();
737 vtxT->GetXYZ(pos); // position
738 vtxT->GetCovMatrix(covVtx); //covariance matrix
739 AliAODVertex * mVTPC = new(vertices[jVertices++])
740 AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
741 mVTPC->SetName(vtxT->GetName());
742 mVTPC->SetTitle(vtxT->GetTitle());
743 mVTPC->SetNContributors(vtxT->GetNContributors());
748 PostData(1, fHistList);
751 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
752 if(!esd)return kFALSE;
753 const AliESDVertex *vtx = esd->GetPrimaryVertex();
754 return IsVertexIn(vtx); // vertex in calls vertex valid
757 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
759 fgAODVertices->Delete();
760 delete fgAODVertices;
763 if(fgAODHeader)delete fgAODHeader;
764 if(fgAODVZERO)delete fgAODVZERO;
773 void AliAnalysisTaskJetServices::SetV0Centroids(TProfile *xa,TProfile *ya,TProfile *xc, TProfile *yc){
776 if(fp1CalibRPXA)delete fp1CalibRPXA;
777 fp1CalibRPXA = (TProfile*)xa->Clone(Form("%sCalib",xa->GetName()));
780 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
783 if(fp1CalibRPYA)delete fp1CalibRPYA;
784 fp1CalibRPYA = (TProfile*)ya->Clone(Form("%sCalib",ya->GetName()));
787 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
790 if(fp1CalibRPXC)delete fp1CalibRPXC;
791 fp1CalibRPXC = (TProfile*)xc->Clone(Form("%sCalib",xc->GetName()));
794 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
797 if(fp1CalibRPYC)delete fp1CalibRPYC;
798 fp1CalibRPYC = (TProfile*)yc->Clone(Form("%sCalib",yc->GetName()));
801 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
805 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
806 if(!aod)return kFALSE;
807 const AliAODVertex *vtx = aod->GetPrimaryVertex();
808 return IsVertexIn(vtx); // VertexIn calls VertexValid
811 Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
813 // We can check the type of the vertex by its name and title
814 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
816 // Tracks PrimaryVertex VertexerTracksNoConstraint
817 // Tracks PrimaryVertex VertexerTracksWithConstraint
818 // TPC TPCVertex VertexerTracksNoConstraint
819 // TPC TPCVertex VertexerTracksWithConstraint
820 // SPD SPDVertex vertexer: 3D
821 // SPD SPDVertex vertexer: Z
824 Printf("%s:%d No ESD vertex found",(char*)__FILE__,__LINE__);
827 Int_t nCont = vtx->GetNContributors();
829 fEventCutInfoESD |= kContributorsCut1;
831 fEventCutInfoESD |= kContributorsCut2;
833 fEventCutInfoESD |= kContributorsCut3;
838 if(nCont<3)return kFALSE;
840 // do not want tpc only primary vertex
841 TString vtxName(vtx->GetName());
842 if(vtxName.Contains("TPCVertex")){
843 fEventCutInfoESD |= kVertexTPC;
846 if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
847 if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
850 TString vtxTitle(vtx->GetTitle());
851 if(vtxTitle.Contains("vertexer: Z")){
852 if(vtx->GetDispersion()>0.02)return kFALSE;
854 fEventCutInfoESD |= kSPDDispersionCut;
859 Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
861 // We can check the type of the vertex by its name and title
862 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
864 // Tracks PrimaryVertex VertexerTracksNoConstraint
865 // TPC TPCVertex VertexerTracksNoConstraint
866 // SPD SPDVertex vertexer: 3D
867 // SPD SPDVertex vertexer: Z
870 Printf("%s:%d No AOD vertex found",(char*)__FILE__,__LINE__);
876 Printf(" n contrib %d",vtx->GetNContributors());
880 // if(vtx->GetNContributors()<3)return kFALSE;
881 // do not want tpc only primary vertex
882 TString vtxName(vtx->GetName());
883 if(vtxName.Contains("TPCVertex"))return kFALSE;
885 // no dispersion yet...
887 TString vtxTitle(vtx->GetTitle());
888 if(vtxTitle.Contains("vertexer: Z")){
889 if(vtx->GetDispersion()>0.02)return kFALSE;
896 Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
898 if(!IsVertexValid(vtx))return kFALSE;
900 Float_t zvtx = vtx->GetZ();
901 Float_t yvtx = vtx->GetY();
902 Float_t xvtx = vtx->GetX();
910 if(TMath::Abs(zvtx)>fVtxZCut){
913 fEventCutInfoESD |= kVertexZCut;
914 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
915 if(r2>(fVtxRCut*fVtxRCut)){
918 fEventCutInfoESD |= kVertexRCut;
923 Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
925 if(!IsVertexValid(vtx))return kFALSE;
927 Float_t zvtx = vtx->GetZ();
928 Float_t yvtx = vtx->GetY();
929 Float_t xvtx = vtx->GetX();
935 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
937 Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
941 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
942 if(!esd)return kFALSE;
943 return esd->IsPileupFromSPD();
946 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
947 if(!esd)return kFALSE;
948 // add track cuts for which we look for cosmics...
950 Bool_t isCosmic = kFALSE;
951 Int_t nTracks = esd->GetNumberOfTracks();
952 Int_t nCosmicCandidates = 0;
954 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
955 AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
956 if (!track1) continue;
957 UInt_t status1 = track1->GetStatus();
958 //If track is ITS stand alone track, skip the track
959 if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
960 if(track1->Pt()<fPtMinCosmic) continue;
961 //Start 2nd track loop to look for correlations
962 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
963 AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
964 if(!track2) continue;
965 UInt_t status2 = track2->GetStatus();
966 //If track is ITS stand alone track, skip the track
967 if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
968 if(track2->Pt()<fPtMinCosmic) continue;
969 //Check if back-to-back
970 Double_t mom1[3],mom2[3];
971 track1->GetPxPyPz(mom1);
972 track2->GetPxPyPz(mom2);
973 TVector3 momv1(mom1[0],mom1[1],mom1[2]);
974 TVector3 momv2(mom2[0],mom2[1],mom2[2]);
975 Float_t theta = (float)(momv1.Phi()-momv2.Phi());
976 if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
978 Float_t deltaPhi = track1->Phi()-track2->Phi();
979 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
981 Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
982 if(rIsol<fRIsolMinCosmic) continue;
984 if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
985 nCosmicCandidates+=1;
992 fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
998 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
1001 if(fCollisionType==kPbPb){
1002 if(esd->GetCentrality()){
1003 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
1005 if(cent<0)cent = 101;
1006 if(cent>80||cent<0)return 5;
1007 if(cent>50)return 4;
1008 if(cent>30)return 3;
1009 if(cent>10)return 2;
1016 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
1018 if(fCollisionType==kPbPb){
1019 Float_t cent = aod->GetHeader()->GetCentrality();
1020 if(cent>80||cent<0)return 5;
1021 if(cent>50)return 4;
1022 if(cent>30)return 3;
1023 if(cent>10)return 2;
1031 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
1033 // Terminate analysis
1037 Bool_t AliAnalysisTaskJetServices::CalculateReactionPlaneAngleVZERO(AliAODEvent *aod){
1039 // 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};
1040 if(!aod)return kFALSE;
1041 static bool bFirst = true;
1042 static Double_t v0phi[64] = {0,};
1046 for(int iArm = 0; iArm<2; iArm++){
1047 for(int iRing = 0; iRing<4 ; iRing++){
1048 for(int iSec = 0; iSec<8 ; iSec++){
1049 v0phi[is] = 22.5 + 45. * iSec;
1050 v0phi[is] *= TMath::Pi()/180;
1051 // cout<< is <<" "<< v0phi[is]<<endl;
1060 const AliAODVZERO *aodVZERO = aod->GetVZEROData();
1061 Double_t numYZNA = 0,numXZNA = 0,sumZNA = 0;
1062 Double_t meanXA = 0,meanYA = 0;
1064 Double_t numYZNC = 0,numXZNC = 0,sumZNC = 0;
1065 Double_t meanXC = 0,meanYC = 0;
1069 static Int_t iOldRun = -1;
1070 static Int_t iFoundBin = -1;
1072 if(aod->GetRunNumber()!=iOldRun&&(fp1CalibRPYA)){
1073 // search only or the bin in case of new runs
1075 Int_t ib = fp1CalibRPYA->FindBin(aod->GetRunNumber());
1076 Float_t err = fp1CalibRPYA->GetBinError(ib);
1077 if(err>0){// value can be zero...
1083 while(iFoundBin<0&&(ibLo>0||ibUp<=fp1CalibRPYA->GetNbinsX())){
1084 err = fp1CalibRPYA->GetBinError(ibLo);
1089 err = fp1CalibRPYA->GetBinError(ibUp);
1090 if(err>0)iFoundBin = ibUp;
1096 iOldRun = aod->GetRunNumber();
1099 if(fDebug)Printf("%s:%d iFoundBin %d",(char*)__FILE__,__LINE__,iFoundBin);
1101 if(iFoundBin>0&&(fp1CalibRPYA)){
1102 meanXA = fp1CalibRPXA->GetBinContent(iFoundBin);
1103 meanYA = fp1CalibRPYA->GetBinContent(iFoundBin);
1104 meanXC = fp1CalibRPXC->GetBinContent(iFoundBin);
1105 meanYC = fp1CalibRPYC->GetBinContent(iFoundBin);
1108 if(fDebug)Printf("%s:%d iFoundBin %1.3E %1.3E %1.3E %1.3E",(char*)__FILE__,__LINE__,meanXA,meanYA,meanXC,meanYC);
1110 for (int i=0; i<64; i++) {
1111 Double_t mult = aodVZERO->GetMultiplicity(i);
1112 Double_t phi = v0phi[i];
1114 if (i<32) { //C-side
1115 Double_t wZNC= mult;
1116 numYZNC += sin(2.*phi)*wZNC;
1117 numXZNC += cos(2.*phi)*wZNC;
1120 else if(i>31){ //A-side
1122 numYZNA += sin(2.*phi)*wZNA;
1123 numXZNA += cos(2.*phi)*wZNA;
1137 XC = numXZNC/sumZNC;
1138 YC = numYZNC/sumZNC;
1139 fPsiVZEROC = 0.5*TMath::ATan2(YC-meanYC, XC-meanXC);
1140 if(fPsiVZEROC>TMath::Pi()){fPsiVZEROC-=TMath::Pi();}
1141 if(fPsiVZEROC<0){fPsiVZEROC+=TMath::Pi();}
1144 fh2XYC->Fill(XC-meanXC,YC-meanYC); // control
1145 fp1RPXC->Fill(aod->GetRunNumber(),XC);
1146 fp1RPYC->Fill(aod->GetRunNumber(),YC);
1152 XA = numXZNA/sumZNA;
1153 YA = numYZNA/sumZNA;
1154 fPsiVZEROA = 0.5*TMath::ATan2(YA-meanYA, XA-meanXA);
1155 if(fPsiVZEROA>TMath::Pi()){fPsiVZEROA-=TMath::Pi();}
1156 if(fPsiVZEROA<0){fPsiVZEROA+=TMath::Pi();}
1158 fh2XYA->Fill(XA-meanXA,YA-meanYA); // control
1159 fp1RPXA->Fill(aod->GetRunNumber(),XA);
1160 fp1RPYA->Fill(aod->GetRunNumber(),YA);
1167 Int_t AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
1169 AliAODEvent *aod = 0;
1170 if(fUseAODInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1171 else aod = AODEvent();
1175 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1176 AliAODTrack *tr = aod->GetTrack(it);
1177 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1178 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1179 if(tr->Pt()<fMinTrackPt)continue;