1 /*************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: A.Abrahantes, E.Lopez, S.Vallero *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
24 #include <TLorentzVector.h>
26 #include <TObjArray.h>
33 #include "AliAnalyseUE.h"
34 #include "AliAnalysisTaskUE.h"
35 #include "AliAnalysisTask.h"
36 #include "AliHistogramsUE.h"
38 #include "AliAnalysisManager.h"
39 #include "AliAODEvent.h"
40 #include "AliESDEvent.h"
41 #include "AliAODHandler.h"
42 #include "AliAODInputHandler.h"
43 #include "AliAODJet.h"
44 #include "AliAODMCParticle.h"
45 #include "AliAODTrack.h"
46 #include "AliESDtrack.h"
47 #include "AliKFVertex.h"
48 #include "AliMCEvent.h"
49 #include "AliMCEventHandler.h"
52 #include "AliAnalysisHelperJetTasks.h"
53 #include "AliGenPythiaEventHeader.h"
54 #include "AliInputEventHandler.h"
58 ////////////////////////////////////////////////
59 //---------------------------------------------
60 // Class for transverse regions analysis
61 //---------------------------------------------
62 ////////////////////////////////////////////////
67 ClassImp(AliAnalyseUE)
69 //-------------------------------------------------------------------
70 AliAnalyseUE::AliAnalyseUE() :
77 fSimulateChJetPt(kFALSE),
80 fAreaReg(1.5393), // Pi*0.7*0.7
84 fUseChargeHadrons(kFALSE),
85 fUseChPartJet(kFALSE),
86 fUsePositiveCharge(kTRUE),
87 fUseSingleCharge(kFALSE),
90 fJet2DeltaPhiCut(2.616), // 150 degrees
96 fSumPtRegionPosit(0.),
97 fSumPtRegionNegat(0.),
98 fSumPtRegionForward(0.),
99 fSumPtRegionBackward(0.),
100 fMaxPartPtRegion(0.),
101 fNTrackRegionPosit(0),
102 fNTrackRegionNegat(0),
103 fNTrackRegionForward(0),
104 fNTrackRegionBackward(0),
113 //-------------------------------------------------------------------
114 AliAnalyseUE::AliAnalyseUE(const AliAnalyseUE & original) :
116 //fTaskUE(original.fTaskUE),
117 fkAOD(original.fkAOD),
119 fkESD(original.fkESD),
120 fDebug(original.fDebug),
121 fSimulateChJetPt(original.fSimulateChJetPt),
122 fStack(original.fStack),
123 fAnaType(original.fAnaType),
124 fAreaReg(original.fAreaReg),
125 fConeRadius(original.fConeRadius),
126 fFilterBit(original.fFilterBit),
127 fRegionType(original.fRegionType),
128 fUseChargeHadrons(original.fUseChargeHadrons),
129 fUseChPartJet(original.fUseChPartJet),
130 fUsePositiveCharge(original.fUsePositiveCharge),
131 fUseSingleCharge(original.fUseSingleCharge),
132 fOrdering(original.fOrdering),
133 fJet1EtaCut(original.fJet1EtaCut),
134 fJet2DeltaPhiCut(original.fJet2DeltaPhiCut),
135 fJet2RatioPtCut(original.fJet2RatioPtCut),
136 fJet3PtCut(original.fJet3PtCut),
137 fTrackEtaCut(original.fTrackEtaCut),
138 fTrackPtCut(original.fTrackPtCut),
139 fHistos(original.fHistos),
140 fSumPtRegionPosit(original.fSumPtRegionPosit),
141 fSumPtRegionNegat(original.fSumPtRegionNegat),
142 fSumPtRegionForward(original.fSumPtRegionForward),
143 fSumPtRegionBackward(original.fSumPtRegionBackward),
144 fMaxPartPtRegion(original.fMaxPartPtRegion),
145 fNTrackRegionPosit(original.fNTrackRegionPosit),
146 fNTrackRegionNegat(original.fNTrackRegionNegat),
147 fNTrackRegionForward(original.fNTrackRegionForward),
148 fNTrackRegionBackward(original.fNTrackRegionBackward),
149 fSettingsTree(original.fSettingsTree),
150 fLtLabel(original.fLtLabel),
151 fLtMCLabel(original.fLtMCLabel)
156 //-------------------------------------------------------------------
157 AliAnalyseUE & AliAnalyseUE::operator = (const AliAnalyseUE & /*source*/)
159 // assignment operator
164 //-------------------------------------------------------------------
165 AliAnalyseUE::~AliAnalyseUE(){
174 //-------------------------------------------------------------------
175 void AliAnalyseUE::AnalyseMC(TVector3 *jetVect,AliMCEvent *mcEvent, AliGenPythiaEventHeader *pythiaGenHeader,Int_t conePosition, Bool_t useAliStack, Bool_t constrainDistance, Double_t minDistance){
177 // Execute the analysis in case of MC input
178 fSumPtRegionPosit = 0.;
179 fSumPtRegionNegat = 0.;
180 fSumPtRegionForward = 0.;
181 fSumPtRegionBackward = 0.;
182 fMaxPartPtRegion = 0.;
183 fNTrackRegionPosit = 0;
184 fNTrackRegionNegat = 0;
185 fNTrackRegionForward = 0;
186 fNTrackRegionBackward = 0;
188 static Double_t const kPI = TMath::Pi();
189 static Double_t const k270rad = 270.*kPI/180.;
191 //Get Jets from MC header
192 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
193 AliAODJet pythiaGenJets[4];
194 TVector3 jetVectnew[4];
196 for(int ip = 0;ip < nPythiaGenJets;++ip){
199 pythiaGenHeader->TriggerJet(ip,p);
200 TVector3 tempVect(p[0],p[1],p[2]);
201 if ( TMath::Abs(tempVect.Eta())>fJet1EtaCut ) continue;
202 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
203 jetVectnew[iCount].SetXYZ(pythiaGenJets[iCount].Px(), pythiaGenJets[iCount].Py(), pythiaGenJets[iCount].Pz());
207 if (!iCount) return;// no jet in eta acceptance
209 //Search the index of the nearest MC jet to the leading jet reconstructed from the input data
211 if (constrainDistance){
214 for (Int_t i=0; i<iCount; i++){
216 dRTemp = jetVectnew[i].DeltaR(jetVect[0]);
220 deltaR = jetVectnew[i].DeltaR(jetVect[0]);
221 if (deltaR < dRTemp){
227 if (jetVectnew[index].DeltaR(jetVect[0]) > minDistance) return;
230 //Let's add some taste to jet and simulate pt of charged alone
231 //eta and phi are kept as original
232 //Play a Normal Distribution
234 if (fSimulateChJetPt){
236 random = gRandom->Gaus(0.6,0.25);
237 if (random > 0. && random < 1. &&
238 (random * jetVectnew[index].Pt()>6.)) break;
242 //Set new Pt & Fill histogram accordingly
243 Double_t maxPtJet1 = random * jetVectnew[index].Pt();
245 fHistos->FillHistogram("hEleadingPt", maxPtJet1 );
247 if (useAliStack){//Try Stack Information to perform UE analysis
249 AliStack* mcStack = mcEvent->Stack();//Load Stack
250 Int_t nTracksMC = mcStack->GetNtrack();
251 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
253 if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
255 TParticle* mctrk = mcStack->Particle(iTracks);
257 Double_t charge = mctrk->GetPDG()->Charge();
258 Double_t pT = mctrk->Pt();
259 Double_t eta = mctrk->Eta();
260 Int_t pdgCode = mctrk->GetPdgCode();
262 if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
264 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
265 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
266 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
267 fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, maxPtJet1 );
269 fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
271 //We are not interested on stack organization but don't loose track of info
273 TVector3 tempVector = jetVectnew[0];
274 jetVectnew[0] = jetVectnew[index];
275 jetVectnew[index] = tempVector;
277 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );
280 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
281 fSumPtRegionPosit += mctrk->Pt();
282 fNTrackRegionPosit++;
283 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
286 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
287 fSumPtRegionNegat += mctrk->Pt();
288 fNTrackRegionNegat++;
289 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
291 if (region == 2){ //forward
292 fSumPtRegionForward += mctrk->Pt();
293 fNTrackRegionForward++;
294 fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
296 if (region == -2){ //backward
297 fSumPtRegionBackward += mctrk->Pt();
298 fNTrackRegionBackward++;
299 fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
301 } //end loop on stack particles
302 }else{//Try mc Particle
304 TClonesArray* farray = (TClonesArray*)fkAOD->FindListObject("mcparticles");
306 Int_t ntrks = farray->GetEntries();
307 if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
308 for (Int_t i =0 ; i < ntrks; i++){
309 AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
311 if (!(mctrk->IsPhysicalPrimary())) continue;
312 //if (!(mctrk->IsPrimary())) continue;
314 Double_t charge = mctrk->Charge();
315 Double_t pT = mctrk->Pt();
316 Double_t eta = mctrk->Eta();
317 Int_t pdgCode = mctrk->GetPdgCode();
319 if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
321 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
323 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
324 if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi();
325 fHistos->FillHistogram("hdNdEtaPhiDist", deltaPhi, maxPtJet1 );
327 fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
329 //We are not interested on stack organization but don't loose track of info
330 TVector3 tempVector = jetVectnew[0];
331 jetVectnew[0] = jetVectnew[index];
332 jetVectnew[index] = tempVector;
334 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );
336 if (region == 1) { //right
337 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
338 fSumPtRegionPosit += mctrk->Pt();
339 fNTrackRegionPosit++;
340 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
342 if (region == -1) { //left
343 if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
344 fSumPtRegionNegat += mctrk->Pt();
345 fNTrackRegionNegat++;
346 fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
348 if (region == 2){ //forward
349 fSumPtRegionForward += mctrk->Pt();
350 fNTrackRegionForward++;
351 fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
353 if (region == -2){ //backward
354 fSumPtRegionBackward += mctrk->Pt();
355 fNTrackRegionBackward++;
356 fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
359 }//end loop AliAODMCParticle tracks
365 //-------------------------------------------------------------------
366 Bool_t AliAnalyseUE::AnaTypeSelection(TVector3 *jetVect ){
368 // Cut events by jets topology
371 // - Jet1 |eta| < jet1EtaCut
372 // 2 = back to back inclusive
374 // - |Jet1.Phi - Jet2.Phi| > jet2DeltaPhiCut
375 // - Jet2.Pt/Jet1Pt > jet2RatioPtCut
376 // 3 = back to back exclusive
378 // - Jet3.Pt < jet3PtCut
380 Double_t eta=jetVect[0].Eta();
381 if( TMath::Abs(eta) > fJet1EtaCut) {
382 if( fDebug > 1 ) AliInfo("\n Skipping Event...Jet1 |eta| > fJet1EtaCut");
385 // back to back inclusive
386 if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() < 0. ) {
387 if( fDebug > 1 ) AliInfo("\n Skipping Event... no second Jet found");
390 if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() > 0. ) {
391 if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
392 jetVect[1].Pt()/jetVect[0].Pt() < fJet2RatioPtCut ) {
393 if( fDebug > 1 ) AliInfo("\n Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
397 // back to back exclusive
398 if( fAnaType > 2 && fAnaType < 4 && jetVect[2].Pt() > 0. ) {
399 if( jetVect[2].Pt() > fJet3PtCut ) {
400 if( fDebug > 1 ) AliInfo("\n Skipping Event... Jet3.Pt > fJet3PtCut ");
408 //-------------------------------------------------------------------
409 void AliAnalyseUE::FillRegions(Bool_t isNorm2Area, TVector3 *jetVect){
411 // Fill the different topological regions
412 Double_t maxPtJet1 = jetVect[0].Pt();
413 static Double_t const kPI = TMath::Pi();
414 static Double_t const k120rad = 120.*kPI/180.;
415 Double_t const kMyTolerance = 0.0000001;
417 //Area for Normalization
418 // Forward and backward
419 Double_t normArea = 1.;
422 SetRegionArea(jetVect);
423 normArea = 2.*fTrackEtaCut*k120rad ;
424 } else fAreaReg = 1.;
426 Double_t avePosRegion = (fNTrackRegionPosit) ? fSumPtRegionPosit/fNTrackRegionPosit : 0.;
427 Double_t aveNegRegion = (fNTrackRegionNegat) ? fSumPtRegionNegat/fNTrackRegionNegat : 0.;
428 if( avePosRegion > aveNegRegion ) {
429 FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
431 FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
434 //How quantities will be sorted before Fill Min and Max Histogram
435 // 1=Plots will be CDF-like
436 // 2=Plots will be Marchesini-like
437 // 3=Minimum zone is selected as the one having lowest pt per track
438 if( fOrdering == 1 ) {
439 if( fSumPtRegionPosit > fSumPtRegionNegat ) {
440 FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
442 FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
444 if (fNTrackRegionPosit > fNTrackRegionNegat ) {
445 FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
447 FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
449 } else if( fOrdering == 2 ) {
450 if (fSumPtRegionPosit > fSumPtRegionNegat) {
451 FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
452 FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
454 FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
455 FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
457 } else if( fOrdering == 3 ){
458 if (avePosRegion > aveNegRegion) {
459 FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
460 FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
462 FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
463 FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
466 fHistos->FillHistogram("hRegionMaxPartPtMaxVsEt",maxPtJet1, fMaxPartPtRegion);
468 // Compute pedestal like magnitudes
469 fHistos->FillHistogram("hRegionDiffSumPtVsEt",maxPtJet1, (TMath::Abs(fSumPtRegionPosit-fSumPtRegionNegat)/(2.0*fAreaReg))+kMyTolerance);
470 fHistos->FillHistogram("hRegionAveSumPtVsEt", maxPtJet1, (fSumPtRegionPosit+fSumPtRegionNegat)/(2.0*fAreaReg));
472 // Transverse as a whole
473 fHistos->FillHistogram("hRegTransMult", maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat, (fNTrackRegionPosit + fNTrackRegionNegat)/(2.0*fAreaReg));
474 fHistos->FillHistogram("hRegTransSumPtVsMult",maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat , (fSumPtRegionNegat + fSumPtRegionPosit)/(2.0 *fAreaReg));
476 // Fill Histograms for Forward and away side w.r.t. leading jet direction
478 //fHistos->FillHistogram("hRegForwardSumPtVsEt",maxPtJet1, fSumPtRegionForward/normArea );
479 //fHistos->FillHistogram("hRegForwardMultVsEt",maxPtJet1, fNTrackRegionForward/normArea );
480 //fHistos->FillHistogram("hRegBackwardSumPtVsEt",maxPtJet1, fSumPtRegionBackward/normArea );
481 //fHistos->FillHistogram("hRegBackwardMultVsEt",maxPtJet1, fNTrackRegionBackward/normArea);
483 // Multiplicity dependence
484 fHistos->FillHistogram("hRegForwardMult", maxPtJet1, fNTrackRegionForward, fNTrackRegionForward/normArea);
485 fHistos->FillHistogram("hRegForwardSumPtvsMult", maxPtJet1, fNTrackRegionForward,fSumPtRegionForward/normArea);
486 fHistos->FillHistogram("hRegBackwardMult", maxPtJet1, fNTrackRegionBackward, fNTrackRegionBackward/normArea );
487 fHistos->FillHistogram("hRegBackwardSumPtvsMult", maxPtJet1, fNTrackRegionBackward,fSumPtRegionBackward/normArea);
491 //-------------------------------------------------------------------
492 void AliAnalyseUE::FindMaxMinRegions(TVector3 *jetVect, Int_t conePosition, Int_t mctrue=0, Int_t eff=0){
494 // If mctrue = 1 consider branch AliAODMCParticles
495 // If eff = 1 track cuts for efficiency studies
497 // Identify the different topological zones
498 fSumPtRegionPosit = 0.;
499 fSumPtRegionNegat = 0.;
500 fSumPtRegionForward = 0.;
501 fSumPtRegionBackward = 0.;
502 fMaxPartPtRegion = 0.;
503 fNTrackRegionPosit = 0;
504 fNTrackRegionNegat = 0;
505 fNTrackRegionForward = 0;
506 fNTrackRegionBackward = 0;
507 static Double_t const kPI = TMath::Pi();
508 static Double_t const kTWOPI = 2.*kPI;
509 static Double_t const k270rad = 270.*kPI/180.;
510 Double_t const kMyTolerance = 0.0000001;
513 TClonesArray *tca = 0;
515 nTracks = fkAOD->GetNTracks();
516 if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
518 tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
520 Printf("%s:%d No AODMC Branch found !!!",(char*)__FILE__,__LINE__);
523 nTracks = tca->GetEntriesFast();
524 if (fDebug > 1) AliInfo(Form(" ==== AOD MC particles = %d \n ",nTracks));
527 //If UE task d0 distribution is not filled
529 if (fHistos->GetHistograms()->FindObject("hDCAxy")) flag_tmp = 1;
531 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
533 AliAODTrack *partRECO=0;
534 AliAODMCParticle *partMC=0;
535 Double_t charge=-999.;
541 // get track reco or true
543 partRECO = dynamic_cast<AliAODTrack*>(fkAOD->GetTrack(ipart));
547 partMC = dynamic_cast<AliAODMCParticle*>(tca->At(ipart));
550 charge = partMC->Charge();
553 pdgcode = partMC->GetPdgCode();
559 if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
562 if (!mctrue && !eff ){
563 if( !TrackSelected(partRECO)) continue; //track selection for data and MC reconstructed
565 if (fkESD && fkESD->GetNumberOfTracks() ){
566 AliInfo("READING ESD *************************************************");
567 Int_t id = partRECO->GetID();
568 AliESDtrack *trackESD;
569 trackESD = (AliESDtrack*)fkESD->GetTrack(id);
572 trackESD->GetImpactParameters(d0,z);
573 fHistos->FillHistogram("hDCAxy", d0, jetVect[0].Pt());
574 }else AliInfo("NO TRACKS ************************************************")
579 if (!TrackSelectedEfficiency(partRECO )) continue; //track selection for MC reconstructed for efficiency studies
580 if(fkESD && fkESD->GetNumberOfTracks()){
581 Int_t id = partRECO->GetID();
582 AliESDtrack * trackESD = (AliESDtrack*) fkESD->GetTrack(id);
585 trackESD->GetImpactParameters(d0,z);
586 fHistos->FillHistogram("hDCAxyPrimary", d0, jetVect[0].Pt());
591 if ( !(TrackMCSelected(charge, pt, eta, pdgcode) && partMC->IsPhysicalPrimary())) continue; //track selection for MC true
594 TVector3 partVect(part->Px(), part->Py(), part->Pz());
595 Bool_t isFlagPart = kTRUE;
596 Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
597 if( deltaPhi > kTWOPI ) deltaPhi-= kTWOPI;
598 if (fAnaType != 4 ) fHistos->FillHistogram(Form("hdNdEtaPhiDist%s",suffix.Data()),deltaPhi, jetVect[0].Pt());
599 else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
600 fHistos->FillHistogram(Form("hdNdEtaPhiDist%s",suffix.Data()),deltaPhi, jetVect[0].Pt());
604 fHistos->FillHistogram(Form("hFullRegPartPtDistVsEt%s",suffix.Data()), part->Pt(), jetVect[0].Pt());
606 Int_t region = IsTrackInsideRegion( jetVect, &partVect, conePosition );
608 if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
609 fSumPtRegionPosit += part->Pt();
610 fNTrackRegionPosit++;
611 fHistos->FillHistogram(Form("hTransRegPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
614 if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
615 fSumPtRegionNegat += part->Pt();
616 fNTrackRegionNegat++;
617 fHistos->FillHistogram(Form("hTransRegPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
619 if (region == 2){ //forward
620 fSumPtRegionForward += part->Pt();
621 fNTrackRegionForward++;
622 fHistos->FillHistogram(Form("hRegForwardPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
624 if (region == -2){ //backward
625 fSumPtRegionBackward += part->Pt();
626 fNTrackRegionBackward++;
627 fHistos->FillHistogram(Form("hRegBackwardPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
629 }//end loop AOD tracks
633 //-------------------------------------------------------------------
634 /*TVector3 AliAnalyseUE::GetLeadingTracksMC(AliMCEvent *mcEvent){
638 Double_t maxPtJet1 = 0.;
643 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
644 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
645 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
649 TObjArray* tracks = SortChargedParticlesMC();
651 nJets = tracks->GetEntriesFast();
654 TParticle* jet = (TParticle*)tracks->At(0);
655 maxPtJet1 = jet->Pt();
656 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
666 //-------------------------------------------------------------------
667 TVector3 AliAnalyseUE::GetLeadingTracksMC(AliMCEvent *mcEvent){
671 Double_t maxPtJet1 = 0.;
676 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
677 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
678 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
682 TObjArray* tracks = SortChargedParticlesMC();
684 nJets = tracks->GetEntriesFast();
687 AliAODMCParticle* jet = (AliAODMCParticle*)tracks->At(0);
688 fLtMCLabel = jet->GetLabel();
689 maxPtJet1 = jet->Pt();
690 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
699 //-------------------------------------------------------------------
700 TVector3 AliAnalyseUE::GetOrderedClusters(TString aodBranch, Bool_t chargedJets, Double_t chJetPtMin){
702 // jets from AOD, on-the-fly or leading particle
703 Double_t maxPtJet1 = 0.;
705 Double_t maxPtJet2 = 0.; // jet 2 need for back to back inclusive
707 Double_t maxPtJet3 = 0.; // jet 3 need for back to back exclusive
711 jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
712 jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
713 jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
716 //TClonesArray* fArrayJets;
717 TObjArray* arrayJets;
718 // 1) JETS FROM AOD BRANCH (standard, non-standard or delta)
719 if (!chargedJets && fAnaType != 4 ) {
720 AliInfo(" ==== Read AODs !");
721 AliInfo(Form(" ==== Reading Branch: %s ", aodBranch.Data()));
722 arrayJets = (TObjArray*)fkAOD->GetList()->FindObject(aodBranch.Data());
724 AliFatal(" No jet-array! ");
728 // Find Leading Jets 1,2,3
729 // (could be skipped if Jets are sort by Pt...)
730 nJets=arrayJets->GetEntries();
731 for( Int_t i=0; i<nJets; ++i ) {
732 AliAODJet* jet = (AliAODJet*)arrayJets->At(i);
733 Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
735 if( jetPt > maxPtJet1 ) {
736 maxPtJet3 = maxPtJet2; index3 = index2;
737 maxPtJet2 = maxPtJet1; index2 = index1;
738 maxPtJet1 = jetPt; index1 = i;
739 } else if( jetPt > maxPtJet2 ) {
740 maxPtJet3 = maxPtJet2; index3 = index2;
741 maxPtJet2 = jetPt; index2 = i;
742 } else if( jetPt > maxPtJet3 ) {
743 maxPtJet3 = jetPt; index3 = i;
748 AliAODJet *jet =(AliAODJet*) arrayJets->At(index1);
749 if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
752 AliAODJet* jet= (AliAODJet*) arrayJets->At(index2);
753 if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
756 AliAODJet* jet = (AliAODJet*) arrayJets->At(index3);
757 if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
763 // 2) ON-THE-FLY CDF ALGORITHM
765 arrayJets = FindChargedParticleJets(chJetPtMin);
767 nJets = arrayJets->GetEntriesFast();
770 AliAODJet* jet = (AliAODJet*)arrayJets->At(0);
771 maxPtJet1 = jet->Pt();
772 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
776 AliAODJet* jet = (AliAODJet*)arrayJets->At(1);
777 maxPtJet2 = jet->Pt();
778 jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
782 AliAODJet* jet = (AliAODJet*)arrayJets->At(2);
783 maxPtJet3 = jet->Pt();
784 jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
793 // 3) LEADING PARTICLE
795 TObjArray* tracks = SortChargedParticles();
797 nJets = tracks->GetEntriesFast();
800 AliAODTrack* jet = (AliAODTrack*)tracks->At(0);
801 fLtLabel = jet->GetLabel();
802 maxPtJet1 = jet->Pt();
803 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
810 fHistos->FillHistogram("hNJets",nJets);
817 //-------------------------------------------------------------------
818 void AliAnalyseUE::Initialize(AliAnalysisTaskUE& taskUE){
819 /*void AliAnalyseUE::Initialize(AliAnalysisTask& task){// when correction task is in trunk
820 if (task->InheritsFrom("AliAnalysisTaskUE")){
821 AliAnalysisTaskUE *taskUE = dynamic_cast<AliAnalysisTaskUE*> task;
823 else if (task->InheritsFrom("AliAnalysisTaskCorrectionsUE")){
824 AliAnalysisTaskCorrectionsUE *taskUE = dynamic_cast<AliAnalysisTaskCorrectionsUE*> task;
828 //Get principal settings from current instance of UE analysis-task
829 fAnaType = taskUE.GetAnaTopology();
830 fkAOD = taskUE.GetAOD();
831 fConeRadius = taskUE.GetConeRadius();
832 fDebug = taskUE.GetDebugLevel();
833 fFilterBit = taskUE.GetFilterBit();
834 fJet1EtaCut = taskUE.GetJet1EtaCut();
835 fJet2DeltaPhiCut = taskUE.GetJet2DeltaPhiCut();
836 fJet2RatioPtCut = taskUE.GetJet2RatioPtCut();
837 fJet3PtCut = taskUE.GetJet3PtCut();
838 fOrdering = taskUE.GetPtSumOrdering() ;
839 fRegionType = taskUE.GetRegionType();
840 fSimulateChJetPt = taskUE.GetSimulateChJetPt();
841 fTrackEtaCut = taskUE.GetTrackEtaCut();
842 fTrackPtCut = taskUE.GetTrackPtCut();
843 fHistos = taskUE.fHistosUE;
844 fUseChargeHadrons = taskUE.GetUseChargeHadrons();
845 fUseChPartJet = taskUE.GetUseChPartJet();
846 fUsePositiveCharge = taskUE.GetUseNegativeChargeType();
847 fUseSingleCharge = taskUE.GetUseSingleCharge();
851 //-------------------------------------------------------------------
852 void AliAnalyseUE::Initialize(Int_t anaType, AliAODEvent* aod,Double_t coneRadius, Int_t debug, Int_t filterBit, Double_t jet1EtaCut, Double_t jet2DeltaPhiCut, Double_t jet2RatioPtCut, Double_t jet3PtCut, Int_t ordering, Int_t regionType,Bool_t simulateChJetPt, Double_t trackEtaCut, Double_t trackPtCut, Bool_t useChargeHadrons, Bool_t useChPartJet, Bool_t usePositiveCharge, Bool_t useSingleCharge, AliHistogramsUE* histos){
856 fConeRadius = coneRadius;
858 fFilterBit = filterBit;
859 fJet1EtaCut = jet1EtaCut;
860 fJet2DeltaPhiCut = jet2DeltaPhiCut;
861 fJet2RatioPtCut = jet2RatioPtCut;
862 fJet3PtCut = jet3PtCut;
863 fOrdering = ordering;
864 fRegionType = regionType;
865 fSimulateChJetPt = simulateChJetPt;
866 fTrackEtaCut = trackEtaCut;
867 fTrackPtCut = trackPtCut;
868 fUseChargeHadrons = useChargeHadrons;
869 fUseChPartJet = useChPartJet;
870 fUsePositiveCharge = usePositiveCharge;
871 fUseSingleCharge = useSingleCharge;
876 //-------------------------------------------------------------------
877 Bool_t AliAnalyseUE::TriggerSelection(AliInputEventHandler* inputHandler){
879 //Use AliPhysicsSelection to select good events
880 if( !inputHandler->InheritsFrom("AliAODInputHandler") ) { // input AOD
881 if (inputHandler->IsEventSelected()) {
882 if (fDebug > 1) AliInfo(" Trigger Selection: event ACCEPTED ... ");
884 if (fDebug > 1) AliInfo(" Trigger Selection: event REJECTED ... ");
894 //-------------------------------------------------------------------
895 Bool_t AliAnalyseUE::VertexSelection(AliAODEvent *aod, Int_t tracks, Double_t zed ){
897 //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
898 Int_t nVertex = aod->GetNumberOfVertices();
899 if( nVertex > 0 ) { // Only one vertex (reject pileup)
900 AliAODVertex* vertex = (AliAODVertex*)aod->GetPrimaryVertex();
901 Int_t nTracksPrim = vertex->GetNContributors();
902 Double_t zVertex = vertex->GetZ();
903 if (fDebug > 1) AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
904 // Select a quality vertex by number of tracks?
905 if( nTracksPrim < tracks || TMath::Abs(zVertex) > zed ) {
906 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
909 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
911 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
918 //-------------------------------------------------------------------
919 Bool_t AliAnalyseUE::VertexSelectionOld(AliAODEvent *aod ){
921 AliKFVertex primVtx(*(aod->GetPrimaryVertex()));
922 Int_t nTracksPrim=primVtx.GetNContributors();
923 if (fDebug > 1) AliInfo(Form(" Primary-vertex Selection: %d",nTracksPrim));
925 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
928 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
933 // PRIVATE METHODS **************************************************
935 TObjArray* AliAnalyseUE::FindChargedParticleJets( Double_t chJetPtMin )
937 // Return a TObjArray of "charged particle jets"
939 // Charged particle jet definition from reference:
940 // "Charged jet evolution and the underlying event
941 // in proton-antiproton collisions at 1.8 TeV"
942 // PHYSICAL REVIEW D 65 092002, CDF Collaboration
944 // We defined "jets" as circular regions in eta-phi space with
945 // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
946 // Our jet algorithm is as follows:
947 // 1- Order all charged particles according to their pT .
948 // 2- Start with the highest pT particle and include in the jet all
949 // particles within the radius R=0.7 considering each particle
950 // in the order of decreasing pT and recalculating the centroid
951 // of the jet after each new particle is added to the jet .
952 // 3- Go to the next highest pT particle not already included in
953 // a jet and add to the jet all particles not already included in
954 // a jet within R=0.7.
955 // 4- Continue until all particles are in a jet.
956 // We defined the transverse momentum of the jet to be
957 // the scalar pT sum of all the particles within the jet, where pT
958 // is measured with respect to the beam axis
960 // 1 - Order all charged particles according to their pT .
961 Int_t nTracks = fkAOD->GetNTracks();
962 if( !nTracks ) return 0;
963 TObjArray tracks(nTracks);
965 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
966 AliAODTrack* part = fkAOD->GetTrack( ipart );
967 if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
968 if( !part->Charge() ) continue;
969 if( part->Pt() < fTrackPtCut ) continue;
970 tracks.AddLast(part);
972 QSortTracks( tracks, 0, tracks.GetEntriesFast() );
974 nTracks = tracks.GetEntriesFast();
975 if( !nTracks ) return 0;
977 TObjArray *jets = new TObjArray(nTracks);
978 TIter itrack(&tracks);
980 // 2- Start with the highest pT particle ...
982 AliAODTrack* track = (AliAODTrack*)itrack.Next();
983 if( !track ) continue;
987 pt = track->Pt(); // Use the energy member to store Pt
988 jets->AddLast( new TLorentzVector(px, py, pz, pt) );
989 tracks.Remove( track );
990 TLorentzVector* jet = (TLorentzVector*)jets->Last();
991 jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
992 // 3- Go to the next highest pT particle not already included...
994 Double_t fPt = jet->E();
995 while ( (track1 = (AliAODTrack*)(itrack.Next())) ) {
996 Double_t tphi = track1->Phi(); // here Phi is from 0 <-> 2Pi
997 if (tphi > TMath::Pi()) tphi -= 2. * TMath::Pi(); // convert to -Pi <-> Pi
998 Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-tphi);
999 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +dphi*dphi );
1000 if( r < fConeRadius ) {
1001 fPt = jet->E()+track1->Pt(); // Scalar sum of Pt
1002 // recalculating the centroid
1003 Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
1004 Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
1005 jet->SetPtEtaPhiE( 1., eta, phi, fPt );
1006 tracks.Remove( track1 );
1011 nTracks = tracks.GetEntries();
1012 // 4- Continue until all particles are in a jet.
1014 } // end while nTracks
1016 // Convert to AODjets....
1017 Int_t njets = jets->GetEntriesFast();
1018 TObjArray* aodjets = new TObjArray(njets);
1019 aodjets->SetOwner(kTRUE);
1020 for(Int_t ijet=0; ijet<njets; ++ijet) {
1021 TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
1022 if (jet->E() < chJetPtMin) continue;
1023 Float_t px, py,pz,en; // convert to 4-vector
1024 px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi)
1025 py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi)
1026 pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
1027 en = TMath::Sqrt(px * px + py * py + pz * pz);
1029 aodjets->AddLast( new AliAODJet(px, py, pz, en) );
1034 // Order jets according to their pT .
1035 QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
1038 if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
1044 //____________________________________________________________________
1045 void AliAnalyseUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
1048 // Fill average particle Pt of control regions
1051 fHistos->FillHistogram("hRegionAvePartPtMaxVsEt", leadingE, ptMax);
1053 fHistos->FillHistogram("hRegionAvePartPtMinVsEt", leadingE, ptMin);
1054 // MAke distributions for UE comparison with MB data
1055 fHistos->FillHistogram("hMinRegAvePt", ptMin);
1059 //____________________________________________________________________
1060 void AliAnalyseUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin )
1063 // Fill Nch multiplicity of control regions
1066 fHistos->FillHistogram("hRegionMultMaxVsEt", leadingE, nTrackPtmax);
1067 fHistos->FillHistogram("hRegionMultMax", nTrackPtmax);
1070 fHistos->FillHistogram("hRegionMultMinVsEt", leadingE, nTrackPtmin );
1071 fHistos->FillHistogram("hRegionMultMin", nTrackPtmin);
1073 // MAke distributions for UE comparison with MB data
1074 fHistos->FillHistogram("hMinRegSumPtvsMult", nTrackPtmin,ptMin);
1078 //____________________________________________________________________
1079 void AliAnalyseUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin )
1081 // Fill sumPt of control regions
1084 fHistos->FillHistogram("hRegionSumPtMaxVsEt", leadingE, ptMax);
1087 fHistos->FillHistogram("hRegionSumPtMinVsEt", leadingE, ptMin);
1089 // MAke distributions for UE comparison with MB data
1090 fHistos->FillHistogram("hMinRegSumPt", ptMin);
1091 fHistos->FillHistogram("hMinRegSumPtJetPtBin", leadingE, ptMin);
1092 fHistos->FillHistogram("hMaxRegSumPtJetPtBin", leadingE, ptMax);
1096 //-------------------------------------------------------------------
1097 Int_t AliAnalyseUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect, Int_t conePosition)
1099 // return de region in delta phi
1100 // -1 negative delta phi
1101 // 1 positive delta phi
1103 static const Double_t k60rad = 60.*TMath::Pi()/180.;
1104 static const Double_t k120rad = 120.*TMath::Pi()/180.;
1107 if( fRegionType == 1 ) {
1108 if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
1109 // transverse regions
1110 if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; //left
1111 if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1; //right
1112 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) < k60rad ) region = 2; //forward
1113 if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) > k120rad ) region = -2; //backward
1115 } else if( fRegionType == 2 ) {
1117 Double_t deltaR = 0.;
1119 TVector3 positVect,negatVect;
1120 if (conePosition==1){
1121 positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
1122 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
1123 }else if (conePosition==2){
1124 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
1125 positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2());
1126 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2());
1127 }else if (conePosition==3){
1128 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
1129 Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) +
1130 jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt());
1131 //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) +
1132 // jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag());
1133 positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2());
1134 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2());
1136 if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
1138 deltaR = positVect.DrEtaPhi(*partVect);
1139 } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
1141 deltaR = negatVect.DrEtaPhi(*partVect);
1144 if (deltaR > fConeRadius) region = 0;
1147 AliError("Unknow region type");
1154 //-------------------------------------------------------------------
1155 void AliAnalyseUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
1157 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1159 static TObject *tmp;
1160 static int i; // "static" to save stack space
1163 while (last - first > 1) {
1167 while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
1169 while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
1185 if (j - first < last - (j + 1)) {
1186 QSortTracks(a, first, j);
1187 first = j + 1; // QSortTracks(j + 1, last);
1189 QSortTracks(a, j + 1, last);
1190 last = j; // QSortTracks(first, j);
1196 //-------------------------------------------------------------------
1197 void AliAnalyseUE::QSortTracksMC(TObjArray &a, Int_t first, Int_t last)
1199 // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1201 static TObject *tmp;
1202 static int i; // "static" to save stack space
1205 while (last - first > 1) {
1209 while (++i < last && ((AliAODMCParticle*)a[i])->Pt() > ((AliAODMCParticle*)a[first])->Pt() )
1211 while (--j > first && ((AliAODMCParticle*)a[j])->Pt() < ((AliAODMCParticle*)a[first])->Pt() )
1227 if (j - first < last - (j + 1)) {
1228 QSortTracksMC(a, first, j);
1229 first = j + 1; // QSortTracks(j + 1, last);
1231 QSortTracksMC(a, j + 1, last);
1232 last = j; // QSortTracks(first, j);
1241 //-------------------------------------------------------------------
1242 void AliAnalyseUE::SetRegionArea(TVector3 *jetVect)
1245 Double_t areaCorrFactor=0.;
1246 Double_t deltaEta = 0.;
1247 if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
1248 else if (fRegionType==2){
1249 deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
1250 if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
1252 areaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
1253 (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
1254 fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-areaCorrFactor;
1256 }else AliWarning("Unknown Region Type");
1257 if (fDebug>10) AliInfo(Form("\n dEta=%5.3f Angle =%5.3f Region Area = %5.3f Corr Factor=%5.4f \n",deltaEta,TMath::ACos(deltaEta/fConeRadius),fAreaReg,areaCorrFactor));
1261 //____________________________________________________________________
1262 TObjArray* AliAnalyseUE::SortChargedParticles()
1264 // return an array with all charged particles ordered according to their pT .
1265 Int_t nTracks = fkAOD->GetNTracks();
1266 if( !nTracks ) return 0;
1267 TObjArray* tracks = new TObjArray(nTracks);
1269 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1270 AliAODTrack* part = fkAOD->GetTrack( ipart );
1271 if( !part->TestFilterBit( fFilterBit ) ) continue; // track cut selection
1272 if( !part->Charge() ) continue;
1273 if( part->Pt() < fTrackPtCut ) continue;
1274 tracks->AddLast( part );
1276 QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
1278 nTracks = tracks->GetEntriesFast();
1279 if( !nTracks ) return 0;
1284 //____________________________________________________________________
1285 /*TObjArray* AliAnalyseUE::SortChargedParticlesMC()
1287 // return an array with all charged particles ordered according to their pT .
1288 AliStack *mcStack = fkMC->Stack();
1289 Int_t nTracks = mcStack->GetNtrack();
1290 if( !nTracks ) return 0;
1291 TObjArray* tracks = new TObjArray(nTracks);
1293 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1294 if (!(mcStack->IsPhysicalPrimary(ipart))) continue;
1295 TParticle *part = mcStack->Particle(ipart);
1297 if( !part->GetPDG()->Charge() ) continue;
1298 if( part->Pt() < fTrackPtCut ) continue;
1299 tracks->AddLast( part );
1301 QSortTracksMC( *tracks, 0, tracks->GetEntriesFast() );
1303 nTracks = tracks->GetEntriesFast();
1304 if( !nTracks ) return 0;
1310 //____________________________________________________________________
1311 TObjArray* AliAnalyseUE::SortChargedParticlesMC()
1313 // return an array with all charged particles ordered according to their pT .
1314 TClonesArray *tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1316 Printf("No branch!!!");
1319 Int_t nTracks = tca->GetEntriesFast();
1320 if( !nTracks ) return 0;
1321 TObjArray* tracks = new TObjArray(nTracks);
1323 for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1324 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(ipart));
1326 if(!part->IsPhysicalPrimary())continue;
1327 if( part->Pt() < fTrackPtCut ) continue;
1328 if(!part->Charge()) continue;
1329 tracks->AddLast( part );
1331 QSortTracksMC( *tracks, 0, tracks->GetEntriesFast() );
1333 nTracks = tracks->GetEntriesFast();
1334 if( !nTracks ) return 0;
1340 //____________________________________________________________________
1341 Bool_t AliAnalyseUE::TrackMCSelected(Double_t charge, Double_t pT, Double_t eta, Int_t pdgCode)const{
1343 // MC track selection
1344 Double_t const kMyTolerance = 0.0000001;
1345 //if (charge == 0. || charge == -99.) return kFALSE;
1346 if (charge < kMyTolerance || charge + 99 < kMyTolerance) return kFALSE;
1348 if ( fUseSingleCharge ) { // Charge selection
1349 if ( fUsePositiveCharge && charge < 0.) return kFALSE; // keep Positives
1350 if ( !fUsePositiveCharge && charge > 0.) return kFALSE; // keep Negatives
1353 //Kinematics cuts on particle
1354 if ((pT < fTrackPtCut) || (TMath::Abs(eta) > fTrackEtaCut )) return kFALSE;
1356 Bool_t isHadron = TMath::Abs(pdgCode)==211 ||
1357 TMath::Abs(pdgCode)==2212 ||
1358 TMath::Abs(pdgCode)==321;
1360 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1366 //____________________________________________________________________
1367 Bool_t AliAnalyseUE::TrackSelected(AliAODTrack* part)const{
1369 // Real track selection
1370 if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1371 if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1372 // PID Selection: Reject everything but hadrons
1373 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
1374 part->GetMostProbablePID()==AliAODTrack::kKaon ||
1375 part->GetMostProbablePID()==AliAODTrack::kProton;
1376 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1378 if ( !part->Charge() ) return kFALSE; //Only charged
1379 if ( fUseSingleCharge ) { // Charge selection
1380 if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1381 if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1384 if ( part->Pt() < fTrackPtCut ) return kFALSE;
1385 if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;
1390 //____________________________________________________________________
1391 Bool_t AliAnalyseUE::TrackSelectedEfficiency(AliAODTrack* part)const{
1393 if (!fStack) AliWarning("Attention: stack not set from mother task");
1394 // Real track selection
1395 if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1396 if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1397 // PID Selection: Reject everything but hadrons
1398 Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
1399 part->GetMostProbablePID()==AliAODTrack::kKaon ||
1400 part->GetMostProbablePID()==AliAODTrack::kProton;
1401 if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1403 if ( !part->Charge() ) return kFALSE; //Only charged
1404 if ( fUseSingleCharge ) { // Charge selection
1405 if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1406 if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1409 /*if ( part->Pt() < fTrackPtCut ) return kFALSE;
1410 if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;*/
1412 //Check if there was a primary fulfilling the required cuts
1413 Double_t charge=-999.;
1418 Int_t label = part->GetLabel();
1419 TClonesArray *tca=0;
1420 tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1421 if(!tca)return kFALSE;
1422 //TParticle *partMC = fStack->ParticleFromTreeK(label);
1423 AliAODMCParticle *partMC=dynamic_cast<AliAODMCParticle*>(tca->At(TMath::Abs(label)));
1424 if(!partMC)return kFALSE;
1425 if (!partMC->IsPhysicalPrimary())return kFALSE;
1426 //charge = ((TParticlePDG*)partMC->GetPDG())->Charge();
1427 charge = partMC->Charge();
1429 eta = partMC->Eta();
1430 pdgcode = partMC->GetPdgCode();
1432 if (!TrackMCSelected(charge, pt, eta, pdgcode)) return kFALSE;