]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/correlationHF/AliAnalysisTaskSEDplusCorrelations.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskSEDplusCorrelations.cxx
CommitLineData
cb7c2594 1
35151011 2/**************************************************************************
3 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
17/* $Id$ */
18
19//*************************************************************************
20// Class AliAnalysisTaskSEDplusCorrelations
21// AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and
22//comparison of heavy-flavour decay candidates
23// to MC truth (kinematics stored in the AOD)
24// Authors: Sadhana Dash (correlation)
25/////////////////////////////////////////////////////////////
26
27#include <TClonesArray.h>
28#include <TNtuple.h>
29#include <TCanvas.h>
30#include <TList.h>
31#include <TString.h>
32#include <TDatabasePDG.h>
33
34#include <AliAnalysisDataSlot.h>
35#include <AliAnalysisDataContainer.h>
36#include "AliAnalysisManager.h"
37#include "AliRDHFCutsDplustoKpipi.h"
38#include "AliAODHandler.h"
39#include "AliAODEvent.h"
40#include "AliAODPidHF.h"
41#include "AliAODVertex.h"
42#include "AliAODTrack.h"
43#include "AliAODRecoDecayHF3Prong.h"
44#include "AliAnalysisVertexingHF.h"
45#include "AliAnalysisTaskSE.h"
cb7c2594 46#include "AliVertexingHFUtils.h"
35151011 47#include "AliAnalysisTaskSEDplusCorrelations.h"
48#include "AliNormalizationCounter.h"
49#include "AliVParticle.h"
50#include "AliHFAssociatedTrackCuts.h"
51#include "AliReducedParticle.h"
52#include "AliHFCorrelator.h"
53
54
55using std::cout;
56using std::endl;
57
58ClassImp(AliAnalysisTaskSEDplusCorrelations)
59
60
61//________________________________________________________________________
62AliAnalysisTaskSEDplusCorrelations::AliAnalysisTaskSEDplusCorrelations():
63AliAnalysisTaskSE(),
64 fOutput(0),
65 fCorrelator(0),
66 fSelect(0),
67 fDisplacement(0),
68 fHistNEvents(0),
cb7c2594 69 fTrig(kTRUE),
35151011 70 fEventMix(0),
35151011 71 fUpmasslimit(1.965),
72 fLowmasslimit(1.765),
73 fNPtBins(0),
74 fBinWidth(0.002),
75 fListCuts(0),
cb7c2594 76//fListCutsAsso(0),
35151011 77 fRDCutsAnalysis(0),
78 fCuts(0),
79 fCounter(0),
80 fReadMC(kFALSE),
81 fUseBit(kTRUE),
82 fMixing(kFALSE),
cb7c2594 83 fSystem(kFALSE),
84 fReco(kTRUE)
35151011 85{
86 // Default constructor
cb7c2594 87
baf235c6 88 for(Int_t i=0;i<3*kMaxPtBins;i++){
cb7c2594 89
5d709b96 90 fMassHistK0S[i]=0;
91 fLeadPt[i]=0;
92 fPtSig[i]=0;
93 fMassHist[i]=0;
94 fMassHistTC[i]=0;
95 fMassHistOrigC[i]=0;
96 fMassHistOrigB[i]=0;
97 fMassHistMC[i]=0;
baf235c6 98 }
cb7c2594 99
35151011 100 for(Int_t i=0;i<kMaxPtBins+1;i++){
101 fArrayBinLimits[i]=0;
102 }
103
104}
105
106//________________________________________________________________________
107AliAnalysisTaskSEDplusCorrelations::AliAnalysisTaskSEDplusCorrelations(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana, AliHFAssociatedTrackCuts *asstrkcuts):
108 AliAnalysisTaskSE(name),
109 fOutput(0),
110 fCorrelator(0),
111 fSelect(0),
112 fDisplacement(0),
113 fHistNEvents(0),
cb7c2594 114 fTrig(kTRUE),
35151011 115 fEventMix(0),
35151011 116 fUpmasslimit(1.965),
117 fLowmasslimit(1.765),
118 fNPtBins(0),
119 fBinWidth(0.002),
120 fListCuts(0),
cb7c2594 121 //fListCutsAsso(0),
35151011 122 fRDCutsAnalysis(dpluscutsana),
123 fCuts(asstrkcuts),
124 fCounter(0),
125 fReadMC(kFALSE),
126 fUseBit(kTRUE),
127 fMixing(kFALSE),
cb7c2594 128 fSystem(kFALSE),
129 fReco(kTRUE)
35151011 130{
131 //
132 // Standrd constructor
133 //
134 fNPtBins=fRDCutsAnalysis->GetNPtBins();
135
cb7c2594 136
baf235c6 137 for(Int_t i=0;i<3*kMaxPtBins;i++){
baf235c6 138
5d709b96 139 fMassHistK0S[i]=0;
140 fLeadPt[i]=0;
141 fPtSig[i]=0;
142 fMassHist[i]=0;
143 fMassHistTC[i]=0;
144 fMassHistOrigC[i]=0;
145 fMassHistOrigB[i]=0;
146 fMassHistMC[i]=0;
cb7c2594 147
148
149
150 }
35151011 151
152 for(Int_t i=0;i<kMaxPtBins+1;i++){
153 fArrayBinLimits[i]=0;
154 }
155
156
157 // Default constructor
158 // Output slot #1 writes into a TList container
159 DefineOutput(1,TList::Class()); //My private output
160 // Output slot #2 writes cut to private output
35151011 161 DefineOutput(2,TList::Class());
162 // Output slot #3 writes cut to private output
163 DefineOutput(3,AliNormalizationCounter::Class());
cb7c2594 164 //DefineOutput(4,AliHFAssociatedTrackCuts::Class());
35151011 165
166
167}
168
169//________________________________________________________________________
170AliAnalysisTaskSEDplusCorrelations::~AliAnalysisTaskSEDplusCorrelations()
171{
172 //
173 // Destructor
174 //
cb7c2594 175 if(fOutput) {delete fOutput; fOutput = 0;}
176 if(fListCuts) {delete fListCuts; fListCuts = 0;}
177 if(fRDCutsAnalysis) {delete fRDCutsAnalysis; fRDCutsAnalysis = 0;}
178 if(fCuts) {delete fCuts; fCuts = 0;}
179 if(fCounter) {delete fCounter; fCounter = 0;}
180 if(fCorrelator) {delete fCorrelator; fCorrelator=0;}
181
35151011 182
183
184}
185//_________________________________________________________________
186void AliAnalysisTaskSEDplusCorrelations::SetMassLimits(Float_t range){
187 // set invariant mass limits
188 Float_t bw=GetBinWidth();
189 fUpmasslimit = 1.865+range;
190 fLowmasslimit = 1.865-range;
191 SetBinWidth(bw);
192}
193//_________________________________________________________________
194void AliAnalysisTaskSEDplusCorrelations::SetMassLimits(Float_t lowlimit, Float_t uplimit){
195 // set invariant mass limits
196 if(uplimit>lowlimit)
197 {
198 Float_t bw=GetBinWidth();
199 fUpmasslimit = lowlimit;
200 fLowmasslimit = uplimit;
201 SetBinWidth(bw);
202 }
203}
204//________________________________________________________________
205void AliAnalysisTaskSEDplusCorrelations::SetBinWidth(Float_t w){
206 // Define width of mass bins
207 Float_t width=w;
208 Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
209 Int_t missingbins=4-nbins%4;
210 nbins=nbins+missingbins;
211 width=(fUpmasslimit-fLowmasslimit)/nbins;
212 if(missingbins!=0){
213 printf("AliAnalysisTaskSEDplusCorrelations::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
214 }
215 else{
216 if(fDebug>1) printf("AliAnalysisTaskSEDplusCorrelations::SetBinWidth: width set to %f\n",width);
217 }
218 fBinWidth=width;
219}
220//_________________________________________________________________
221Int_t AliAnalysisTaskSEDplusCorrelations::GetNBinsHistos(){
222 // Compute number of mass bins
223 return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
224}
225
226
227//__________________________________________
228void AliAnalysisTaskSEDplusCorrelations::Init(){
229 //
230 // Initialization
231 //
232 if(fDebug > 1) printf("AnalysisTaskSEDplusCorrelations::Init() \n");
233
cb7c2594 234
35151011 235 fListCuts=new TList();
cb7c2594 236
35151011 237 // fListCutsAsso=new TList();
238
239 AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
240 analysis->SetName("AnalysisCuts");
241
242 // AliHFAssociatedTrackCuts *trkcuts = new AliHFAssociatedTrackCuts(*fCuts);
243 //trkcuts->SetName("Assotrkcuts");
244
245 fListCuts->Add(analysis);
cb7c2594 246 //fListCutsAsso->Add(trkcuts);
35151011 247
248
249
250 PostData(2,fListCuts);
cb7c2594 251 //ŧ PostData(4,fListCutsAsso);
35151011 252
253 return;
254}
255
256//________________________________________________________________________
257void AliAnalysisTaskSEDplusCorrelations::UserCreateOutputObjects()
258{
259 // Create the output container
260 //
cb7c2594 261
35151011 262 if(fDebug > 1) printf("AnalysisTaskSEDplusCorrelations::UserCreateOutputObjects() \n");
cb7c2594 263 // Correlator creation and definition
35151011 264
265 Double_t Pi = TMath::Pi();
cb7c2594 266
267 // Set up the correlator
268
269
35151011 270 fCorrelator = new AliHFCorrelator("Correlator",fCuts,fSystem); // fCuts is the hadron cut object, fSystem to switch between pp or PbPb
cb7c2594 271
35151011 272 fCorrelator->SetDeltaPhiInterval((-0.5-1./32)*Pi,(1.5-1./32)*Pi); // set correct phi interval
cb7c2594 273
35151011 274 fCorrelator->SetEventMixing(fMixing); //set kFALSE/kTRUE for mixing Off/On
cb7c2594 275
35151011 276 fCorrelator->SetAssociatedParticleType(fSelect); // set 1/2/3 for hadron/kaons/kzeros
cb7c2594 277
35151011 278 fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
cb7c2594 279
35151011 280 fCorrelator->SetUseMC(fReadMC);
34538691 281 fCorrelator->SetPIDmode(2);
cb7c2594 282 fCorrelator->SetUseReco(fReco); // to analyse reco/kine
283
284 // For Event Mixing
285
35151011 286 Bool_t pooldef = fCorrelator->DefineEventPool();
287
288 if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
289
290
291 // Several histograms are more conveniently managed in a TList
292 fOutput = new TList();
293 fOutput->SetOwner();
294 fOutput->SetName("OutputHistos");
295
296 TString hisname;
297 Int_t index=0;
298 Int_t nbins=GetNBinsHistos();
8d35b368 299
300
35151011 301
302 for(Int_t i=0;i<fNPtBins;i++){
303
304 index=GetHistoIndex(i);
35151011 305
35151011 306 hisname.Form("hMassPtK0S%d",i);
307 fMassHistK0S[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.3,0.8);
308 fMassHistK0S[index]->Sumw2();
309
310 hisname.Form("hLeadPt%d",i);
311 fLeadPt[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.0,50);
312 fLeadPt[index]->Sumw2();
313
cb7c2594 314
35151011 315 hisname.Form("hMassPt%d",i);
316 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
317 fMassHist[index]->Sumw2();
cb7c2594 318
319
35151011 320 hisname.Form("hMassPt%dTC",i);
321 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
322 fMassHistTC[index]->Sumw2();
323
cb7c2594 324
35151011 325 index=GetSignalHistoIndex(i);
35151011 326
35151011 327 hisname.Form("hSigPtK0S%d",i);
328 fMassHistK0S[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.3,0.8);
329 fMassHistK0S[index]->Sumw2();
330
331 hisname.Form("hSigLeadPt%d",i);
332 fLeadPt[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.0,50);
333 fLeadPt[index]->Sumw2();
334
335 hisname.Form("hSigPt%d",i);
336 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
337 fMassHist[index]->Sumw2();
338
cb7c2594 339 hisname.Form("hSigOrigCPt%d",i);
340 fMassHistOrigC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
341 fMassHistOrigC[index]->Sumw2();
35151011 342
cb7c2594 343 hisname.Form("hSigOrigBPt%d",i);
344 fMassHistOrigB[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
345 fMassHistOrigB[index]->Sumw2();
35151011 346
cb7c2594 347
348 hisname.Form("hSigMCPt%d",i);
349 fMassHistMC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
350 fMassHistMC[index]->Sumw2();
35151011 351
cb7c2594 352 }
35151011 353
354
355 for(Int_t i=0; i<3*fNPtBins; i++){
cb7c2594 356
35151011 357 fOutput->Add(fMassHistK0S[i]);
358 fOutput->Add(fLeadPt[i]);
359 fOutput->Add(fMassHist[i]);
360 fOutput->Add(fMassHistTC[i]);
cb7c2594 361 fOutput->Add(fMassHistOrigC[i]);
362 fOutput->Add(fMassHistOrigB[i]);
363 fOutput->Add(fMassHistMC[i]);
364
35151011 365 }
35151011 366
35151011 367
cb7c2594 368 TH1D * EventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
369 if(fReadMC) fOutput->Add(EventTypeMC);
370
371
372 //Event Mixing histos (for control plots)
373
374 Int_t NumberEvents = fCuts->GetMaxNEventsInPool();
375 Int_t NofTracks = fCuts->GetMinNTracksInPool();
376 Int_t NofCentBins = fCuts->GetNCentPoolBins();
377 Double_t * CentBins = fCuts->GetCentPoolBins();
378 Int_t NumberZVtxBins = fCuts->GetNZvtxPoolBins();
379 Double_t *ZVtxBins = fCuts->GetZvtxPoolBins();
380
381
382
383 fEventMix = new TH2D("EventMixingCheck","EventMixingCheck",100,0,600,100,-15,15);
384
385 fEventMix->GetXaxis()->SetTitle("Multiplicity ");
386 fEventMix->GetYaxis()->SetTitle("Z vertex [cm]");
387
35151011 388 if(fMixing)fOutput->Add(fEventMix);
35151011 389
cb7c2594 390 TH3D * EventsInPool = new TH3D("EventsInPool","Number of events in pool",NofCentBins,0,600,NumberZVtxBins,-15,15,1000000,0,1000000);
391
392 EventsInPool->GetXaxis()->SetTitle("Multiplicity ");
393 EventsInPool->GetYaxis()->SetTitle("Z vertex [cm]");
394 EventsInPool->GetZaxis()->SetTitle("Number of events in pool");
395 if(fMixing) fOutput->Add(EventsInPool);
396
397 Int_t MaxNofTracks = (NumberEvents+1)*NofTracks;
398
399
400 TH3D * NofTracksInPool = new TH3D("NofTracksInPool","Number of tracks in pool",NofCentBins,0,500,NumberZVtxBins,-15,15,MaxNofTracks,0,MaxNofTracks);
401 NofTracksInPool->GetXaxis()->SetTitle("Multiplicity ");
402 NofTracksInPool->GetYaxis()->SetTitle("Z vertex [cm]");
403 NofTracksInPool->GetZaxis()->SetTitle("Number of tracks");
404
405 if(fMixing) fOutput->Add(NofTracksInPool);
406
407 TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Calls per pool bin",NofCentBins,CentBins,NumberZVtxBins,ZVtxBins);
408 NofPoolBinCalls->GetXaxis()->SetTitle("Multiplicity ");
409 NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]");
410 if(fMixing) fOutput->Add(NofPoolBinCalls);
411
412
413
414 fHistNEvents = new TH1F("fHistNEvents", "number of events ",10,-0.5,10.5);
35151011 415 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
416 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents accepted");
cb7c2594 417 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvent with good vertex");
418 fHistNEvents->GetXaxis()->SetBinLabel(4,"Total number of candidates");
419 fHistNEvents->GetXaxis()->SetBinLabel(5,"Number without bitmask");
420 fHistNEvents->GetXaxis()->SetBinLabel(6,"Number after Physics Selection");
421 fHistNEvents->GetXaxis()->SetBinLabel(7,"Total number in Fiducial Accpt");
422 fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of good candidates");
423
35151011 424
425 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
426 fHistNEvents->Sumw2();
427 fHistNEvents->SetMinimum(0);
428 fOutput->Add(fHistNEvents);
cb7c2594 429
35151011 430 // Counter for Normalization
431 TString normName="NormalizationCounter";
432 AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
cb7c2594 433
35151011 434 if(cont)normName=(TString)cont->GetName();
435 fCounter = new AliNormalizationCounter(normName.Data());
436 fCounter->Init();
437
cb7c2594 438 CreateCorrelationObjs();
35151011 439
440 PostData(1,fOutput);
441 PostData(3,fCounter);
442 return;
443}
444
445//________________________________________________________________________
446void AliAnalysisTaskSEDplusCorrelations::UserExec(Option_t */*option*/)
447{
448 // Do the analysis
449 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
450 TClonesArray *array3Prong = 0;
451
cb7c2594 452 if(fReco) std::cout << "USING RECONSTRUCTION" << std::endl;
453 if(!fReco) std::cout << "USING MC TRUTH" << std::endl;
454
35151011 455 if(!fMixing){
456 if(fSelect==1) cout << "TASK::Correlation with hadrons on SE "<< endl;
457 if(fSelect==2) cout << "TASK::Correlation with kaons on SE "<< endl;
458 if(fSelect==3) cout << "TASK::Correlation with kzeros on SE "<< endl;
cb7c2594 459 }
35151011 460 if(fMixing){
461 if(fSelect==1) cout << "TASK::Correlation with hadrons on ME "<< endl;
462 if(fSelect==2) cout << "TASK::Correlation with kaons on ME "<< endl;
463 if(fSelect==3) cout << "TASK::Correlation with kzeros on ME "<< endl;
464 }
465
466
467 if(!aod && AODEvent() && IsStandardAOD()) {
468
469 aod = dynamic_cast<AliAODEvent*> (AODEvent());
470 AliAODHandler* aodHandler = (AliAODHandler*)
471 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
472 if(aodHandler->GetExtensions()) {
473 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
474 AliAODEvent *aodFromExt = ext->GetAOD();
475 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
476 }
477 } else if(aod) {
478 array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
479 }
480
481 if(!aod || !array3Prong) {
482 printf("AliAnalysisTaskSEDplusCorrelations::UserExec: Charm3Prong branch not found!\n");
483 return;
484 }
485
486
487 // the AODs with null vertex pointer didn't pass the PhysSel
cb7c2594 488
35151011 489 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
490 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC);
cb7c2594 491
35151011 492 fHistNEvents->Fill(0); // count event
493
494
495 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
496
497
498
499 PostData(1,fOutput);
cb7c2594 500
35151011 501 if(!isEvSel)return;
cb7c2594 502
35151011 503 fHistNEvents->Fill(1);
504
34538691 505 // set PIDResponse for associated tracks
cb7c2594 506
34538691 507 fCorrelator->SetPidAssociated();
35151011 508
509 TClonesArray *arrayMC=0;
510 AliAODMCHeader *mcHeader=0;
cb7c2594 511
35151011 512 // AOD primary vertex
513 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
cb7c2594 514
515
35151011 516 TString primTitle = vtx1->GetTitle();
cb7c2594 517
518 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0){
519 fHistNEvents->Fill(2);
520 }
521
522
523
524
35151011 525 // load MC particles
526 if(fReadMC){
527
528 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
529 if(!arrayMC) {
530 printf("AliAnalysisTaskSEDplusCorrelations::UserExec: MC particles branch not found!\n");
531 return;
532 }
533
534 // load MC header
535 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
536 if(!mcHeader) {
537 printf("AliAnalysisTaskSEDplusCorrelations::UserExec: MC header branch not found!\n");
538 return;
539 }
cb7c2594 540
541 Bool_t isMCeventgood = kFALSE;
542
543
544 Int_t eventType = mcHeader->GetEventType();
545 Int_t NMCevents = fCuts->GetNofMCEventType();
546
547
548 for(Int_t k=0; k<NMCevents; k++){
549 Int_t * MCEventType = fCuts->GetMCEventType();
550
551 if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
552
553 }
554 ((TH1D*)fOutput->FindObject("EventTypeMC"))->Fill(eventType);
555
556 if(NMCevents && !isMCeventgood){
557
558 std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
559 return;
560 }
561
35151011 562 }
563
564 //HFCorrelators initialization (for this event)
565
cb7c2594 566 fCorrelator->SetAODEvent(aod); // set the event to be processed
567
35151011 568 Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing
cb7c2594 569
35151011 570 if(!correlatorON) {
571 AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
572 return;
573 }
574 if(fReadMC) fCorrelator->SetMCArray(arrayMC);
575
cb7c2594 576 Int_t n3Prong = 0;
35151011 577
cb7c2594 578 if(fReco) n3Prong = array3Prong->GetEntriesFast();
579
580
581 if(!fReco) n3Prong = arrayMC->GetEntriesFast(); //for MC kine
582
583
8d35b368 584 printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
cb7c2594 585
35151011 586 Int_t index;
587 Int_t pdgDgDplustoKpipi[3]={321,211,211};
cb7c2594 588
35151011 589 Int_t nSelectedloose=0,nSelectedtight=0;
590
cb7c2594 591 Double_t ptCand;
592 Double_t phiCand;
593 Double_t etaCand;
594 Double_t invMass = -1;
595
596 Int_t nIDs[3] = {-9999999};
35151011 597
cb7c2594 598 Bool_t isDplus = kFALSE;
599
600 Int_t labDp = -1;
601
35151011 602 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
cb7c2594 603
604 AliAODMCParticle* DplusMC;
605
606 if(fReco){
35151011 607 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
cb7c2594 608
609 fHistNEvents->Fill(3);
610
611
612
613 if(d->Pt()<2.0) continue; // Dplus below 2.0 not considered.
614
615
616
617 // cout<<multipli<<" multi"<<endl;
618
35151011 619
cb7c2594 620
35151011 621 if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
cb7c2594 622 fHistNEvents->Fill(4);
35151011 623 continue;
624 }
625
626 Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
cb7c2594 627
628 if(!fRDCutsAnalysis->GetIsSelectedCuts()) continue;
629
630 fHistNEvents->Fill(5);
631
632
35151011 633 if(fReadMC){
634 labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
cb7c2594 635
636 // cout<<labDp<<"***"<<endl;
35151011 637 if(labDp>=0){
638 isDplus = kTRUE;
35151011 639
cb7c2594 640 }
641 }
642
643 invMass=d->InvMassDplus();
35151011 644 Double_t rapid=d->YDplus();
cb7c2594 645 etaCand = d->Eta();
646 phiCand = d->Phi();
647 ptCand = d->Pt();
648
35151011 649 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
cb7c2594 650
35151011 651 if(isFidAcc){
cb7c2594 652 fHistNEvents->Fill(6);
653 nSelectedloose++;}
654 if(!isFidAcc) continue;
655 if(!passTightCuts)continue;
656
657 fHistNEvents->Fill(7);
658
659 nSelectedtight++;
660
661
662 nIDs[0] = ((AliAODTrack*)d->GetDaughter(0))->GetID();
663 nIDs[1] = ((AliAODTrack*)d->GetDaughter(1))->GetID();
664 nIDs[2] = ((AliAODTrack*)d->GetDaughter(2))->GetID();
665
35151011 666 }
cb7c2594 667
668 Int_t labDplus=-1;
669
670 if(!fReco){
671
672 DplusMC = dynamic_cast<AliAODMCParticle*>(arrayMC->At(i3Prong));
673
674 if (!DplusMC) {
675 AliWarning("Dplus MC Particle not found ");
676
677 continue;
678 }
679 labDplus = DplusMC->GetLabel();
680
681 Int_t PDG =TMath::Abs(DplusMC->PdgCode());
682
683 if(PDG !=411) continue; // not a Dplus
684 ptCand = DplusMC->Pt();
685 phiCand = DplusMC->Phi();
686 etaCand = DplusMC->Eta();
687
688 Bool_t isFidAccMC =fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,DplusMC->Y());
689
690 if(!isFidAccMC)continue;
35151011 691
cb7c2594 692 }
693
694 //cout << "PHI D+ = " << phiCand << endl;
695 // cout << "ETA D+ = " << etaCand << endl;
696 // cout << "PT D+ = " << ptCand << endl;
697
698
699
700
701 Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
702
35151011 703 if(iPtBin>=0){
cb7c2594 704
35151011 705 index=GetHistoIndex(iPtBin);
cb7c2594 706
707 cout<<"*****"<<iPtBin<<endl;
708 Double_t effTrig;
709 if (fTrig){
35151011 710
cb7c2594 711 Double_t multipli = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1,1.);
712
713 effTrig = fCuts->GetTrigWeight(ptCand,multipli);
35151011 714
cb7c2594 715 // cout<<"*****"<<effTrig<<" "<<multipli<<" "<<iPtBin<<endl;
716
717 }
718 else
719 {
720
721 effTrig = 1.0;
722
723 }
724
725
726 Double_t trigweight = 1/effTrig ;
727
728 // cout<<"****"<<trigweight<<" ***"<<endl;
729
730 if(fReco && !fMixing){
731
732 fMassHist[index]->Fill(invMass); //without trig weight
733
734 fMassHistTC[index]->Fill(invMass,trigweight); //with trig wt
735 }
736
737 if( fReco && fReadMC && isDplus && !fMixing) {
738
739 index=GetSignalHistoIndex(iPtBin);
740
741 fMassHistMC[index]->Fill(invMass,trigweight);
35151011 742
cb7c2594 743 if(labDp>=0){
744
745 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
746 Int_t DplusSource = CheckOrigin(arrayMC,partDp);
747
748 if(DplusSource == 4){ // is from charm
35151011 749
cb7c2594 750 fMassHistOrigC[index]->Fill(invMass,trigweight);
751 }
752
753 if(DplusSource == 5){ // is from beauty
35151011 754
cb7c2594 755 fMassHistOrigB[index]->Fill(invMass,trigweight);
35151011 756
757
cb7c2594 758 }
759 }
760 }
761
762
763 if(!fReco && fReadMC && !fMixing){
35151011 764
cb7c2594 765
766 if(labDplus>=0){
767
768 fMassHist[index]->Fill(1.869);
769
770 AliAODMCParticle *kineDp = (AliAODMCParticle*)arrayMC->At(labDplus);
771
772 Int_t DplusSource = CheckOrigin(arrayMC,kineDp);
773
774
775 if(DplusSource==4){ // is from charm
776
777 ((TH1F*)fOutput->FindObject(Form("histOrig_Dplus_Bin%d",iPtBin)))->Fill(0.);
35151011 778 }
35151011 779
cb7c2594 780 if(DplusSource==5){ // is from beauty
35151011 781
cb7c2594 782
783 ((TH1F*)fOutput->FindObject(Form("histOrig_Dplus_Bin%d",iPtBin)))->Fill(1.);
784 }
785 }
786 }
35151011 787
cb7c2594 788
789 //Dplus info
35151011 790
cb7c2594 791 Double_t phiDplus = fCorrelator->SetCorrectPhiRange(phiCand);
792 fCorrelator->SetTriggerParticleProperties(ptCand,phiDplus,etaCand);
793
794
795
796 Double_t ptlead = 0;
797 Double_t philead = 0;
798 Double_t etalead = 0;
799 Double_t refpt = 0;
800
801 Int_t LeadSource = 0;
802
803
804
805 Bool_t execPool = fCorrelator->ProcessEventPool();
806
807 //printf("*************: %d\n",execPool);
808
809 if(fMixing && !execPool) {
810 AliInfo("Mixed event analysis: pool is not ready");
811 continue;
812 }
813
814 Int_t NofEventsinPool = 1;
815 if(fMixing) {
816 NofEventsinPool = fCorrelator->GetNofEventsInPool();
35151011 817
cb7c2594 818 //cout<<"*********"<<NofEventsinPool<<endl;
819 }
820
821
822
823 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, it stops at 1
824
825 Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix);
826 if(!analyzetracks) {
827 AliInfo("AliHFCorrelator::Cannot process the track array");
828 continue;
829 }
830
35151011 831
cb7c2594 832 //Int_t NofTracks = fCorrelator->GetNofTracks();
833
834 //cout<<"***number of tracks****"<<fCorrelator->GetNofTracks()<<endl;
835
836 for (Int_t iTrack = 0;iTrack<fCorrelator->GetNofTracks();iTrack++) {
837 Bool_t runcorrelation = fCorrelator->Correlate(iTrack);
838
839 if(!runcorrelation) continue;
840
841 Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
842
843
844 Double_t DeltaEta = fCorrelator->GetDeltaEta();
845
846 AliReducedParticle* redpart = fCorrelator->GetAssociatedParticle();
847 Double_t phiHad=redpart->Phi();
848 Double_t ptHad=redpart->Pt();
849 Double_t etaHad=redpart->Eta();
850 Int_t label = redpart->GetLabel();
851 Int_t trackid = redpart->GetID();
852
853 Double_t effi = redpart->GetWeight();
854 Double_t weight = (1/effi)*(trigweight);
855
856 phiHad = fCorrelator->SetCorrectPhiRange(phiHad);
857 // cout<<effi<<"******"<<endl;
858
35151011 859
35151011 860
cb7c2594 861 if (!fMixing && fReco){
862 if( trackid == nIDs[0] || trackid == nIDs[1] || trackid == nIDs[2]) continue; //discard the Dplus daughters
863 }
864
865 if(fReco && !fReadMC)
866 {
867 FillCorrelations(ptHad,invMass,DeltaPhi,DeltaEta,iPtBin,fSelect,weight);
868 }
869
870
871 if(fReco && fReadMC && isDplus){
872
873 if(label<0) continue;
874
875 AliAODMCParticle* mchad = (AliAODMCParticle*)arrayMC->At(label);
876
877 if (!mchad) continue;
878
879 if (!mchad->IsPhysicalPrimary()) continue; //reject the Reco track if correspondent Kine track is not primary
880
881 if(labDp < 0) continue;
882
883 AliAODMCParticle *reDp = (AliAODMCParticle*)arrayMC->At(labDp);
884
885 Int_t RDSource = CheckOrigin(arrayMC,reDp);
886
887 Int_t MCSource = CheckTrackOrigin(arrayMC,(AliAODMCParticle*)arrayMC->At(label)) ;// check source of associated particle (hadron/kaon/K0)
8d35b368 888
35151011 889
cb7c2594 890 FillMCRCorrelations(ptHad,invMass,DeltaPhi,DeltaEta,iPtBin,MCSource,RDSource,weight);
891
35151011 892 }
cb7c2594 893
894 // montecarlo kine
895 if( !fReco && fReadMC){
35151011 896
cb7c2594 897 if(label<0) continue;
35151011 898
cb7c2594 899 if(TMath::Abs(etaHad) > 0.8) continue;
35151011 900
cb7c2594 901 AliAODMCParticle *hadMC = (AliAODMCParticle*)arrayMC->At(label);
902 if(!hadMC) continue;
903 if (!hadMC->IsPhysicalPrimary()) continue;
904 if(labDplus<0) continue;
905
906 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDplus);
907 if(IsDDaughter(partDp,hadMC,arrayMC)) continue;
908
35151011 909
cb7c2594 910 AliAODMCParticle *trueDp = (AliAODMCParticle*)arrayMC->At(labDplus);
911 Int_t DSource = CheckOrigin(arrayMC,trueDp);
912
913 Int_t PartSource = CheckTrackOrigin(arrayMC,(AliAODMCParticle*)arrayMC->At(label)) ;// check source of associated particle (hadron/kaon/K0)
35151011 914
915
cb7c2594 916 FillMCTruthCorrelations(ptHad,DeltaPhi,DeltaEta,iPtBin,PartSource,DSource,fSelect);
917 }
918
919
920
921 // For leading particle
922
923 if (ptHad > refpt) {
924 //refpt = ptHad; ptlead = ptHad;
925
926 if(fReadMC){
927 if(label<0) continue;
928 if(!(AliAODMCParticle*)arrayMC->At(label))continue;
929
930 LeadSource = CheckTrackOrigin(arrayMC,(AliAODMCParticle*)arrayMC->At(label)) ;// check source of associated particle (hadron/kaon/K0)
931 }
932
933 philead = phiCand - phiHad;
934 etalead = etaCand - etaHad;
935 if (philead < (-1)*TMath::Pi()/2) philead += 2*TMath::Pi();
936
937 if (philead > 3*TMath::Pi()/2) philead -= 2*TMath::Pi();
938
939 refpt = ptHad; ptlead = ptHad;
940 }
941 }
942
943 /*if (fReco){
944 Double_t fillSparseLeadDplus[3] = {philead,invMass,etalead};
945 }
946 else{
947 Double_t fillSparseLeadDplus[3] = {philead,1.87,etalead};
948 }
949
950
951 if(fReco && !fReadMC && !fMixing){
952
953 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Lead_Bin%d",iPtBin)))->Fill(fillSparseLeadDplus,eweight);
954
955 }
956
957 if(fReco && fReadMC && isDplus && !fMixing){
958
959 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Lead_Bin%d",iPtBin)))->Fill(fillSparseLeadDplus,eweight);
960
961
962
963 if(LeadSource>=1&&LeadSource<=4){ // is from charm
964
965 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Lead_From_c_Bin%d",iPtBin)))->Fill(fillSparseLeadDplus,eweight);
966
967
968 }
969 if(LeadSource>=4&&LeadSource<=8){ // is from beauty
970
971
972 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Lead_From_b_Bin%d",iPtBin)))->Fill(fillSparseLeadDplus,eweight);
973
974
975
976 }
977
978 // from NHF
979 if(LeadSource==0){ // is from charm ->D
980
981
982 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Lead_From_nhf_Bin%d",iPtBin)))->Fill(fillSparseLeadDplus,eweight);
983
984 }
985
986 fLeadPt[index]->Fill(ptlead);
987 }
988
989 if(!fReco && fReadMC && !fMixing){
990
991 index=GetSignalHistoIndex(iPtBin);
992
993
994 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Lead_Bin%d",iPtBin)))->Fill(fillSparseLeadDplus,eweight);
995
996
997 if(LeadSource>=1&&LeadSource<=4){ // is from charm
998
999 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Lead_From_c_Bin%d",iPtBin)))->Fill(fillSparseLeadDplus,eweight);
1000
1001
1002
1003 }
1004 if(LeadSource>=4&&LeadSource<=8){ // is from beauty
1005
1006((THnSparseF*)fOutput->FindObject(Form("hPhi_Lead_From_b_Bin%d",iPtBin)))->Fill(fillSparseLeadDplus,eweight);
1007
1008 }
1009
1010 // from NHF
1011 if(LeadSource==0){ // is from charm ->D
1012
1013((THnSparseF*)fOutput->FindObject(Form("hPhi_Lead_From_nhf_Bin%d",iPtBin)))->Fill(fillSparseLeadDplus,eweight);
1014
1015
1016 }
1017 fLeadPt[index]->Fill(ptlead);
1018
1019 }*/
1020
1021 } //jmix
35151011 1022
cb7c2594 1023
1024 } //pt bin
1025
1026
35151011 1027
1028 }
1029
1030 if(fMixing){
1031 Bool_t updated = fCorrelator->PoolUpdate();
1032
1033 if(!updated) AliInfo("Pool was not updated");
1034 }
1035 fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1036 fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1037
1038
1039 PostData(1,fOutput);
1040 PostData(2,fListCuts);
1041 PostData(3,fCounter);
cb7c2594 1042 // PostData(4,fListCutsAsso);
35151011 1043 return;
1044}
1045
1046//________________________________________________________________________
cb7c2594 1047
35151011 1048void AliAnalysisTaskSEDplusCorrelations::Terminate(Option_t */*option*/)
1049{
1050 // Terminate analysis
1051 //
1052 if(fDebug > 1) printf("AnalysisTaskSEDplusCorrelations: Terminate() \n");
1053
1054 fOutput = dynamic_cast<TList*> (GetOutputData(1));
1055 if (!fOutput) {
1056 printf("ERROR: fOutput not available\n");
1057 return;
1058 }
35151011 1059
1060 return;
1061}
1062
1063
1064//________________________________________________________________________
cb7c2594 1065void AliAnalysisTaskSEDplusCorrelations::FillCorrelations(Double_t ptTrack, Double_t mass, Double_t deltaPhi, Double_t deltaEta, Int_t ind, Int_t sel, Double_t eweight) const{
35151011 1066
35151011 1067
cb7c2594 1068 //Filling THnSparse for correlations , only for charged tracks
1069
1070
35151011 1071 if(sel==1){
cb7c2594 1072
1073 if(!fMixing){
1074
1075 Double_t ptLim_Sparse = ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d",ind)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1076 Double_t EtaLim_Sparse = ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d",ind)))->GetAxis(3)->GetXmax();
1077
1078 if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01; if(deltaEta > EtaLim_Sparse) deltaEta = EtaLim_Sparse-0.01;
1079 if(deltaEta < -EtaLim_Sparse) deltaEta = -EtaLim_Sparse+0.01;
1080
1081 //variables for filling histos
1082
1083 Double_t fillSparseDplus[4] = {deltaPhi,mass,ptTrack,deltaEta};
1084
1085
1086
1087
1088 //sparse fill for data (tracks) + weighted
1089
1090 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d",ind)))->Fill(fillSparseDplus,eweight);
1091
1092 }
1093
1094 if(fMixing) {
1095
1096
1097 //variables for filling histos
1098
1099 Double_t fillSparseDplusMix[3] = {deltaPhi,mass,deltaEta};
1100
1101
1102 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ind)))->Fill(fillSparseDplusMix,eweight);
1103 }
1104
35151011 1105 }
cb7c2594 1106
35151011 1107 if(sel==2){
cb7c2594 1108
35151011 1109 }
1110 if(sel==3){
cb7c2594 1111
35151011 1112 }
cb7c2594 1113
35151011 1114 return;
1115}
1116
1117
1118
1119//________________________________________________________________________
35151011 1120
cb7c2594 1121
1122
1123void AliAnalysisTaskSEDplusCorrelations::FillMCTruthCorrelations( Double_t ptTrack,Double_t deltaPhi, Double_t deltaEta, Int_t ind, Int_t mcSource, Int_t origDplus, Int_t sel) const{
1124
1125 // Filling histos with MC Kine information
1126
1127 Double_t invm = 1.876;
1128
1129 Double_t ptLim_Sparse = ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d",ind)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1130 Double_t EtaLim_Sparse = ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d",ind)))->GetAxis(3)->GetXmax();
1131 if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
1132 if(deltaEta > EtaLim_Sparse) deltaEta = EtaLim_Sparse-0.01;
1133 if(deltaEta < -EtaLim_Sparse) deltaEta = -EtaLim_Sparse+0.01;
1134
1135
35151011 1136
cb7c2594 1137 //variables for filling histos
1138
1139 Double_t fillSparseDplus[4] = {deltaPhi,invm,ptTrack,deltaEta};
1140
35151011 1141
1142 if(sel==1){
cb7c2594 1143
1144 if(!fMixing){
1145
1146 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d",ind)))->Fill(fillSparseDplus);
1147
1148 if(origDplus==4&&mcSource>=1&&mcSource<=3){ // is from charm
1149
1150 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_HFC_Bin%d",ind)))->Fill(fillSparseDplus);
1151 }
1152 if(origDplus==5&&mcSource>=4&&mcSource<=8){ // is from beauty
1153
1154 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_HFB_Bin%d",ind)))->Fill(fillSparseDplus);
1155 }
1156
1157 if(mcSource==0){ // is from NHF
1158
1159 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_NHF_Bin%d",ind)))->Fill(fillSparseDplus);
1160 }
1161
1162 }
1163
1164 if(fMixing) {
35151011 1165
cb7c2594 1166
1167 //variables for filling histos
1168 Double_t fillSparseDplusMix[3] = {deltaPhi,invm,deltaEta};
1169
1170 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ind)))->Fill(fillSparseDplusMix);
1171 }
1172
1173 }
1174
1175 if(sel==2){
1176
1177 }
1178
1179 if(sel==3){
1180
1181 }
1182
1183
1184 return;
1185}
1186
1187
1188void AliAnalysisTaskSEDplusCorrelations::FillMCRCorrelations(Double_t ptTrack, Double_t mass, Double_t deltaPhi, Double_t deltaEta, Int_t ind, Int_t mcSource, Int_t origDplus, Double_t eweight) const{
1189
1190 // Filling histos with MC information
1191
1192
1193 Double_t ptLim_Sparse = ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d",ind)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1194 Double_t EtaLim_Sparse = ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d",ind)))->GetAxis(3)->GetXmax();
1195 if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
1196 if(deltaEta > EtaLim_Sparse) deltaEta = EtaLim_Sparse-0.01;
1197 if(deltaEta < -EtaLim_Sparse) deltaEta = -EtaLim_Sparse+0.01;
1198
1199
1200 if(!fMixing){
1201
1202 Double_t fillSparseDplus[4] = {deltaPhi,mass,ptTrack,deltaEta};
1203
1204
1205
1206 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d",ind)))->Fill(fillSparseDplus,eweight);
1207
1208 if(origDplus==4&&mcSource>=1&&mcSource<=3){ // is from charm
1209
1210 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_HFC_Bin%d",ind)))->Fill(fillSparseDplus,eweight);
1211 }
1212 if(origDplus==5&&mcSource>=4&&mcSource<=8){ // is from beauty
1213
1214 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_HFB_Bin%d",ind)))->Fill(fillSparseDplus,eweight);
1215 }
1216
1217 if(mcSource==0){ // is from NHF
1218
1219 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_NHF_Bin%d",ind)))->Fill(fillSparseDplus,eweight);
1220 }
1221 }
1222
1223 if(fMixing) {
1224
1225 //variables for filling histos
1226 Double_t fillSparseDplusMix[3] = {deltaPhi,mass,deltaEta};
1227
1228 ((THnSparseF*)fOutput->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ind)))->Fill(fillSparseDplusMix,eweight);
1229 }
1230
1231
1232 return;
1233}
1234
1235Bool_t AliAnalysisTaskSEDplusCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* mcArray) const {
1236
1237 //Daughter removal in MCKine case
1238 Bool_t isDaughter = kFALSE;
1239 Int_t labelDplus = d->GetLabel();
1240
1241 Int_t mother = track->GetMother();
1242
1243
1244 while (mother > 0){
1245 AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1246 if (mcMoth){
1247 if (mcMoth->GetLabel() == labelDplus) isDaughter = kTRUE;
1248 mother = mcMoth->GetMother(); //goes back to last generation
1249 } else{
1250 AliError(" mother particle not found!");
1251 break;
1252 }
1253 }
1254
1255 return isDaughter;
1256}
1257
1258
1259Int_t AliAnalysisTaskSEDplusCorrelations::CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1260 //
1261 // checks on particle (#) origin:
1262 // 0) Not HF
1263 // 1) D->#
1264 // 2) D->X->#
1265 // 3) c hadronization
1266 // 4) B->#
1267 // 5) B->X-># (X!=D)
1268 // 6) B->D->#
1269 // 7) B->D->X->#
1270 // 8) b hadronization
1271 //
1272 if(fDebug>2) printf("AliAnalysisTaskSEDplusCorrelations::CheckTrkOrigin() \n");
35151011 1273
cb7c2594 1274 Int_t pdgGranma = 0;
1275 Int_t mother = 0;
1276 mother = mcPartCandidate->GetMother();
1277 Int_t istep = 0;
1278 Int_t abspdgGranma =0;
1279 Bool_t isFromB=kFALSE;
1280 Bool_t isDdaugh=kFALSE;
1281 Bool_t isDchaindaugh=kFALSE;
1282 Bool_t isBdaugh=kFALSE;
1283 Bool_t isBchaindaugh=kFALSE;
1284 Bool_t isQuarkFound=kFALSE;
1285
1286 while (mother > 0){
1287 istep++;
1288 AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1289 if (mcMoth){
1290 pdgGranma = mcMoth->GetPdgCode();
1291 abspdgGranma = TMath::Abs(pdgGranma);
1292 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1293 isBchaindaugh=kTRUE;
1294 if(istep==1) isBdaugh=kTRUE;
1295 }
1296 if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)){
1297 isDchaindaugh=kTRUE;
1298 if(istep==1) isDdaugh=kTRUE;
1299 }
1300 if(abspdgGranma==4 || abspdgGranma==5) {isQuarkFound=kTRUE; if(abspdgGranma==5) isFromB = kTRUE;}
1301 mother = mcMoth->GetMother();
1302 }else{
1303 AliError("Failed casting the mother particle!");
1304 break;
35151011 1305 }
cb7c2594 1306 }
1307
1308 //decides what to return based on the flag status
1309 if(isQuarkFound) {
1310 if(!isFromB) { //charm
1311 if(isDdaugh) return 1; //charm immediate
1312 else if(isDchaindaugh) return 2; //charm chain
1313 else return 3; //charm hadronization
1314 }
1315 else { //beauty
1316 if(isBdaugh) return 4; //b immediate
1317 else if(isBchaindaugh) { //b chain
1318 if(isDchaindaugh) {
1319 if(isDdaugh) return 6; //d immediate
1320 return 7; //d chain
1321 }
1322 else return 5; //b, not d
1323 }
1324 else return 8; //b hadronization
35151011 1325 }
1326 }
cb7c2594 1327 else return 0; //no HF quark
1328}
35151011 1329
cb7c2594 1330Int_t AliAnalysisTaskSEDplusCorrelations::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1331 //
1332 // checking whether the mother of the particles come from a charm or a bottom quark
1333 //
1334
1335 Int_t pdgGranma = 0;
1336 Int_t mother = 0;
1337 mother = mcPartCandidate->GetMother();
1338 Int_t istep = 0;
1339 Int_t abspdgGranma =0;
1340 Bool_t isFromB=kFALSE;
1341 Bool_t isQuarkFound=kFALSE;
1342 while (mother >0 ){
1343 istep++;
1344 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1345 if (mcGranma){
1346 pdgGranma = mcGranma->GetPdgCode();
1347 abspdgGranma = TMath::Abs(pdgGranma);
1348 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1349 isFromB=kTRUE;
1350 }
1351 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1352 mother = mcGranma->GetMother();
1353 }else{
1354 AliError("Failed casting the mother particle!");
1355 break;
35151011 1356 }
1357 }
cb7c2594 1358
1359 if(isFromB) return 5;
1360 else return 4;
1361}
1362
1363
1364void AliAnalysisTaskSEDplusCorrelations::CreateCorrelationObjs() {
1365//
1366
1367 TString namePlot = "";
1368 Int_t nbinsmass=GetNBinsHistos();
1369 //These for limits in THnSparse (one per bin, same limits).
1370 //Vars: DeltaPhi, InvMass, PtTrack, Displacement, DeltaEta
1371 Int_t nBinsPhi[4] = {32,nbinsmass,30,16};
1372 Double_t binMinPhi[4] = {-TMath::Pi()/2.+TMath::Pi()/32.,fLowmasslimit,0.,-1.6}; //is the minimum for all the bins
1373 Double_t binMaxPhi[4] = {3*TMath::Pi()/2.+TMath::Pi()/32.,fUpmasslimit,30.0,1.6}; //is the maximum for all the bins
1374
1375 //Vars: DeltaPhi, InvMass, DeltaEta
1376 Int_t nBinsMix[3] = {32,nbinsmass,16};
1377 Double_t binMinMix[3] = {-TMath::Pi()/2.+TMath::Pi()/32.,fLowmasslimit,-1.6}; //is the minimum for all the bins
1378 Double_t binMaxMix[3] = {3*TMath::Pi()/2.+TMath::Pi()/32.,fUpmasslimit,1.6}; //is the maximum for all the bins
1379
1380 for(Int_t i=0;i<fNPtBins;i++){
1381
1382 if(!fMixing) {
1383 //THnSparse plots: correlations for various invariant mass (MC and data)
1384 namePlot="hPhi_Charg_Bin";
1385 namePlot+=i;
1386
1387 THnSparseF *hPhi = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1388 hPhi->Sumw2();
1389 fOutput->Add(hPhi);
1390
1391 namePlot="hPhi_Kaon_Bin";
1392 namePlot+=i;
1393
1394 THnSparseF *hPhiK = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1395 hPhiK->Sumw2();
1396 fOutput->Add(hPhiK);
1397
1398 namePlot="hPhi_K0_Bin";
1399 namePlot+=i;
1400
1401 THnSparseF *hPhiK0 = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1402 hPhiK0->Sumw2();
1403 fOutput->Add(hPhiK0);
1404
1405 //histos for c/b origin for D+ (MC only)
1406 if (fReadMC) {
1407
1408 //generic origin for tracks
1409 namePlot="hPhi_Charg_HFC_Bin";
1410 namePlot+=i;
1411
1412 THnSparseF *hPhiH_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1413 hPhiH_c->Sumw2();
1414 fOutput->Add(hPhiH_c);
1415
1416 namePlot="hPhi_Kaon_HFC_Bin";
1417 namePlot+=i;
1418
1419 THnSparseF *hPhiK_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1420 hPhiK_c->Sumw2();
1421 fOutput->Add(hPhiK_c);
1422
1423 namePlot="hPhi_K0_HFC_Bin";
1424 namePlot+=i;
1425
1426 THnSparseF *hPhiK0_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1427 hPhiK0_c->Sumw2();
1428 fOutput->Add(hPhiK0_c);
1429
1430 namePlot="hPhi_Charg_HFB_Bin";
1431 namePlot+=i;
1432
1433 THnSparseF *hPhiH_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1434 hPhiH_b->Sumw2();
1435 fOutput->Add(hPhiH_b);
1436
1437 namePlot="hPhi_Kaon_HFB_Bin";
1438 namePlot+=i;
1439
1440 THnSparseF *hPhiK_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1441 hPhiK_b->Sumw2();
1442 fOutput->Add(hPhiK_b);
1443
1444 namePlot="hPhi_K0_HFB_Bin";
1445 namePlot+=i;
1446
1447 THnSparseF *hPhiK0_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1448 hPhiK0_b->Sumw2();
1449 fOutput->Add(hPhiK0_b);
1450
1451 namePlot="hPhi_Charg_NHF_Bin";
1452 namePlot+=i;
1453
1454 THnSparseF *hPhiH_NHF = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1455 hPhiH_NHF->Sumw2();
1456 fOutput->Add(hPhiH_NHF);
1457
1458 namePlot="hPhi_Kaon_NHF_Bin";
1459 namePlot+=i;
1460
1461 THnSparseF *hPhiK_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1462 hPhiK_Non->Sumw2();
1463 fOutput->Add(hPhiK_Non);
1464
1465 namePlot="hPhi_K0_NHF_Bin";
1466 namePlot+=i;
1467
1468 THnSparseF *hPhiK0_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
1469 hPhiK0_Non->Sumw2();
1470 fOutput->Add(hPhiK0_Non);
1471 }
1472
1473 //leading hadron correlations
1474 namePlot="hPhi_Lead_Bin";
1475 namePlot+=i;
1476
1477 THnSparseF *hCorrLead = new THnSparseF(namePlot.Data(), "Leading particle correlations; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1478 hCorrLead->Sumw2();
1479 fOutput->Add(hCorrLead);
1480
1481 if (fReadMC) {
1482 namePlot="hPhi_Lead_From_c_Bin";
1483 namePlot+=i;
1484
1485 THnSparseF *hCorrLead_c = new THnSparseF(namePlot.Data(), "Leading particle correlations - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1486 hCorrLead_c->Sumw2();
1487 fOutput->Add(hCorrLead_c);
1488
1489 namePlot="hPhi_Lead_From_b_Bin";
1490 namePlot+=i;
1491
1492 THnSparseF *hCorrLead_b = new THnSparseF(namePlot.Data(), "Leading particle correlations - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1493 hCorrLead_b->Sumw2();
1494 fOutput->Add(hCorrLead_b);
1495
1496
1497
1498 namePlot="hPhi_Lead_From_nhf_Bin";
1499 namePlot+=i;
1500
1501 THnSparseF *hCorrLead_Non = new THnSparseF(namePlot.Data(), "Leading particle correlations - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1502 hCorrLead_Non->Sumw2();
1503 fOutput->Add(hCorrLead_Non);
1504 }
1505
1506
1507
1508 //pT distribution histos
1509 namePlot = "hist_Pt_Charg_Bin"; namePlot+=i;
1510 TH1F *hPtHad = new TH1F(namePlot.Data(), "Charged track pT (in D+ events); p_{T} (GeV/c)",240,0.,20.);
1511 hPtHad->SetMinimum(0);
1512 fOutput->Add(hPtHad);
1513
1514 namePlot = "hist_Pt_Kcharg_Bin"; namePlot+=i;
1515 TH1F *hPtH = new TH1F(namePlot.Data(), "Hadrons pT (in D+ events); p_{T} (GeV/c)",240,0.,12.);
1516 hPtH->SetMinimum(0);
1517 fOutput->Add(hPtH);
1518
1519 namePlot = "hist_Pt_K0_Bin"; namePlot+=i;
1520 TH1F *hPtK = new TH1F(namePlot.Data(), "Kaons pT (in D+ events); p_{T} (GeV/c)",240,0.,12.);
1521 hPtK->SetMinimum(0);
1522 fOutput->Add(hPtK);
1523
1524
35151011 1525 }
cb7c2594 1526
1527 if(fMixing) {
1528 //THnSparse plots for event mixing!
1529
1530
1531 namePlot="hPhi_K0_Bin";
1532 namePlot+=i;namePlot+="_EvMix";
1533
1534 THnSparseF *hPhiK_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1535 hPhiK_EvMix->Sumw2();
1536 fOutput->Add(hPhiK_EvMix);
1537
1538 namePlot="hPhi_Kcharg_Bin";
1539 namePlot+=i;namePlot+="_EvMix";
1540
1541 THnSparseF *hPhiH_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1542 hPhiK_EvMix->Sumw2();
1543 fOutput->Add(hPhiH_EvMix);
1544
1545 namePlot="hPhi_Charg_Bin";
1546 namePlot+=i;namePlot+="_EvMix";
1547
1548 THnSparseF *hPhiC_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1549 hPhiC_EvMix->Sumw2();
1550 fOutput->Add(hPhiC_EvMix);
35151011 1551 }
35151011 1552 }
1553
cb7c2594 1554 //out of bin loop
1555 if(!fMixing) {
1556 TH1F *hCountC = new TH1F("hist_Count_Charg", "Charged track counter; # Tracks",100,0.,100.);
1557 hCountC->SetMinimum(0);
1558 fOutput->Add(hCountC);
1559
1560 TH1F *hCountH = new TH1F("hist_Count_Kcharg", "Hadrons counter; # Tracks",20,0.,20.);
1561 hCountH->SetMinimum(0);
1562 fOutput->Add(hCountH);
1563
1564 TH1F *hCountK = new TH1F("hist_Count_K0", "Kaons counter; # Tracks",20,0.,20.);
1565 hCountK->SetMinimum(0);
1566 fOutput->Add(hCountK);
1567 }
1568
1569 if (fReadMC) {
1570 TH1D *hEventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
1571 fOutput->Add(hEventTypeMC);
1572 }
1573
1574 // if (fFillGlobal) { //all-events plots
1575 //pt distributions
35151011 1576
cb7c2594 1577 if(!fMixing) {
1578 TH1F *hPtHAll = new TH1F("hist_Pt_Charg_AllEv", "Charged track pT (All); p_{T} (GeV/c)",240,0.,12.);
1579 hPtHAll->SetMinimum(0);
1580 fOutput->Add(hPtHAll);
1581
1582 TH1F *hPtKAll = new TH1F("hist_Pt_Kaons_AllEv", "Kaons pT (All); p_{T} (GeV/c)",240,0.,12.);
1583 hPtKAll->SetMinimum(0);
1584 fOutput->Add(hPtKAll);
1585
1586 TH1F *hPtK0All = new TH1F("hist_Pt_K0_AllEv", "K0 pT (All); p_{T} (GeV/c)",240,0.,12.);
1587 hPtK0All->SetMinimum(0);
1588 fOutput->Add(hPtK0All);
1589
1590
1591 //phi distributions
1592 TH1F *hPhiDistHAll = new TH1F("hist_PhiDistr_Charg", "Charged track phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1593 hPhiDistHAll->SetMinimum(0);
1594 fOutput->Add(hPhiDistHAll);
1595
1596 TH1F *hPhiDistKAll = new TH1F("hist_PhiDistr_Kcharg", "Kaons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1597 hPhiDistKAll->SetMinimum(0);
1598 fOutput->Add(hPhiDistKAll);
1599
1600 TH1F *hPhiDistK0All = new TH1F("hist_PhiDistr_K0", "K0 phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1601 hPhiDistK0All->SetMinimum(0);
1602 fOutput->Add(hPhiDistK0All);
1603
1604 TH1F *hPhiDistDAll = new TH1F("hist_PhiDistr_Dplus", "D^{+} phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1605 hPhiDistDAll->SetMinimum(0);
1606 fOutput->Add(hPhiDistDAll);
1607 }
1608 //}
1609
1610 //for MC analysis only
1611 for(Int_t i=0;i<fNPtBins;i++) {
1612
1613 if (fReadMC && !fMixing) {
1614
1615
1616 //origin of tracks histos
1617 namePlot="histOrig_Had_Bin"; namePlot+=i;
1618 TH1F *hOrigin_had = new TH1F(namePlot.Data(), "Origin of associated tracksin MC",10,-0.5,9.5);
1619 hOrigin_had->SetMinimum(0);
1620
1621 hOrigin_had->GetXaxis()->SetBinLabel(1,"All ");
1622 hOrigin_had->GetXaxis()->SetBinLabel(2," from Heavy flavour");
1623 hOrigin_had->GetXaxis()->SetBinLabel(3," from c->D");
1624 hOrigin_had->GetXaxis()->SetBinLabel(4," from b->D");
1625 hOrigin_had->GetXaxis()->SetBinLabel(5," from b->B");
1626 hOrigin_had->GetXaxis()->SetBinLabel(6," from NHF");
1627 if(fReadMC) fOutput->Add(hOrigin_had);
1628
1629 namePlot="histOrig_Kaon_Bin"; namePlot+=i;
1630 TH1F *hOrigin_kaon = new TH1F(namePlot.Data(), "Origin of kaons in MC",10,-0.5,9.5);
1631 hOrigin_kaon->SetMinimum(0);
1632 hOrigin_kaon->GetXaxis()->SetBinLabel(1,"All Kaons");
1633 hOrigin_kaon->GetXaxis()->SetBinLabel(2,"Kaons from Heavy flavour");
1634 hOrigin_kaon->GetXaxis()->SetBinLabel(3,"Kaons from Charm");
1635 hOrigin_kaon->GetXaxis()->SetBinLabel(4,"Kaons from Beauty");
1636 hOrigin_kaon->GetXaxis()->SetBinLabel(5,"Kaons from NHF");
1637 if(fReadMC) fOutput->Add(hOrigin_kaon);
1638
1639
1640 namePlot="histOrig_K0_Bin"; namePlot+=i;
1641 TH1F *hOrigin_K0 = new TH1F(namePlot.Data(), "Origin of Kshorts in MC",10,-0.5,9.5);
1642 hOrigin_K0->SetMinimum(0);
1643
1644 hOrigin_K0->GetXaxis()->SetBinLabel(1,"All K0s");
1645 hOrigin_K0->GetXaxis()->SetBinLabel(2,"K0s from Heavy flavour");
1646 hOrigin_K0->GetXaxis()->SetBinLabel(3,"K0s from Charm");
1647 hOrigin_K0->GetXaxis()->SetBinLabel(4,"K0s from Beauty");
1648 hOrigin_K0->GetXaxis()->SetBinLabel(5,"K0s from NHF");
1649 if(fReadMC) fOutput->Add(hOrigin_K0);
1650 }
1651
1652 if (fReadMC) {
1653 //origin of Dplus histos
1654 namePlot="histOrig_Dplus_Bin"; namePlot+=i;
1655 TH1F *hOrigin_Dplus = new TH1F(namePlot.Data(),"Origin of D+ in MC",2,0.,2.);
1656 hOrigin_Dplus->SetMinimum(0);
1657 hOrigin_Dplus->GetXaxis()->SetBinLabel(1,"From c");
1658 hOrigin_Dplus->GetXaxis()->SetBinLabel(2,"From b");
1659 fOutput->Add(hOrigin_Dplus);
1660 }
1661 }
1662}
35151011 1663
1664
1665
1666
1667
1668
1669
1670
1671