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),
99 fMaxCosmicAngle(0.01),
101 fTrackRecEtaWindow(0.9),
114 fh1PtHardTrials(0x0),
115 fh1SelectionInfoESD(0x0),
116 fh1EventCutInfoESD(0),
120 fh2ReducedTrigger(0),
121 fh2CentralityTriggerESD(0),
122 fh2CentralityTrigger(0),
123 fh2TriggerCount(0x0),
124 fh2ESDTriggerCount(0x0),
126 fh2ESDTriggerVtx(0x0),
127 fh2ESDTriggerRun(0x0),
129 fh1NCosmicsPerEvent(0x0),
143 fh2RPCentrality(0x0),
144 fh2RPACentrality(0x0),
145 fh2RPCCentrality(0x0),
146 fTriggerAnalysis(0x0),
149 fRunRange[0] = fRunRange[1] = 0;
152 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
153 AliAnalysisTaskSE(name),
154 fUseAODInput(kFALSE),
155 fUsePhysicsSelection(kFALSE),
157 fFilterAODCollisions(kFALSE),
158 fPhysicsSelectionFlag(AliVEvent::kMB),
159 fSelectionInfoESD(0),
163 fCollisionType(kPbPb),
174 fMaxCosmicAngle(0.01),
176 fTrackRecEtaWindow(0.9),
189 fh1PtHardTrials(0x0),
190 fh1SelectionInfoESD(0x0),
191 fh1EventCutInfoESD(0),
195 fh2ReducedTrigger(0),
196 fh2CentralityTriggerESD(0),
197 fh2CentralityTrigger(0),
198 fh2TriggerCount(0x0),
199 fh2ESDTriggerCount(0x0),
201 fh2ESDTriggerVtx(0x0),
202 fh2ESDTriggerRun(0x0),
204 fh1NCosmicsPerEvent(0x0),
218 fh2RPCentrality(0x0),
219 fh2RPACentrality(0x0),
220 fh2RPCCentrality(0x0),
221 fTriggerAnalysis(0x0),
224 fRunRange[0] = fRunRange[1] = 0;
225 DefineOutput(1,TList::Class());
230 Bool_t AliAnalysisTaskJetServices::Notify()
233 // Implemented Notify() to read the cross sections
234 // and number of trials from pyxsec.root
237 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
238 Float_t xsection = 0;
243 TFile *curfile = tree->GetCurrentFile();
245 Error("Notify","No current file");
248 if(!fh1Xsec||!fh1Trials){
249 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
252 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
253 fh1Xsec->Fill("<#sigma>",xsection);
254 // construct a poor man average trials
255 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
256 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
257 fh1Trials->Fill("#sum{avg ntrials}",ftrials);
262 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
266 // Create the output container
269 fRandomizer = new TRandom3(0);
270 if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
273 if(!fHistList)fHistList = new TList();
274 fHistList->SetOwner();
275 PostData(1, fHistList); // post data in any case once
277 Bool_t oldStatus = TH1::AddDirectoryStatus();
278 TH1::AddDirectory(kFALSE);
279 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
280 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
281 fHistList->Add(fh1Xsec);
283 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
284 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
285 fHistList->Add(fh1Trials);
287 fh1AvgTrials = new TH1F("fh1AvgTrials","avg trials root file",1,0,1);
288 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
289 fHistList->Add(fh1AvgTrials);
291 const Int_t nBinPt = 125;
292 Double_t binLimitsPt[nBinPt+1];
293 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
295 binLimitsPt[iPt] = 0.0;
298 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
303 fh1CentralityESD = new TH1F("fh1CentralityESD","cent",103,-1,102);
304 fHistList->Add(fh1CentralityESD);
306 fh1Centrality = new TH1F("fh1Centrality","cent",103,-1,102);
307 fHistList->Add(fh1Centrality);
309 fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
310 fHistList->Add(fh1RP);
312 fh2ReducedTrigger = new TH2F("fh2ReducedTrigger","red trigger;red trigger;centrality",1<<fNTrigger,-0.5,(1<<fNTrigger)-0.5,100,-0.5,99.5);
313 fHistList->Add(fh2ReducedTrigger);
315 fh2CentralityTriggerESD = new TH2F("fh2CentralityTriggerESD",";cent;trigger no",103,-1,102,fNTrigger,-0.5,fNTrigger-0.5);
316 fHistList->Add(fh2CentralityTriggerESD);
318 fh2CentralityTrigger = new TH2F("fh2CentralityTrigger",";cent;trigger no",103,-1,102,fNTrigger,-0.5,fNTrigger-0.5);
319 fHistList->Add(fh2CentralityTrigger);
321 for(int i = 0;i<fNTrigger;++i){
322 fh2CentralityTriggerESD->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
323 fh2CentralityTrigger->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
327 fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
328 fHistList->Add(fh2TriggerCount);
330 fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
331 fHistList->Add(fh2ESDTriggerCount);
332 const Int_t nBins = 6*kConstraints;
333 fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
334 fHistList->Add(fh2TriggerVtx);
336 fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
337 fHistList->Add(fh2ESDTriggerVtx);
340 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
341 fHistList->Add(fh1PtHard);
342 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
343 fHistList->Add(fh1PtHardTrials);
344 fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
345 AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
346 fHistList->Add(fh1SelectionInfoESD);
348 fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
349 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
350 fHistList->Add(fh1EventCutInfoESD);
352 // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
353 // 3 triggers BB BE/EB EE
355 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);
356 fHistList->Add(fh2ESDTriggerRun);
358 fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
359 fHistList->Add(fh2VtxXY);
360 // =========== Switch on Sumw2 for all histos ===========
361 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
362 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
367 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
371 fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
372 fHistList->Add(fh1NCosmicsPerEvent),
375 fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle from tracks" , 20, 0.,100., 180, 0, TMath::Pi());
376 fHistList->Add(fh2RPCentrality);
379 fh2RPACentrality = new TH2F("fh2RPACentrality" ,"Reaction Plane Angle from vzero A" , 20, 0.,100., 180, 0, TMath::Pi());
380 fHistList->Add(fh2RPACentrality);
382 fh2RPCCentrality = new TH2F("fh2RPCCentrality" ,"Reaction Plane Angle from vzero C" , 20, 0.,100., 180, 0, TMath::Pi());
383 fHistList->Add(fh2RPCCentrality);
385 fh2RPAC = new TH2F("fh2RPAC" ,"Reaction Plane Angle vzero a vs c" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
386 fHistList->Add( fh2RPAC);
388 fh2RPAT = new TH2F("fh2RPAT" ,"Reaction Plane Angle vzero a vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
389 fHistList->Add( fh2RPAT);
391 fh2RPCT = new TH2F("fh2RPCT" ,"Reaction Plane Angle vzero c vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
392 fHistList->Add( fh2RPCT);
394 fh2XYA = new TH2F("fh2XYA" ,"XY vzeroa ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
395 fHistList->Add(fh2XYA);
396 fh2XYC = new TH2F("fh2XYC" ,"XY vzeroc ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
397 fHistList->Add(fh2XYC);
400 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);
401 fHistList->Add(fp1RPXA);
403 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);
404 fHistList->Add(fp1RPYA);
407 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);
408 fHistList->Add(fp1RPXC);
410 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);
411 fHistList->Add(fp1RPYC);
414 TH1::AddDirectory(oldStatus);
416 // Add an AOD branch for replication
417 if(fNonStdFile.Length()){
418 if (fDebug > 1) AliInfo("Replicating header");
419 fgAODHeader = new AliAODHeader;
420 AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
421 if (fDebug > 1) AliInfo("Replicating vzeros");
422 fgAODVZERO = new AliAODVZERO;
423 AddAODBranch("AliAODVZERO",&fgAODVZERO,fNonStdFile.Data());
424 if (fDebug > 1) AliInfo("Replicating primary vertices");
425 fgAODVertices = new TClonesArray("AliAODVertex",3);
426 fgAODVertices->SetName("vertices");
427 AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
431 void AliAnalysisTaskJetServices::Init()
436 if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
440 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
444 // Execute analysis for current event
447 AliAODEvent *aod = 0;
448 AliESDEvent *esd = 0;
452 AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
453 AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
454 AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,0); // set slection to false
455 fSelectionInfoESD = 0; // reset
456 fEventCutInfoESD = 0; // reset
457 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
460 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
465 aod = dynamic_cast<AliAODEvent*>(InputEvent());
467 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
473 // assume that the AOD is in the general output...
476 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
479 esd = dynamic_cast<AliESDEvent*>(InputEvent());
483 Printf("Vertices %d",aod->GetNumberOfVertices());
484 Printf("tracks %d",aod->GetNumberOfTracks());
485 Printf("jets %d",aod->GetNJets());
487 fSelectionInfoESD |= kNoEventCut;
488 fEventCutInfoESD |= kNoEventCut;
490 Bool_t esdVtxValid = false;
491 Bool_t esdVtxIn = false;
492 Bool_t aodVtxValid = false;
493 Bool_t aodVtxIn = false;
495 // Apply additional constraints
496 Bool_t esdEventSelected = IsEventSelected(esd);
497 Bool_t esdEventPileUp = IsEventPileUp(esd);
498 Bool_t esdEventCosmic = IsEventCosmic(esd);
500 Bool_t aodEventSelected = IsEventSelected(aod);
502 Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
505 fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
506 fh1EventCutInfoESD->Fill (fEventCutInfoESD);
508 if(esdEventSelected) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kVertexIn;
509 if(esdEventPileUp) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsPileUp;
510 if(esdEventCosmic) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsCosmic;
511 if(physicsSelection) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kPhysicsSelection;
516 // here we have all selection information, fill histogram
517 for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
518 if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
520 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD);
522 if(esd&&aod&&fDebug){
523 if(esdEventSelected&&!aodEventSelected){
524 Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
525 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
526 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
527 Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
529 Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
534 // loop over all possible triggers for esd
537 if(fCollisionType==kPbPb){
538 if(aod)cent = aod->GetHeader()->GetCentrality();
539 if(fDebug)Printf("%s:%d %3.3f",(char*)__FILE__,__LINE__,cent);
540 if(cent<0)cent = 101;
549 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
550 const AliESDVertex *vtxESDSPD = esd->GetPrimaryVertexSPD();
551 esdVtxValid = IsVertexValid(vtxESD);
552 esdVtxIn = IsVertexIn(vtxESD,vtxESDSPD);
553 if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
554 if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
555 if(cent<=80&&esdVtxIn){
556 aodH->SetFillAOD(kTRUE);
557 aodH->SetFillExtension(kTRUE);
562 Float_t zvtx = vtxESD->GetZ();
563 Int_t iCl = GetEventClass(esd);
564 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
565 Bool_t cand = physicsSelection;
567 if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
568 fh2ESDTriggerCount->Fill(0.,kAllTriggered);
569 fh2ESDTriggerCount->Fill(iCl,kAllTriggered);
571 fh2ESDTriggerCount->Fill(0.,kSelectedALICE);
572 fh2ESDTriggerCount->Fill(iCl,kSelectedALICE);
573 fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
575 // if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
577 fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
578 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
579 fh2ESDTriggerVtx->Fill(iCl,zvtx);
581 fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
582 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
583 fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
586 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
587 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
588 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
592 if(cand&&esdVtxIn&&iCl<5){
593 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
594 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
595 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
596 fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
597 fh2ESDTriggerCount->Fill(iCl,kSelected);
598 fh2ESDTriggerCount->Fill(0.,kSelected);
599 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
600 if(esd->GetCentrality()){
601 Float_t tmpCent = 100;
602 tmpCent = esd->GetCentrality()->GetCentralityPercentile("V0M");
603 if(tmpCent<0)tmpCent = 101;
604 fh1CentralityESD->Fill(tmpCent);
606 for(int it = 0;it<fNTrigger;++it){
607 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
608 fh2CentralityTriggerESD->Fill(tmpCent,it);
612 fh2ReducedTrigger->Fill(ir,cent);
620 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
621 const AliAODVertex *vtxAODSPD = aod->GetPrimaryVertexSPD();
622 aodVtxValid = IsVertexValid(vtxAOD);
623 aodVtxIn = IsVertexIn(vtxAOD,vtxAODSPD);
624 Float_t zvtx = vtxAOD->GetZ();
625 Int_t iCl = GetEventClass(aod);
626 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
627 Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
628 if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
629 fh2TriggerCount->Fill(0.,kAllTriggered);
630 fh2TriggerCount->Fill(iCl,kAllTriggered);
632 fh2TriggerCount->Fill(0.,kSelectedALICE);
633 fh2TriggerCount->Fill(iCl,kSelectedALICE);
634 fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
637 fh2TriggerCount->Fill(0.,kTriggeredVertex);
638 fh2TriggerCount->Fill(iCl,kTriggeredVertex);
639 fh2TriggerVtx->Fill(iCl,zvtx);
641 fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
642 fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
643 fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
646 fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
647 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
648 fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
652 if(cand&&aodVtxIn&&iCl<5){
653 fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
654 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
655 fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
656 fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
657 fh2TriggerCount->Fill(iCl,kSelected);
658 fh2TriggerCount->Fill(0.,kSelected);
659 fh1Centrality->Fill(cent);
660 for(int it = 0;it<fNTrigger;++it){
661 if(fInputHandler->IsEventSelected()&fTriggerBit[it])fh2CentralityTrigger->Fill(cent,it);
663 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
664 if(aodH&&cand&&fFilterAODCollisions&&!esd){
665 if(fCentrality<=80&&aodVtxIn){
666 aodH->SetFillAOD(kTRUE);
667 aodH->SetFillExtension(kTRUE);
672 GetListOfTracks(&recTracks);
673 CalculateReactionPlaneAngleVZERO(aod);
674 fRPAngle = aod->GetHeader()->GetEventplane();
675 fh1RP->Fill(fRPAngle);
676 fh2RPCentrality->Fill(fCentrality,fRPAngle);
677 fh2RPACentrality->Fill(fCentrality,fPsiVZEROA);
678 fh2RPCCentrality->Fill(fCentrality,fPsiVZEROC);
679 fh2RPAC->Fill(fPsiVZEROA,fPsiVZEROC);
680 fh2RPAT->Fill(fPsiVZEROA,fRPAngle);
681 fh2RPCT->Fill(fPsiVZEROC,fRPAngle);
682 if(fRPMethod==kRPTracks)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fRPAngle); // set slection to false
683 else if(fRPMethod==kRPVZEROA)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
684 else if(fRPMethod==kRPVZEROC)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
686 if(fUseAODInput&&fCentrality<=80){
687 if(fFilterAODCollisions&&aod){
688 aodH->SetFillAOD(kTRUE);
689 aodH->SetFillExtension(kTRUE);
695 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
699 Double_t nTrials = 1; // Trials for MC trigger
701 fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
702 AliMCEvent* mcEvent = MCEvent();
703 // AliStack *pStack = 0;
705 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
707 nTrials = pythiaGenHeader->Trials();
708 ptHard = pythiaGenHeader->GetPtHard();
709 int iProcessType = pythiaGenHeader->ProcessType();
711 // 12 f+barf -> f+barf
716 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
717 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
718 fh1PtHard->Fill(ptHard);
719 fh1PtHardTrials->Fill(ptHard,nTrials);
721 }// if pythia gen header
727 if(fNonStdFile.Length()&&aod){
729 *fgAODHeader = *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
730 Double_t q[kRPMethods] = {fRPAngle,fPsiVZEROA,fPsiVZEROC};
731 fgAODHeader->SetQTheta(q,kRPMethods);
734 *fgAODVZERO = *(dynamic_cast<AliAODVZERO*>(aod->GetVZEROData()));
737 fgAODVertices->Delete();
738 TClonesArray &vertices = *fgAODVertices;
739 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
740 // we only use some basic information,
746 vtxAOD->GetXYZ(pos); // position
747 vtxAOD->GetCovMatrix(covVtx); //covariance matrix
749 AliAODVertex * vtx = new(vertices[jVertices++])
750 AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
751 vtx->SetName(vtxAOD->GetName());
752 vtx->SetTitle(vtxAOD->GetTitle());
753 TString vtitle = vtxAOD->GetTitle();
754 vtx->SetNContributors(vtxAOD->GetNContributors());
756 // Add SPD "main" vertex
757 const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
758 vtxS->GetXYZ(pos); // position
759 vtxS->GetCovMatrix(covVtx); //covariance matrix
760 AliAODVertex * mVSPD = new(vertices[jVertices++])
761 AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
762 mVSPD->SetName(vtxS->GetName());
763 mVSPD->SetTitle(vtxS->GetTitle());
764 mVSPD->SetNContributors(vtxS->GetNContributors());
766 // Add tpc only vertex
768 const AliESDVertex *vtxT = esd->GetPrimaryVertexTPC();
769 vtxT->GetXYZ(pos); // position
770 vtxT->GetCovMatrix(covVtx); //covariance matrix
771 AliAODVertex * mVTPC = new(vertices[jVertices++])
772 AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
773 mVTPC->SetName(vtxT->GetName());
774 mVTPC->SetTitle(vtxT->GetTitle());
775 mVTPC->SetNContributors(vtxT->GetNContributors());
780 PostData(1, fHistList);
783 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
784 if(!esd)return kFALSE;
785 const AliESDVertex *vtx = esd->GetPrimaryVertex();
786 const AliESDVertex *vtxSPD = esd->GetPrimaryVertex();
787 return IsVertexIn(vtx,vtxSPD); // vertex in calls vertex valid
790 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
792 fgAODVertices->Delete();
793 delete fgAODVertices;
796 if(fgAODHeader)delete fgAODHeader;
797 if(fgAODVZERO)delete fgAODVZERO;
802 delete [] fTriggerBit;
803 delete [] fTriggerName;
808 void AliAnalysisTaskJetServices::SetV0Centroids(TProfile *xa,TProfile *ya,TProfile *xc, TProfile *yc){
811 if(fp1CalibRPXA)delete fp1CalibRPXA;
812 fp1CalibRPXA = (TProfile*)xa->Clone(Form("%sCalib",xa->GetName()));
815 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
818 if(fp1CalibRPYA)delete fp1CalibRPYA;
819 fp1CalibRPYA = (TProfile*)ya->Clone(Form("%sCalib",ya->GetName()));
822 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
825 if(fp1CalibRPXC)delete fp1CalibRPXC;
826 fp1CalibRPXC = (TProfile*)xc->Clone(Form("%sCalib",xc->GetName()));
829 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
832 if(fp1CalibRPYC)delete fp1CalibRPYC;
833 fp1CalibRPYC = (TProfile*)yc->Clone(Form("%sCalib",yc->GetName()));
836 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
840 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
841 if(!aod)return kFALSE;
842 const AliAODVertex *vtx = aod->GetPrimaryVertex();
843 const AliAODVertex *vtxSPD = aod->GetPrimaryVertexSPD();
844 return IsVertexIn(vtx,vtxSPD); // VertexIn calls VertexValid
847 Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
849 // We can check the type of the vertex by its name and title
850 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
852 // Tracks PrimaryVertex VertexerTracksNoConstraint
853 // Tracks PrimaryVertex VertexerTracksWithConstraint
854 // TPC TPCVertex VertexerTracksNoConstraint
855 // TPC TPCVertex VertexerTracksWithConstraint
856 // SPD SPDVertex vertexer: 3D
857 // SPD SPDVertex vertexer: Z
860 Printf("%s:%d No ESD vertex found",(char*)__FILE__,__LINE__);
863 Int_t nCont = vtx->GetNContributors();
865 fEventCutInfoESD |= kContributorsCut1;
867 fEventCutInfoESD |= kContributorsCut2;
869 fEventCutInfoESD |= kContributorsCut3;
874 if(nCont<3)return kFALSE;
876 // do not want tpc only primary vertex
877 TString vtxName(vtx->GetName());
878 if(vtxName.Contains("TPCVertex")){
879 fEventCutInfoESD |= kVertexTPC;
882 if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
883 if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
886 TString vtxTitle(vtx->GetTitle());
887 if(vtxTitle.Contains("vertexer: Z")){
888 if(vtx->GetDispersion()>0.02)return kFALSE;
890 fEventCutInfoESD |= kSPDDispersionCut;
895 Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
897 // We can check the type of the vertex by its name and title
898 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
900 // Tracks PrimaryVertex VertexerTracksNoConstraint
901 // TPC TPCVertex VertexerTracksNoConstraint
902 // SPD SPDVertex vertexer: 3D
903 // SPD SPDVertex vertexer: Z
906 Printf("%s:%d No AOD vertex found",(char*)__FILE__,__LINE__);
910 TString vtxName(vtx->GetName());
912 Printf(" n contrib %d",vtx->GetNContributors());
913 Printf("vtxName: %s",vtxName.Data());
917 // if(vtx->GetNContributors()<3)return kFALSE;
918 // do not want tpc only primary vertex
920 if(vtxName.Contains("TPCVertex"))return kFALSE;
922 // no dispersion yet...
924 TString vtxTitle(vtx->GetTitle());
925 if(vtxTitle.Contains("vertexer: Z")){
926 if(vtx->GetDispersion()>0.02)return kFALSE;
933 Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx,const AliESDVertex* vtxSPD) {
935 if(!IsVertexValid(vtx))return kFALSE;
937 Float_t zvtx = vtx->GetZ();
938 Float_t yvtx = vtx->GetY();
939 Float_t xvtx = vtx->GetX();
940 Float_t deltaZ = zvtx - vtxSPD->GetZ();
941 if(fDebug)Printf("%s:%d deltaz %3.3f ",__FILE__,__LINE__,deltaZ);
949 if(TMath::Abs(zvtx)>fVtxZCut){
952 fEventCutInfoESD |= kVertexZCut;
953 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
954 if(r2>(fVtxRCut*fVtxRCut)){
957 fEventCutInfoESD |= kVertexRCut;
962 if(TMath::Abs(deltaZ)>fVtxDeltaZCut)return kFALSE;
963 fEventCutInfoESD |= kVertexDeltaZCut;
968 Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx,const AliAODVertex* vtxSPD) const {
970 if(!IsVertexValid(vtx))return kFALSE;
972 Float_t zvtx = vtx->GetZ();
973 Float_t yvtx = vtx->GetY();
974 Float_t xvtx = vtx->GetX();
975 Float_t deltaZ = zvtx - vtxSPD->GetZ();
976 if(fDebug)Printf("%s:%d deltaz %3.3f ",__FILE__,__LINE__,deltaZ);
982 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
984 Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
987 if(TMath::Abs(deltaZ)>fVtxDeltaZCut)vertexIn = kFALSE;
988 if(fDebug)Printf("%s:%d vertexIn %d ",__FILE__,__LINE__,vertexIn);
994 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
995 if(!esd)return kFALSE;
996 return esd->IsPileupFromSPD();
999 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
1000 if(!esd)return kFALSE;
1001 // add track cuts for which we look for cosmics...
1003 Bool_t isCosmic = kFALSE;
1004 Int_t nTracks = esd->GetNumberOfTracks();
1005 Int_t nCosmicCandidates = 0;
1007 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
1008 AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
1009 if (!track1) continue;
1010 UInt_t status1 = track1->GetStatus();
1011 //If track is ITS stand alone track, skip the track
1012 if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
1013 if(track1->Pt()<fPtMinCosmic) continue;
1014 //Start 2nd track loop to look for correlations
1015 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
1016 AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
1017 if(!track2) continue;
1018 UInt_t status2 = track2->GetStatus();
1019 //If track is ITS stand alone track, skip the track
1020 if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
1021 if(track2->Pt()<fPtMinCosmic) continue;
1022 //Check if back-to-back
1023 Double_t mom1[3],mom2[3];
1024 track1->GetPxPyPz(mom1);
1025 track2->GetPxPyPz(mom2);
1026 TVector3 momv1(mom1[0],mom1[1],mom1[2]);
1027 TVector3 momv2(mom2[0],mom2[1],mom2[2]);
1028 Float_t theta = (float)(momv1.Phi()-momv2.Phi());
1029 if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
1031 Float_t deltaPhi = track1->Phi()-track2->Phi();
1032 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
1034 Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
1035 if(rIsol<fRIsolMinCosmic) continue;
1037 if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
1038 nCosmicCandidates+=1;
1045 fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
1051 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
1054 if(fCollisionType==kPbPb){
1055 if(esd->GetCentrality()){
1056 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
1058 if(cent<0)cent = 101;
1059 if(cent>80||cent<0)return 5;
1060 if(cent>50)return 4;
1061 if(cent>30)return 3;
1062 if(cent>10)return 2;
1069 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
1071 if(fCollisionType==kPbPb){
1072 Float_t cent = aod->GetHeader()->GetCentrality();
1073 if(cent>80||cent<0)return 5;
1074 if(cent>50)return 4;
1075 if(cent>30)return 3;
1076 if(cent>10)return 2;
1084 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
1086 // Terminate analysis
1090 Bool_t AliAnalysisTaskJetServices::CalculateReactionPlaneAngleVZERO(AliAODEvent *aod){
1092 // 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};
1093 if(!aod)return kFALSE;
1094 static bool bFirst = true;
1095 static Double_t v0phi[64] = {0,};
1099 for(int iArm = 0; iArm<2; iArm++){
1100 for(int iRing = 0; iRing<4 ; iRing++){
1101 for(int iSec = 0; iSec<8 ; iSec++){
1102 v0phi[is] = 22.5 + 45. * iSec;
1103 v0phi[is] *= TMath::Pi()/180;
1104 // cout<< is <<" "<< v0phi[is]<<endl;
1113 const AliAODVZERO *aodVZERO = aod->GetVZEROData();
1114 Double_t numYZNA = 0,numXZNA = 0,sumZNA = 0;
1115 Double_t meanXA = 0,meanYA = 0;
1117 Double_t numYZNC = 0,numXZNC = 0,sumZNC = 0;
1118 Double_t meanXC = 0,meanYC = 0;
1122 static Int_t iOldRun = -1;
1123 static Int_t iFoundBin = -1;
1125 if(aod->GetRunNumber()!=iOldRun&&(fp1CalibRPYA)){
1126 // search only or the bin in case of new runs
1128 Int_t ib = fp1CalibRPYA->FindBin(aod->GetRunNumber());
1129 Float_t err = fp1CalibRPYA->GetBinError(ib);
1130 if(err>0){// value can be zero...
1136 while(iFoundBin<0&&(ibLo>0||ibUp<=fp1CalibRPYA->GetNbinsX())){
1137 err = fp1CalibRPYA->GetBinError(ibLo);
1142 err = fp1CalibRPYA->GetBinError(ibUp);
1143 if(err>0)iFoundBin = ibUp;
1149 iOldRun = aod->GetRunNumber();
1152 if(fDebug)Printf("%s:%d iFoundBin %d",(char*)__FILE__,__LINE__,iFoundBin);
1154 if(iFoundBin>0&&(fp1CalibRPYA)){
1155 meanXA = fp1CalibRPXA->GetBinContent(iFoundBin);
1156 meanYA = fp1CalibRPYA->GetBinContent(iFoundBin);
1157 meanXC = fp1CalibRPXC->GetBinContent(iFoundBin);
1158 meanYC = fp1CalibRPYC->GetBinContent(iFoundBin);
1161 if(fDebug)Printf("%s:%d iFoundBin %1.3E %1.3E %1.3E %1.3E",(char*)__FILE__,__LINE__,meanXA,meanYA,meanXC,meanYC);
1163 for (int i=0; i<64; i++) {
1164 Double_t mult = aodVZERO->GetMultiplicity(i);
1165 Double_t phi = v0phi[i];
1167 if (i<32) { //C-side
1168 Double_t wZNC= mult;
1169 numYZNC += sin(2.*phi)*wZNC;
1170 numXZNC += cos(2.*phi)*wZNC;
1173 else if(i>31){ //A-side
1175 numYZNA += sin(2.*phi)*wZNA;
1176 numXZNA += cos(2.*phi)*wZNA;
1190 XC = numXZNC/sumZNC;
1191 YC = numYZNC/sumZNC;
1192 fPsiVZEROC = 0.5*TMath::ATan2(YC-meanYC, XC-meanXC);
1193 if(fPsiVZEROC>TMath::Pi()){fPsiVZEROC-=TMath::Pi();}
1194 if(fPsiVZEROC<0){fPsiVZEROC+=TMath::Pi();}
1197 fh2XYC->Fill(XC-meanXC,YC-meanYC); // control
1198 fp1RPXC->Fill(aod->GetRunNumber(),XC);
1199 fp1RPYC->Fill(aod->GetRunNumber(),YC);
1205 XA = numXZNA/sumZNA;
1206 YA = numYZNA/sumZNA;
1207 fPsiVZEROA = 0.5*TMath::ATan2(YA-meanYA, XA-meanXA);
1208 if(fPsiVZEROA>TMath::Pi()){fPsiVZEROA-=TMath::Pi();}
1209 if(fPsiVZEROA<0){fPsiVZEROA+=TMath::Pi();}
1211 fh2XYA->Fill(XA-meanXA,YA-meanYA); // control
1212 fp1RPXA->Fill(aod->GetRunNumber(),XA);
1213 fp1RPYA->Fill(aod->GetRunNumber(),YA);
1220 Int_t AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
1222 AliAODEvent *aod = 0;
1223 if(fUseAODInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1224 else aod = AODEvent();
1228 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1229 AliAODTrack *tr = aod->GetTrack(it);
1230 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1231 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1232 if(tr->Pt()<fMinTrackPt)continue;
1241 void AliAnalysisTaskJetServices::SetNTrigger(Int_t n){
1244 delete [] fTriggerName;
1245 fTriggerName = new TString [fNTrigger];
1246 delete [] fTriggerBit;fTriggerBit = 0;
1247 fTriggerBit = new UInt_t [fNTrigger];
1254 void AliAnalysisTaskJetServices::SetTrigger(Int_t i,UInt_t it,const char* c){
1259 Printf("%p",&fTriggerName[i]);
1261 fTriggerBit[i] = it;
1262 fTriggerName[i] = c; // placement new
1263 // new(&fTriggerName[i]) TString(c); // placement new
1264 Printf("%s",fTriggerName[i].Data());