]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGUD/multVScentPbPb/AliTrackletTaskMultipp.cxx
Fix typo
[u/mrichter/AliRoot.git] / PWGUD / multVScentPbPb / AliTrackletTaskMultipp.cxx
CommitLineData
c15ec07e 1/*************************************************************************
2* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
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**************************************************************************/
15
16///////////////////////////////////////////////////////////////////////////
c79a7c1c 17// Class AliTrackletTaskMultipp //
c15ec07e 18// Analysis task to produce data and MC histos needed for tracklets //
19// dNdEta extraction in multiple bins in one go //
20// Author: ruben.shahoyan@cern.ch //
21///////////////////////////////////////////////////////////////////////////
22/*
23 Important parameters to set:
24 1) make sure to initialize correct geometry in UserCreateOutputObjects
25 2) The cut on signal selection variable (delta, dphi ...) should be decided beforehand
26...
27*/
28
29#include "TChain.h"
30#include "TTree.h"
31#include "TRandom.h"
32#include "TH1F.h"
33#include "TH2F.h"
34#include "TH3F.h"
35#include "THnSparse.h"
36#include "TList.h"
37#include "TNtuple.h"
38#include "TObjArray.h"
39#include "TGeoGlobalMagField.h"
40
41#include "AliAnalysisManager.h"
42
43#include "AliMultiplicity.h"
44#include "AliESDEvent.h"
45#include "AliESDInputHandler.h"
46#include "AliESDInputHandlerRP.h"
47#include "../ANALYSIS/EventMixing/AliMixEventInputHandler.h"
48#include "AliCDBPath.h"
49#include "AliCDBManager.h"
50#include "AliCDBEntry.h"
51#include "AliCDBStorage.h"
52#include "AliGeomManager.h"
53#include "AliMagF.h"
54#include "AliESDVZERO.h"
55#include "AliESDZDC.h"
56#include "AliRunLoader.h"
57#include "AliMCEventHandler.h"
58#include "AliMCEvent.h"
59#include "AliMCParticle.h"
60#include "AliStack.h"
61#include "AliGenEventHeader.h"
62#include "AliCentrality.h"
63#include "../ITS/AliITSRecPoint.h"
64#include "../ITS/AliITSgeomTGeo.h"
65#include "../ITS/AliITSMultReconstructor.h"
c79a7c1c 66#include "../ITS/AliITSsegmentationSPD.h"
c15ec07e 67
68#include "AliLog.h"
69
70#include "AliPhysicsSelection.h"
c79a7c1c 71#include "AliTrackletTaskMultipp.h"
c15ec07e 72#include "AliITSMultRecBg.h"
73#include "AliGenEventHeader.h"
74#include "AliGenHijingEventHeader.h"
75#include "AliGenDPMjetEventHeader.h"
76#include "AliESDtrackCuts.h"
77
c79a7c1c 78ClassImp(AliTrackletTaskMultipp)
c15ec07e 79
80// centrality percentile (inverted: 100% - most central)
c79a7c1c 81const Float_t AliTrackletTaskMultipp::fgkCentPerc[] = {0,100};//{0,5,10,20,30};
82//const Float_t AliTrackletTaskMultipp::fgkCentPerc[] = {0,5,10,20,30,40};
c15ec07e 83 //{0,10,20,30,40,50,60,70,80,90,95,101};
84
c79a7c1c 85const char* AliTrackletTaskMultipp::fgkCentSelName[] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZENvsZDC"};
c15ec07e 86
c79a7c1c 87const char* AliTrackletTaskMultipp::fgkPDGNames[] = {
c15ec07e 88"#pi^{+}",
89"p",
90"K^{+}",
91"K^{*+}",
92"e^{-}",
93"#mu^{-}",
94"#rho^{+}",
95"D^{+}",
96"D^{*+}",
97"D_{s}^{+}",
98"D_{s}^{*+}",
99"#Delta^{-}",
100"#Delta^{+}",
101"#Delta^{++}",
102"#Sigma^{-}",
103"#Sigma^{+}",
104"#Sigma^{*-}",
105"#Sigma^{*+}",
106"#Sigma^{*+}_{c}",
107"#Sigma^{*++}_{c}",
108"#Xi^{-}",
109"#Xi^{*-}",
110"#Lambda^{+}_{c}",
111"n",
112"#Delta^{0}",
113"#gamma",
114"K^{0}_{S}",
115"K^{0}_{L}",
116"K^{0}",
117"K^{*}",
118"#eta",
119"#pi^{0}",
120"#rho^{0}",
121"#varphi",
122"#eta'",
123"#omega",
124"#Lambda",
125"#Sigma^{0}",
126"#Sigma^{*0}_{c}",
127"#Sigma^{*0}",
128"D^{0}",
129"D^{*0}",
130"#Xi_{0}",
131"#Xi^{*0}",
132"#Xi^{0}_{c}",
133"#Xi^{*0}_{c}",
134"Nuclei",
135"Others"
136};
137
c79a7c1c 138const int AliTrackletTaskMultipp::fgkPDGCodes[] = {
c15ec07e 139 211,
140 2212,
141 321,
142 323,
143 11,
144 13,
145 213,
146 411,
147 413,
148 431,
149 433,
150 1114,
151 2214,
152 2224,
153 3112,
154 3222,
155 3114,
156 3224,
157 4214,
158 4224,
159 3312,
160 3314,
161 4122,
162 2112,
163 2114,
164 22,
165 310,
166 130,
167 311,
168 313,
169 221,
170 111,
171 113,
172 333,
173 331,
174 223,
175 3122,
176 3212,
177 4114,
178 3214,
179 421,
180 423,
181 3322,
182 3324,
183 4132,
184 4314
185// nuclei
186// unknown
187};
188
189//________________________________________________________________________
190/*//Default constructor
c79a7c1c 191AliTrackletTaskMultipp::AliTrackletTaskMultipp(const char *name)
c15ec07e 192 : AliAnalysisTaskSE(name),
193*/
194//________________________________________________________________________
c79a7c1c 195AliTrackletTaskMultipp::AliTrackletTaskMultipp(const char *name)
c15ec07e 196 : AliAnalysisTaskSE(name),
197//
198 fOutput(0),
199//
200 fDoNormalReco(kFALSE),
201 fDoInjection(kFALSE),
202 fDoRotation(kFALSE),
203 fDoMixing(kFALSE),
204 //
205 fUseMC(kFALSE),
206 fCheckReconstructables(kFALSE),
207//
208 fHistosTrData(0),
209 fHistosTrInj(0),
210 fHistosTrRot(0),
211 fHistosTrMix(0),
212//
213 fHistosTrPrim(0),
214 fHistosTrSec(0),
215 fHistosTrComb(0),
216 fHistosTrCombU(0),
217//
218 fHistosTrRcblPrim(0),
219 fHistosTrRcblSec(0),
220 fHistosCustom(0),
221//
222 fEtaMin(-3.0),
223 fEtaMax(3.0),
224 fZVertexMin(-20),
225 fZVertexMax( 20),
226//
227 fScaleDTBySin2T(kFALSE),
228 fCutOnDThetaX(kFALSE),
229 fNStdDev(1.),
230 fDPhiWindow(0.08),
231 fDThetaWindow(0.025),
232 fDPhiShift(0.0045),
233 fPhiOverlapCut(0.005),
234 fZetaOverlap(0.05),
235 fPhiRot(0.),
236 fInjScale(1.),
237 fRemoveOverlaps(kFALSE),
238//
239 fDPhiSCut(0.06),
240 fNStdCut(1.),
241 fMCV0Scale(0.7520),
242//
243 fMultReco(0),
244 fRPTree(0),
245 fRPTreeMix(0),
246 fStack(0),
c79a7c1c 247 fMCevent(0),
c15ec07e 248 //
249 fNPart(0),
250 fNBColl(0),
251 fCurrCentBin(-1),
252 fNCentBins(0),
253 fUseCentralityVar(kCentV0M)
254{
255 // Constructor
256
257 DefineOutput(1, TList::Class());
258 //
259 SetScaleDThetaBySin2T();
260 SetNStdDev();
261 SetPhiWindow();
262 SetThetaWindow();
263 SetPhiShift();
264 SetPhiOverlapCut();
265 SetZetaOverlapCut();
266 SetPhiRot();
267 SetRemoveOverlaps();
268 //
269}
270
271//________________________________________________________________________
c79a7c1c 272AliTrackletTaskMultipp::~AliTrackletTaskMultipp()
c15ec07e 273{
274 // Destructor
275 // histograms are in the output list and deleted when the output
276 // list is deleted by the TSelector dtor
c79a7c1c 277 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { //RRR
c15ec07e 278 printf("Deleteing output\n");
279 delete fOutput;
c15ec07e 280 }
281 //
282 delete fMultReco;
283 //
284 delete fHistosTrData;
285 delete fHistosTrPrim;
286 delete fHistosTrSec;
287 delete fHistosTrComb;
288 delete fHistosTrCombU;
289 delete fHistosTrInj;
290 delete fHistosTrRot;
291 delete fHistosTrMix;
292 delete fHistosTrRcblPrim;
293 delete fHistosTrRcblSec;
294 delete fHistosCustom;
295 //
296}
297
298//________________________________________________________________________
c79a7c1c 299void AliTrackletTaskMultipp::UserCreateOutputObjects()
c15ec07e 300{
c79a7c1c 301 // create output
c15ec07e 302 fOutput = new TList();
303 fOutput->SetOwner();
304 //
305 //
306 Bool_t needGeom = GetDoNormalReco() || GetDoInjection() || GetDoRotation() || GetDoMixing();
307 if (needGeom) {
308 AliCDBManager *man = AliCDBManager::Instance();
309 if (fUseMC) {
310 Bool_t newGeom = kTRUE;
311 man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
312 if (newGeom) {
313 // new geom
314 AliCDBEntry* obj = man->Get("GRP/Geometry/Data",130844,8);
315 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
316 if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,6,-1)) AliFatal("Failed to misalign geometry");
317 }
318 else {
319 // old geom
320 AliCDBEntry* obj = man->Get("GRP/Geometry/Data",130845,7);
321 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
322 if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130845,5,-1)) AliFatal("Failed to misalign geometry");
323 }
324 }
325 else {
326 man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB"); //man->SetRun(137366);
327 AliCDBEntry* obj = man->Get("GRP/Geometry/Data",137366);
328 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
329 if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137366,-1,-1)) AliFatal("Failed to misalign geometry");
330 }
331 }
332 //
333 // Create histograms
334 fNCentBins = sizeof(fgkCentPerc)/sizeof(Float_t)-1;
335 //---------------------------------------------Standard histos per tracklet type--->>
336 UInt_t hPattern = 0xffffffff;
337 fHistosTrData = BookHistosSet("TrData",hPattern);
338 hPattern &= ~(BIT(kHEtaZvSPD1)); // fill single clusters for "data" only
339 if (GetDoInjection()) fHistosTrInj = BookHistosSet("TrInj",hPattern);
340 if (GetDoRotation()) fHistosTrRot = BookHistosSet("TrRot",hPattern);
341 if (GetDoMixing()) fHistosTrMix = BookHistosSet("TrMix",hPattern);
342 if (fUseMC) {
343 fHistosTrPrim = BookHistosSet("TrPrim",hPattern);
344 fHistosTrSec = BookHistosSet("TrSec",hPattern);
345 fHistosTrComb = BookHistosSet("TrComb",hPattern);
346 fHistosTrCombU = BookHistosSet("TrCombU",hPattern);
347 if (fCheckReconstructables) {
348 fHistosTrRcblPrim = BookHistosSet("TrRcblPrim",hPattern);
349 fHistosTrRcblSec = BookHistosSet("TrRcblSec",hPattern);
350 }
351 }
352 //---------------------------------------------Standard histos per tracklet type---<<
353 //
354 //---------------------------------------------Custom Histos----------------------->>
355 // put here any non standard histos
356 fHistosCustom = BookCustomHistos();
357 //
358 //---------------------------------------------Custom Histos-----------------------<<
359 int nhist = fOutput->GetEntries();
360 for (int i=0;i<nhist;i++) {
361 TObject* hst = fOutput->At(i);
362 if (!hst || !(hst->InheritsFrom(TH1::Class()))) continue;
363 ((TH1*)hst)->Sumw2();
364 }
365 //
366 if (fUseCentralityVar<0||fUseCentralityVar>kNCentTypes) {
367 printf("Wrong centrality type %d\n",fUseCentralityVar);
368 exit(1);
369 }
c79a7c1c 370 AliInfo(Form("Centrality type selected: %s\n",fgkCentSelName[fUseCentralityVar]));
c15ec07e 371 PostData(1, fOutput);
372 //
373}
374
375//________________________________________________________________________
c79a7c1c 376void AliTrackletTaskMultipp::UserExec(Option_t *)
c15ec07e 377{
378 // Main loop
379 //
380 Bool_t needRecPoints = GetDoNormalReco() || GetDoInjection() || GetDoRotation() || GetDoMixing();
381 //
382 AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
383 fRPTree = fRPTreeMix = 0;
384 AliESDInputHandler *handler = (AliESDInputHandler*)anMan->GetInputEventHandler();
385 AliESDInputHandlerRP *handRP = 0;
386 if (needRecPoints) {
387 handRP = (AliESDInputHandlerRP*)handler;
388 if (!handRP) { printf("No RP handler\n"); return; }
389 }
390 AliESDEvent *esd = handler->GetEvent();
391 if (!esd) { printf("No AliESDEvent\n"); return; }
392 //
393 // do we need to initialize the field?
394 AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
395 if (!field && !esd->InitMagneticField()) {printf("Failed to initialize the B field\n");return;}
396 //
397 //
398 TH1* hstat = (TH1*)fHistosCustom->UncheckedAt(kHStat);
399 hstat->Fill(kEvTot0); // RS
400 const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
401 const AliMultiplicity* multESD = esd->GetMultiplicity();
402 //
403 if (vtxESD->GetNContributors()<1) return;
404 TString vtxTyp = vtxESD->GetTitle();
405 if (vtxTyp.Contains("vertexer: Z")) {
406 if (vtxESD->GetDispersion()>0.04) return;
407 if (vtxESD->GetZRes()>0.25) return;
408 }
409 //
410 AliCentrality *centrality = esd->GetCentrality();
411 //
412 hstat->Fill(kEvTot); // RS
413 //
414 Double_t esdvtx[3];
415 vtxESD->GetXYZ(esdvtx);
416 for (int i=3;i--;) fESDVtx[i] = esdvtx[i];
417 //
418 float vtxf[3] = {vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ()};
419 //
420 // registed Ntracklets and ZVertex of the event
421 ((TH1*)fHistosCustom->UncheckedAt(kHZVtxNoSel))->Fill(esdvtx[2]);
422 //
423 if(vtxf[2] < fZVertexMin || vtxf[2] > fZVertexMax) return;
424 //
425 // centrality->Dump();
426 Float_t multV0=0;
427 AliESDVZERO* esdV0 = esd->GetVZEROData();
428 if (esdV0) {
429 multV0 = esdV0->GetMTotV0A()+esdV0->GetMTotV0C();
430 if (fUseMC) multV0 *= fMCV0Scale;
431 }
432 ((TH1*)fHistosCustom->UncheckedAt(kHV0NoSel))->Fill(multV0);
433 //
434 float nSPD2 = multESD->GetNumberOfITSClusters(1);
435 ((TH1*)fHistosCustom->UncheckedAt(kHNClSPD2NoSel))->Fill(nSPD2);
436 //
437 //------------------------------------------------------
438 AliESDZDC *esdZDC = esd->GetESDZDC();
439 float zdcEnergy=0,zemEnergy=0;
440 if (esdZDC) {
441 zdcEnergy = (esdZDC->GetZDCN1Energy() + esdZDC->GetZDCP1Energy() + esdZDC->GetZDCN2Energy()+ esdZDC->GetZDCP2Energy());
442 zemEnergy = (esdZDC->GetZDCEMEnergy(0)+esdZDC->GetZDCEMEnergy(1))/8.;
443 }
444 ((TH2*)fHistosCustom->UncheckedAt(kHZDCZEMNoSel))->Fill(zemEnergy,zdcEnergy);
445 //
c79a7c1c 446 Float_t centPercentile = centrality->GetCentralityPercentileUnchecked(fgkCentSelName[fUseCentralityVar]);
c15ec07e 447
448 // temporary >>>>>>>>>>>>>>>>>>>>>>>>
449 if (fUseCentralityVar==kCentZEMvsZDC) {
450 float zdcEn = zdcEnergy;
451 float zemEn = zemEnergy;
452 Float_t slope;
453 Float_t zdcPercentile;
454 if (zemEn > 295.) {
455 slope = (zdcEn + 15000.)/(zemEn - 295.);
456 slope += 2.23660e+02;
457 zdcPercentile = (TMath::ATan(slope) - 1.56664)/8.99571e-05;
458 if (zdcPercentile<0) zdcPercentile = 0;
459 }
460 else zdcPercentile = 100;
461 centPercentile = zdcPercentile;
462 }
463 // temporary >>>>>>>>>>>>>>>>>>>>>>>>
464
465 fCurrCentBin = GetCentralityBin(centPercentile);
466 //
467 // printf("CentPerc: %f : Bin %d\n",centPercentile, fCurrCentBin);
468 if (fCurrCentBin<0) {
469 //printf("Reject: %.1f : V0:%.1f V0Cor:%.1f V0CR:%.1f SPD2c:%.1f\n",mltTst, multV0,multV0Corr,multV0CorrResc,nSPD2Corr);
470 return;
471 }
472 //
473 ((TH1*)fHistosCustom->UncheckedAt(kHStatCentBin))->Fill(fCurrCentBin);
474 ((TH1*)fHistosCustom->UncheckedAt(kHStatCent))->Fill(centPercentile);
475 //
476 AliMCEventHandler* eventHandler = 0;
c79a7c1c 477 fMCevent = 0;
c15ec07e 478 fStack = 0;
479 //
480 if (fUseMC) {
481 eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
482 if (!eventHandler) { printf("ERROR: Could not retrieve MC event handler\n"); return; }
c79a7c1c 483 fMCevent = eventHandler->MCEvent();
484 if (!fMCevent) { printf("ERROR: Could not retrieve MC event\n"); return; }
485 fStack = fMCevent->Stack();
c15ec07e 486 if (!fStack) { printf("Stack not available\n"); return; }
487 }
488 //
489 if (needRecPoints) {
490 fRPTree = handRP->GetTreeR("ITS");
491 if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
492 }
493 //
494 // =============================================================================>>>
495 // MC Generator info
496 AliGenEventHeader* mcGenH = 0;
497 fNPart = 0;
498 fNBColl = 0;
499 if (fUseMC) {
c79a7c1c 500 mcGenH = fMCevent->GenEventHeader();
c15ec07e 501 if (mcGenH->InheritsFrom(AliGenHijingEventHeader::Class())) {
502 AliGenHijingEventHeader* hHijing = (AliGenHijingEventHeader*)mcGenH;
503 fNPart = (hHijing->ProjectileParticipants()+hHijing->TargetParticipants())/2.;
504 fNBColl = hHijing->NN()+hHijing->NNw()+hHijing->NwN()+hHijing->NwNw();
505 }
506 else if (mcGenH->InheritsFrom(AliGenDPMjetEventHeader::Class())) {
507 AliGenDPMjetEventHeader* hDpmJet = (AliGenDPMjetEventHeader*)mcGenH;
508 fNPart = (hDpmJet->ProjectileParticipants()+hDpmJet->TargetParticipants())/2.;
509 fNBColl = hDpmJet->NN()+hDpmJet->NNw()+hDpmJet->NwN()+hDpmJet->NwNw();
510 }
511 else {} // unknown generator
512 }
513 //
514 // register Ntracklets and ZVertex of the event
515 ((TH2*)fHistosCustom->UncheckedAt(kHZVtx))->Fill(esdvtx[2],fCurrCentBin);
516 ((TH2*)fHistosCustom->UncheckedAt(kHV0))->Fill(multV0,fCurrCentBin);
517 ((TH2*)fHistosCustom->UncheckedAt(kHNClSPD2))->Fill(nSPD2,fCurrCentBin);
518 ((TH3*)fHistosCustom->UncheckedAt(kHZDCZEM))->Fill(zemEnergy,zdcEnergy,fCurrCentBin);
519 //
520 if (fUseMC) FillMCPrimaries();
521 //
522 // normal reconstruction
523 hstat->Fill(kBinEntries+kEvProcData + kEntriesPerBin*fCurrCentBin);
524 //
525 if (GetDoNormalReco() || GetDoInjection()) { // for the injection the normal reco should be done
526 InitMultReco();
527 fMultReco->Run(fRPTree, vtxf);
528 printf("Multiplicity Reconstructed:\n");
529 AliMultiplicity* mlt = fMultReco->GetMultiplicity();
530 if (mlt) mlt->Print();
531 if (GetDoNormalReco()) FillHistos(kData,mlt);
532 FillClusterInfo();
533 //
534 }
535 if (!GetDoNormalReco()) {
536 FillHistos(kData,multESD); // fill data histos from ESD
537 FillClusterInfoFromMult(multESD, vtxf[2] );
538 FillClusterAutoCorrelationFromMult(multESD, vtxf[2]);
539 }
540 //
541 // Injection: it must come right after the normal reco since needs its results
542 if (GetDoInjection()) {
543 if (!fMultReco) InitMultReco(); // in principle, not needed, the reco is created above
544 fMultReco->SetRecType(AliITSMultRecBg::kBgInj);
545 fMultReco->Run(fRPTree, vtxf);
546 printf("Multiplicity from Injection:\n");
547 AliMultiplicity* mlt = fMultReco->GetMultiplicity();
548 if (mlt) mlt->Print();
549 hstat->Fill(kBinEntries + kEvProcInj + kEntriesPerBin*fCurrCentBin);
550 FillHistos(kBgInj,mlt);
551 }
552 //
553 // Rotation
554 if (GetDoRotation()) {
555 InitMultReco();
556 fMultReco->SetRecType(AliITSMultRecBg::kBgRot);
557 fMultReco->SetPhiRotationAngle(fPhiRot);
558 fMultReco->Run(fRPTree, vtxf);
559 printf("Multiplicity from Rotation:\n");
560 AliMultiplicity* mlt = fMultReco->GetMultiplicity();
561 if (mlt) mlt->Print();
562 hstat->Fill(kBinEntries + kEvProcRot + kEntriesPerBin*fCurrCentBin);
563 FillHistos(kBgRot,mlt);
564 }
565 //
566 if (GetDoMixing()) {
567 /*
568 AliMixEventInputHandler* handToMix = (AliMixEventInputHandler*)handRP->MixingHandler();
569 if (!handToMix) { printf("No Mixing handler\n"); return; }
570 handToMix->GetEntry();
571 if(handToMix->MixedEventNumber()<1) {printf("Mixing: No enough events in pool\n"); return;}
572 AliESDInputHandlerRP* handRPMix = (AliESDInputHandlerRP*) handToMix->InputEventHandler(0);
573
574 if (!handRPMix) { printf("No Mixing RP handler\n"); return; }
575 fRPTreeMix = handRPMix->GetTreeR("ITS");
576 if (!fRPTreeMix) { AliError(" Invalid ITS cluster tree of the 2nd event!\n"); return; }
577 //
578 AliESDEvent *esdMix = handRPMix->GetEvent();
579 const AliESDVertex* vtxESDmix = esdMix->GetVertex();
580 ((TH2*)fHistosCustom->UncheckedAt(kHZVtxMixDiff))->Fill(vtxESDmix->GetZ()-esdvtx[2],fCurrCentBin);
581 ((TH2*)fHistosCustom->UncheckedAt(kHNTrMixDiff) )->
582 Fill(esdMix->GetMultiplicity()->GetNumberOfTracklets() - multESD->GetNumberOfTracklets(),fCurrCentBin);
583 //
584 InitMultReco();
585 fMultReco->SetRecType(AliITSMultRecBg::kBgMix);
586 fMultReco->Run(fRPTree, vtxf,fRPTreeMix);
587 printf("Multiplicity from Mixing:\n");
588 AliMultiplicity* mlt = fMultReco->GetMultiplicity();
589 if (mlt) mlt->Print();
590 hstat->Fill(kBinEntries + kEvProcMix + kEntriesPerBin*fCurrCentBin);
591 FillHistos(kBgMix,mlt);
592 */
593 //
594 }
595 // =============================================================================<<<
596 //
597 if (fMultReco) delete fMultReco;
598 fMultReco = 0;
599 //
600}
601
602
603//________________________________________________________________________
c79a7c1c 604void AliTrackletTaskMultipp::Terminate(Option_t *)
c15ec07e 605{
c79a7c1c 606 // finish processing
c15ec07e 607 Printf("Terminating...");
608 TH1* hstat;
609 TList *lst = dynamic_cast<TList*>(GetOutputData(1));
610 printf("Term: %p %p %p\n",fOutput,lst,fHistosCustom);
611 if (lst && (hstat=(TH1*)lst->FindObject("hStat"))) {
612 Info("Terminate","Registering used settings");
613 // fill used settings
614 hstat->Fill(kOneUnit,1.);
615 hstat->Fill(kCentVar,fUseCentralityVar);
616 hstat->Fill(kDPhi,fDPhiWindow);
617 hstat->Fill(kDTht,fDThetaWindow);
618 hstat->Fill(kNStd,fNStdDev);
619 hstat->Fill(kPhiShift,fDPhiShift);
620 hstat->Fill(kThtS2,fScaleDTBySin2T);
621 hstat->Fill(kThtCW,fCutOnDThetaX);
622 hstat->Fill(kPhiOvl,fPhiOverlapCut);
623 hstat->Fill(kZEtaOvl,fZetaOverlap);
624 hstat->Fill(kNoOvl,fRemoveOverlaps);
625 //
626 hstat->Fill(kPhiRot,fPhiRot);
627 hstat->Fill(kInjScl,fInjScale);
628 hstat->Fill(kEtaMin,fEtaMin);
629 hstat->Fill(kEtaMax,fEtaMax);
630 hstat->Fill(kZVMin,fZVertexMin);
631 hstat->Fill(kZVMax,fZVertexMax);
632 //
633 hstat->Fill(kDPiSCut,fDPhiSCut);
634 hstat->Fill(kNStdCut,fNStdCut);
635 hstat->Fill(kMCV0Scale, fMCV0Scale);
636 //
637 }
638 //
639 // AliAnalysisTaskSE::Terminate();
640}
641
642
643//_________________________________________________________________________
c79a7c1c 644void AliTrackletTaskMultipp::InitMultReco()
c15ec07e 645{
646 // create mult reconstructor
647 if (fMultReco) delete fMultReco;
648 fMultReco = new AliITSMultRecBg();
649 fMultReco->SetCreateClustersCopy(kTRUE);
650 fMultReco->SetScaleDThetaBySin2T(fScaleDTBySin2T);
651 fMultReco->SetNStdDev(fNStdDev);
652 fMultReco->SetPhiWindow( fDPhiWindow );
653 fMultReco->SetThetaWindow( fDThetaWindow );
654 fMultReco->SetPhiShift( fDPhiShift );
655 fMultReco->SetRemoveClustersFromOverlaps(fRemoveOverlaps);
656 fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
657 fMultReco->SetZetaOverlapCut(fZetaOverlap);
658 fMultReco->SetHistOn(kFALSE);
659 fMultReco->SetRecType( AliITSMultRecBg::kData );
660}
661
662//_________________________________________________________________________
c79a7c1c 663TObjArray* AliTrackletTaskMultipp::BookCustomHistos()
c15ec07e 664{
665 // book custom histos, not related to specific tracklet type
666 TObjArray* histos = new TObjArray();
667 TH1F* hstat;
668 //
669 // ------------ job parameters, statistics ------------------------------>>>
670 int nbs = kBinEntries + fNCentBins*kEntriesPerBin;
671 hstat = new TH1F("hStat","Run statistics",nbs,0.5,nbs+0.5);
672 //
673 hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot0");
674 hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot");
675 hstat->GetXaxis()->SetBinLabel(kOneUnit,"ScaleMerge");
676 hstat->GetXaxis()->SetBinLabel(kNWorkers,"Workers");
677 //
c79a7c1c 678 hstat->GetXaxis()->SetBinLabel(kCentVar, fgkCentSelName[fUseCentralityVar]);
c15ec07e 679 hstat->GetXaxis()->SetBinLabel(kDPhi, "#Delta#varphi");
680 hstat->GetXaxis()->SetBinLabel(kDTht, "#Delta#theta");
681 hstat->GetXaxis()->SetBinLabel(kNStd, "N.std");
682 hstat->GetXaxis()->SetBinLabel(kPhiShift,"#delta#varphi");
683 hstat->GetXaxis()->SetBinLabel(kThtS2,"scale #Delta#theta");
684 hstat->GetXaxis()->SetBinLabel(kPhiOvl,"#varpho_{Ovl}");
685 hstat->GetXaxis()->SetBinLabel(kZEtaOvl,"#z_{Ovl}");
686 hstat->GetXaxis()->SetBinLabel(kNoOvl, "rem.ovl");
687 //
688 hstat->GetXaxis()->SetBinLabel(kPhiRot,"#varphi_{rot}");
689 hstat->GetXaxis()->SetBinLabel(kInjScl,"inj");
690 hstat->GetXaxis()->SetBinLabel(kEtaMin,"#eta_{min}");
691 hstat->GetXaxis()->SetBinLabel(kEtaMax,"#eta_{max}");
692 hstat->GetXaxis()->SetBinLabel(kZVMin,"ZV_{min} cut");
693 hstat->GetXaxis()->SetBinLabel(kZVMax,"ZV_{max} cut");
694 //
695 hstat->GetXaxis()->SetBinLabel(kDPiSCut,"#Delta#varphi-#delta_{#phi} cut");
696 hstat->GetXaxis()->SetBinLabel(kNStdCut,"#Delta cut");
697 //
698 hstat->GetXaxis()->SetBinLabel(kMCV0Scale,"MC V0 scale");
699 //
700 for (int i=0;i<fNCentBins;i++) {
701 TString bnt = "b"; bnt+= i;
702 int offs = kBinEntries + i*kEntriesPerBin;
703 hstat->GetXaxis()->SetBinLabel(offs + kEvProcData, bnt+" Ev.ProcData");
704 hstat->GetXaxis()->SetBinLabel(offs + kEvProcInj, bnt+" Ev.ProcInj");
705 hstat->GetXaxis()->SetBinLabel(offs + kEvProcRot, bnt+" Ev.ProcRot");
706 hstat->GetXaxis()->SetBinLabel(offs + kEvProcMix, bnt+" Ev.ProcMix");
707 //
708 }
709 //
710 hstat->Fill(kNWorkers);
711 //
712 AddHisto(histos,hstat,kHStat);
713 //
714 // ------------------------ events per centrality bin ----------------------
715 TH1D* hCentAx = new TH1D("EvCentr","Events per centrality",fNCentBins,fgkCentPerc);
716 hCentAx->GetXaxis()->SetTitle("Centrality parameter");
717 AddHisto(histos,hCentAx,kHStatCent);
718 //
719 TH1D* hCentBin = new TH1D("EvCentrBin","Events per centrality bin",fNCentBins,-0.5,fNCentBins-0.5);
720 hCentBin->GetXaxis()->SetTitle("Centrality Bin");
721 AddHisto(histos,hCentBin,kHStatCentBin);
722 //
723 // ------------ job parameters, statistics ------------------------------<<<
724 //
725 double etaMn=-3,etaMx=3;
726 double zMn=-30, zMx=30;
727 int nEtaBins = int((etaMx-etaMn)/0.1);
728 if (nEtaBins<1) nEtaBins = 1;
729 //
730 int nZVBins = int(zMx-zMn);
731 if (nZVBins<1) nZVBins = 1;
732 //
733 // Z vertex distribution for events before selection
734 TH1F* hzvns = new TH1F("zvNoSel","Z vertex before selection",nZVBins,zMn,zMx);
735 hzvns->GetXaxis()->SetTitle("Zvertex");
736 AddHisto(histos,hzvns,kHZVtxNoSel);
737 //
738 // V0 for events before selection
739 int nbmltV0 = 250;
740 double maxmltV0 = 25000;
741 //
742 TH1F* hnV0ns = new TH1F("V0NoSel","V0 signal Before Cent. Selection",nbmltV0,0,maxmltV0);
743 hnV0ns->GetXaxis()->SetTitle("V0 signal");
744 AddHisto(histos,hnV0ns,kHV0NoSel);
745 //
746 // N SPD2 clusters
747 int nbmltSPD2 = 175;
748 double maxmltSPD2 = 7000;
749 TH1F* hncl2ns = new TH1F("NClustersSPD2NoSel","N Clusters on SPD2 Before Cent Selection",nbmltSPD2,0,maxmltSPD2);
750 hncl2ns->GetXaxis()->SetTitle("N Clus SPD2");
751 AddHisto(histos,hncl2ns,kHNClSPD2NoSel);
752 //
753 int nbzdc=50;
754 double maxZDC=6000., maxZEM=2500.;
755 TH2F* hzdczemns = new TH2F("ZDCZEMNoSel","ZDC vs ZEM Before Cent Selection",
756 nbzdc,0,maxZEM,nbzdc,0,maxZDC);
757 hzdczemns->GetXaxis()->SetTitle("ZEM");
758 hzdczemns->GetXaxis()->SetTitle("ZDC");
759 AddHisto(histos,hzdczemns,kHZDCZEMNoSel);
760 //
761 TH2F* hzv = new TH2F("zv","Z vertex after Selection per Cent.Bin",nZVBins,zMn,zMx, fNCentBins, -0.5,fNCentBins-0.5);
762 hzv->GetXaxis()->SetTitle("Zvertex");
763 hzv->GetYaxis()->SetTitle("Cent.Bin ID");
764 AddHisto(histos,hzv,kHZVtx);
765 //
766 // V0
767 TH2F* hnV0 = new TH2F("V0","V0 signal per Cent.Bin ",nbmltV0,0,maxmltV0, fNCentBins, -0.5,fNCentBins-0.5);
768 hnV0->GetXaxis()->SetTitle("V0 signal");
769 hnV0->GetYaxis()->SetTitle("Cent.Bin ID");
770 AddHisto(histos,hnV0,kHV0);
771 //
772 // N SPD2 clusters
773 TH2F* hncl2 = new TH2F("NClustersSPD2","N Clusters on SPD2 per Cent.Bin",nbmltSPD2,0,maxmltSPD2, fNCentBins, -0.5,fNCentBins-0.5);
774 hncl2->GetXaxis()->SetTitle("N Clus SPD2");
775 hncl2->GetYaxis()->SetTitle("Cent.Bin ID");
776 AddHisto(histos,hncl2,kHNClSPD2);
777 //
778 // ZDCZEM
779 TH3F* hzdczem = new TH3F("ZDCZEM","ZDC vs ZEM per Cent.Bin",nbzdc,0,maxZEM,nbzdc,0,maxZDC, fNCentBins, -0.5,fNCentBins-0.5);
780 hzdczem->GetXaxis()->SetTitle("ZEM");
781 hzdczem->GetYaxis()->SetTitle("ZDC");
782 AddHisto(histos,hzdczem,kHZDCZEM);
783 //
784 //----------------------------------------------------------------------
785 int nEtaBinsS = int((fEtaMax-fEtaMin)/0.1);
786 if (nEtaBinsS<1) nEtaBins = 1;
787 //
788 int nZVBinsS = int(fZVertexMax-fZVertexMin);
789 if (nZVBinsS<1) nZVBinsS = 1;
790
791 if (fUseMC) {
792 // Z vertex vs Eta distribution for primaries
793 char buffn[100],bufft[500];
794 for (int ib=0;ib<fNCentBins;ib++) {
795 sprintf(buffn,"b%d_zvEtaPrimMC",ib);
796 sprintf(bufft,"bin%d Zvertex vs #eta PrimMC",ib);
797 TH2F* hzvetap = new TH2F(buffn,bufft, nEtaBinsS,fEtaMin,fEtaMax,nZVBinsS,fZVertexMin,fZVertexMax);
798 hzvetap->GetXaxis()->SetTitle("#eta");
799 hzvetap->GetYaxis()->SetTitle("Zvertex");
800 AddHisto(histos,hzvetap,kHZVEtaPrimMC+ib);
801 }
802 //
803 // <n> primaries according to MC generator
804 TH1F* hnprimM = new TH1F("nPrimMean","<N> primaries",fNCentBins, -0.5,fNCentBins-0.5);
805 hnprimM->GetXaxis()->SetTitle("Cent.Bin ID");
806 AddHisto(histos,hnprimM,kHNPrimMeanMC);
807 //
808 // <n> primaries per part.pair according to MC generator
809 TH1F* hnprim2partM = new TH1F("nPrim2Part","<N> primaries per part.pair",fNCentBins, -0.5,fNCentBins-0.5);
810 hnprim2partM->GetXaxis()->SetTitle("Cent.Bin ID");
811 AddHisto(histos,hnprim2partM,kHNPrim2PartMC);
812 //
813 // <n> primaries per part.pair vs npart.pair according to MC generator
814 TH2F* hnprim2partNp = new TH2F("nPrim2Part_vs_NPart","<N> primaries per part.pair vs N part.pairs",105,0,210,200,0,40);
815 hnprim2partNp->GetXaxis()->SetTitle("N.part.pairs");
816 hnprim2partNp->GetYaxis()->SetTitle("N.prim/N.part.pairs");
817 AddHisto(histos,hnprim2partNp,kHNPrim2PartNpMC);
818 //
819 // <n> primaries per b.coll vs npart.pair according to MC generator
820 TH2F* hnprim2BCollNp = new TH2F("nPrim2BColl_vs_NPart","<N> primaries per bin.coll vs N part.pairs",105,0,210,200,0,40);
821 hnprim2BCollNp->GetXaxis()->SetTitle("N.part.pairs");
822 hnprim2BCollNp->GetYaxis()->SetTitle("N.prim/N.bin.coll.");
823 AddHisto(histos,hnprim2BCollNp,kHNPrim2BCollNpMC);
824 //
825 // <n> primaries per bin.coll. according to MC generator
826 TH1F* hnprim2BCollM = new TH1F("nPrim2BColl","<N> primaries per bin.coll",fNCentBins, -0.5,fNCentBins-0.5);
827 hnprim2BCollM->GetXaxis()->SetTitle("Cent.Bin ID");
828 AddHisto(histos,hnprim2BCollM,kHNPrim2BCollMC);
829 //
830 // n participants according to MC generator
831 TH2F* hnpart = new TH2F("nPart","N participant pairs",210,0,210,fNCentBins, -0.5,fNCentBins-0.5);
832 hnpart->GetXaxis()->SetTitle("N part. pairs");
833 hnpart->GetYaxis()->SetTitle("Cent.Bin ID");
834 AddHisto(histos,hnpart,kHNPartMC);
835 //
836 // <n> participants according to MC generator
837 TH1F* hnpartM = new TH1F("nPartMean","<N> participant pairs",fNCentBins, -0.5,fNCentBins-0.5);
838 hnpartM->GetXaxis()->SetTitle("Cent.Bin ID");
839 AddHisto(histos,hnpartM,kHNPartMeanMC);
840 //
841 // n bin coll. according to MC generator
842 TH2F* hnbcoll = new TH2F("nBColl","N bin. coll",2000,0,2000,fNCentBins, -0.5,fNCentBins-0.5);
843 hnbcoll->GetXaxis()->SetTitle("N bin. coll");
844 hnbcoll->GetYaxis()->SetTitle("Cent.Bin ID");
845 AddHisto(histos,hnbcoll,kHNBCollMC);
846 //
847 // <n> bin col according to MC generator
848 TH1F* hnbcollM = new TH1F("nBCollMean","<N> bin.colls",fNCentBins, -0.5,fNCentBins-0.5);
849 hnbcollM->GetXaxis()->SetTitle("Cent.Bin ID");
850 AddHisto(histos,hnbcollM,kHNBCollMeanMC);
851 //
852 }
853 //
854 if (GetDoMixing()) {
855 //
856 // Difference in Z vertex for mixed events
857 TH2F* hzdiff = new TH2F("MixSPDVertexDiff","SPD #Delta Z Vertex distribution per mult bin ",100,-5,5, fNCentBins, -0.5,fNCentBins-0.5);
858 hzdiff->GetXaxis()->SetTitle("#Delta Z Vertex [cm]");
859 hzdiff->GetYaxis()->SetTitle(Form("Entries / %1.2f [cm] per mult bin",10./100.));
860 AddHisto(histos,hzdiff,kHZVtxMixDiff);
861 //
862 // Difference in N tracklets for mixed events
863 TH2F* hntdiff = new TH2F("MixNTrackletsDiff"," SPD tracklets Diff ",200,-1000,1000, fNCentBins, -0.5,fNCentBins-0.5);
864 hntdiff->GetXaxis()->SetTitle("# tracklet diff");
865 AddHisto(histos,hntdiff,kHNTrMixDiff);
866 }
867 //
868 // --------------------------------------------------
869 if (fUseMC) {
870 int npdg = sizeof(fgkPDGNames)/sizeof(char*);
871 TH2F* hpdgP = new TH2F("pdgPrim","primary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
872 AddHisto(histos,hpdgP,kHPrimPDG);
873 TH2F* hpdgS = new TH2F("pdgSec","secondary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
874 AddHisto(histos,hpdgS,kHSecPDG);
875 TH2F* hpdgPP = new TH2F("pdgPrimPar","primary parent PDG ",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
876 AddHisto(histos,hpdgPP,kHPrimParPDG);
877 TH2F* hpdgSP = new TH2F("pdgSecPar","secondary parent PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
878 AddHisto(histos,hpdgSP,kHSecParPDG);
879 for (int i=0;i<npdg;i++) {
880 hpdgP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
881 hpdgS->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
882 hpdgPP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
883 hpdgSP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
884 }
885 }
886 //
887 // -------------------------------------------------
888 TH2F* hclinf=0;
889 hclinf = new TH2F("cl0InfoUsed","#phi vs Z of used clusters, Lr0",60,-15,15, 80,0,2*TMath::Pi());
890 AddHisto(histos,hclinf,kHClUsedInfoL0);
891 hclinf = new TH2F("cl1InfoUsed","#phi vs Z of used clusters, Lr1",60,-15,15, 2*80,0,2*TMath::Pi());
892 AddHisto(histos,hclinf,kHClUsedInfoL1);
893 hclinf = new TH2F("cl0InfoAll","#phi vs Z of all clusters, Lr0",60,-15,15, 80,0,2*TMath::Pi());
894 AddHisto(histos,hclinf,kHClAllInfoL0);
895 hclinf = new TH2F("cl1InfoAll","#phi vs Z of all clusters, Lr1",60,-15,15, 2*80,0,2*TMath::Pi());
896 AddHisto(histos,hclinf,kHClAllInfoL1);
897 //
898 // -------------------------------------------------
899 // correlations between SPD1 clusters
900 hclinf = new TH2F("clDstVsZall","Distance between SPD1 clusters vs Z, all cl",30,-15,15, 50,0.,5.);
901 AddHisto(histos,hclinf,kHclDstZAll);
902 hclinf = new TH2F("clDstVsPhiall","Distance between SPD1 clusters vs #phi, all cl",30,0,2*TMath::Pi(),50,0.,5.);
903 AddHisto(histos,hclinf,kHclDstPhiAll);
904 //
905 hclinf = new TH2F("clDstVsZused","Distance between SPD1 clusters vs Z, used-unused cl",30,-15,15, 50,0.,5.);
906 AddHisto(histos,hclinf,kHclDstZUsed);
907 hclinf = new TH2F("clDstVsPhiused","Distance between SPD1 clusters vs #phi, used-unused cl",30,0,2*TMath::Pi(),50,0.,5.);
908 AddHisto(histos,hclinf,kHclDstPhiUsed);
909 //
910 histos->SetOwner(kFALSE);
911 //
912 return histos;
913}
914
915//_________________________________________________________________________
c79a7c1c 916TObjArray* AliTrackletTaskMultipp::BookHistosSet(const char* pref, UInt_t selHistos)
c15ec07e 917{
918 // book standard set of histos attaching the pref in front of the name/title
919 //
920 const int kNDPhiBins = 100;
921 const int kNDThtBins = 100;
922 int nDistBins = int(fNStdDev)*5;
923 //
924 int nEtaBins = int((fEtaMax-fEtaMin)/0.1);
925 if (nEtaBins<1) nEtaBins = 1;
926 //
927 int nZVBins = int(fZVertexMax-fZVertexMin);
928 if (nZVBins<1) nZVBins = 1;
929 float dphir = fDPhiWindow*TMath::Sqrt(fNStdDev);
930 float dthtr = fDThetaWindow*TMath::Sqrt(fNStdDev);
931 //
932 TObjArray* histos = new TObjArray();
933 TH2F* h2;
934 TH1F* h1;
935 char buffn[100],bufft[500];
936 //
937 for (int ib=0;ib<fNCentBins;ib++) {
938 //
939 int offs = ib*kNStandardH;
940 if (selHistos & (0x1<<kHEtaZvCut) ) {
941 sprintf(buffn,"b%d_%s_ZvEtaCutT",ib,pref);
942 sprintf(bufft,"bin%d (%s) Zv vs Eta with tracklet cut",ib,pref);
943 h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
944 h2->GetXaxis()->SetTitle("#eta");
945 h2->GetYaxis()->SetTitle("Zv");
946 AddHisto(histos,h2,offs+kHEtaZvCut);
947 }
948 //
949 if (selHistos & (0x1<<kHDPhiDTheta) ) {
950 sprintf(buffn,"b%d_%s_dPhidTheta",ib,pref);
951 sprintf(bufft,"bin%d (%s) #Delta#theta vs #Delta#varphi",ib,pref);
952 h2 = new TH2F(buffn,bufft,kNDPhiBins,-dphir,dphir,kNDThtBins,-dthtr,dthtr);
953 h2->GetXaxis()->SetTitle("#Delta#varphi [rad]");
954 h2->GetYaxis()->SetTitle("#Delta#theta [rad]");
955 AddHisto(histos,h2,offs+kHDPhiDTheta);
956 }
957 //
958 if (selHistos & (0x1<<kHDPhiSDThetaX) ) {
959 sprintf(buffn,"b%d_%s_dPhiSdThetaX",ib,pref);
960 sprintf(bufft,"bin%d (%s) #Delta#theta%s vs #Delta#varphi-#delta_{#varphi}",ib,pref,fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
961 h2 = new TH2F(buffn,bufft,kNDPhiBins,-dphir,dphir,kNDThtBins,-dthtr,dthtr);
962 h2->GetXaxis()->SetTitle("#Delta#varphi-#delta_{#varphi} [rad]");
963 sprintf(bufft,"#Delta#theta%s",fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
964 h2->GetYaxis()->SetTitle(bufft);
965 AddHisto(histos,h2,offs+kHDPhiSDThetaX);
966 }
967 //
968 if (selHistos & (0x1<<kHWDist) ) {
969 sprintf(buffn,"b%d_%s_WDist",ib,pref);
970 sprintf(bufft,"bin%d #Delta=[(#Delta#varphi-#delta_{#varphi})/#sigma#varphi]^{2}+"
971 "[#Delta#theta%s/#sigma#theta]^{2}",ib,fScaleDTBySin2T ? "*sin^{-2}(#theta)":"");
972 h1 = new TH1F(buffn,bufft,nDistBins,0,fNStdDev);
973 sprintf(bufft,"#Delta=[(#Delta#varphi-#delta_{#varphi})/#sigma#varphi]^{2}+"
974 "[#Delta#theta%s/#sigma#theta]^{2}",fScaleDTBySin2T ? "*sin^{-2}(#theta)":"");
975 h1->GetXaxis()->SetTitle(bufft);
976 AddHisto(histos,h1,offs+kHWDist);
977 }
978 //
979 if (selHistos & (0x1<<kHEtaZvSPD1) ) {
980 sprintf(buffn,"b%d_%s_ZvEtaSPD1",ib,pref);
981 sprintf(bufft,"bin%d (%s) Zv vs Eta SPD1 clusters",ib,pref);
982 h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
983 h2->GetXaxis()->SetTitle("#eta");
984 h2->GetYaxis()->SetTitle("Zv");
985 AddHisto(histos,h2,offs+kHEtaZvSPD1);
986 }
987 //
988 }
989 //
990 histos->SetOwner(kFALSE);
991 return histos;
992}
993
994//_________________________________________________________________________
c79a7c1c 995void AliTrackletTaskMultipp::AddHisto(TObjArray* histos, TObject* h, Int_t at)
c15ec07e 996{
997 // add single histo to the set
998 if (at>=0) histos->AddAtAndExpand(h,at);
999 else histos->Add(h);
1000 fOutput->Add(h);
1001}
1002
1003//_________________________________________________________________________
c79a7c1c 1004void AliTrackletTaskMultipp::FillHistos(Int_t type, const AliMultiplicity* mlt)
c15ec07e 1005{
1006 // fill histos of given type
1007 if (!mlt) return;
1008 //
1009 TObjArray* histos = 0;
1010 if (type == kData) histos = fHistosTrData;
1011 else if (type == kBgInj) histos = fHistosTrInj;
1012 else if (type == kBgRot) histos = fHistosTrRot;
1013 else if (type == kBgMix) histos = fHistosTrMix;
1014 //
1015 Bool_t fillMC = (type==kData) && fUseMC && fStack;
1016 //
1017 //
1018 //---------------------------------------- CHECK ------------------------------>>>
1019 TArrayF vtxMC;
1020 AliGenHijingEventHeader* pyHeader = 0;
1021 //
1022 if (fUseMC) {
c79a7c1c 1023 pyHeader = (AliGenHijingEventHeader*) fMCevent->GenEventHeader();//header->GenEventHeader();
c15ec07e 1024 pyHeader->PrimaryVertex(vtxMC);
1025 }
1026 //---------------------------------------- CHECK ------------------------------<<<
1027 //
1028 if (!histos) return;
1029 int ntr = mlt->GetNumberOfTracklets();
1030 for (int itr=ntr;itr--;) {
1031 //
1032 //---------------------------------------- CHECK ------------------------------>>>
1033 /*
1034 if (fUseMC) {
1035 Bool_t reject = kFALSE;
1036 while(1) {
1037 int lab0 = mlt->GetLabel(itr,0);
1038 int lab1 = mlt->GetLabel(itr,1);
1039 if (lab0!=lab1) break;
1040 if (!fStack->IsPhysicalPrimary(lab0)) break;
1041 //
1042 TParticle* part = fStack->Particle(lab0);
1043 Float_t dz = part->Vz() - vtxMC[2];
1044 if (TMath::Abs(dz)<1e-6) break;
1045 reject = kTRUE;
1046 break;
1047 }
1048 if (reject) continue;
1049 }
1050 */
1051 //---------------------------------------- CHECK ------------------------------<<<
1052 //
1053 double theta = mlt->GetTheta(itr);
1054 double eta = -TMath::Log(TMath::Tan(theta/2));
1055 if (eta<fEtaMin || eta>fEtaMax) continue;
1056 //
1057 double dtheta = mlt->GetDeltaTheta(itr);
1058 double dThetaX = dtheta;
1059 if (fScaleDTBySin2T) {
1060 double sint = TMath::Sin(theta);
1061 dThetaX /= (sint*sint);
1062 }
1063 if (fCutOnDThetaX && TMath::Abs(dThetaX)>fDThetaWindow) continue;
1064 // double phi = mlt->GetPhi(itr);
1065 double dphi = mlt->GetDeltaPhi(itr);
1066 double dist = mlt->CalcDist(itr);
1067 //
1068 FillHistosSet(histos,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist);
1069 // special handling for mc info
1070 if (fillMC && fStack) {
1071 int lab0 = mlt->GetLabel(itr,0);
1072 int lab1 = mlt->GetLabel(itr,1);
1073 int typeMC = 2; // comb.bg.
1074 if (lab0 == lab1) typeMC = fStack->IsPhysicalPrimary(lab0) ? 0:1; // prim or sec
1075 if (typeMC==0) FillHistosSet(fHistosTrPrim,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // primary
1076 else if (typeMC==1) FillHistosSet(fHistosTrSec, eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // secondary
1077 else {
1078 FillHistosSet(fHistosTrComb,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // comb
1079 // for combinatorals fill also the uncorrelated part
1080 if (fMultReco) {
1081 float *trl = fMultReco->GetTracklet(itr);
1082 int clId0 = (int)trl[AliITSMultReconstructor::kClID1];
1083 int clId1 = (int)trl[AliITSMultReconstructor::kClID2];
1084 float *clLabs0 = fMultReco->GetClusterOfLayer(0,clId0) + AliITSMultReconstructor::kClMC0;
1085 float *clLabs1 = fMultReco->GetClusterOfLayer(1,clId1) + AliITSMultReconstructor::kClMC0;
1086 if (!HaveCommonParent(clLabs0,clLabs1))
1087 FillHistosSet(fHistosTrCombU,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist);
1088 }
1089 } // combinatorials
1090
1091 if (dist<fNStdCut) {
1092 double dphiS = TMath::Abs(dphi) - fDPhiShift; if (dphi<0) dphiS = -dphiS;
1093 if (dphiS<fDPhiSCut) FillSpecies(typeMC, lab0);
1094 }
1095 if (fCheckReconstructables) CheckReconstructables();
1096 }
1097 }
1098 //
1099 //-------------------------------------------------------------TMP RS - singles ------->>>
1100 int offsH = fCurrCentBin*kNStandardH;
1101 TH2* hSingles = (TH2*)histos->UncheckedAt(offsH+kHEtaZvSPD1);
1102 if (hSingles) {
1103 int nclS = mlt->GetNumberOfSingleClusters();
1104 double *thtS = mlt->GetThetaSingle();
1105 for (int ics=nclS;ics--;) {
1106 double etaS = -TMath::Log(TMath::Tan(thtS[ics]/2));
1107 if (etaS<fEtaMin || etaS>fEtaMax) continue;
1108 hSingles->Fill(etaS,fESDVtx[2]);
1109 }
1110 }
1111 //-------------------------------------------------------------TMP RS - singles -------<<<
1112 //
1113}
1114
1115//_________________________________________________________________________
c79a7c1c 1116void AliTrackletTaskMultipp::FillMCPrimaries()
c15ec07e 1117{
1118 // fill all MC primaries Zv vs Eta
c79a7c1c 1119 if (!fStack || !fMCevent) return;
c15ec07e 1120
1121 //---------------------------------------- CHECK ------------------------------>>>
1122 TArrayF vtxMC;
1123 AliGenHijingEventHeader* pyHeader = 0;
1124 //
1125 if (fUseMC) {
c79a7c1c 1126 pyHeader = (AliGenHijingEventHeader*) fMCevent->GenEventHeader();//header->GenEventHeader();
c15ec07e 1127 pyHeader->PrimaryVertex(vtxMC);
1128 }
1129 //---------------------------------------- CHECK ------------------------------<<<
1130 //
1131 int ntr = fStack->GetNtrack();
1132 TH2* hprimEtaZ = (TH2F*)fHistosCustom->UncheckedAt(kHZVEtaPrimMC+fCurrCentBin);
1133 int nprim = 0;
1134 for (int itr=ntr;itr--;) {
1135 if (!fStack->IsPhysicalPrimary(itr)) continue;
c79a7c1c 1136 AliMCParticle *part = (AliMCParticle*)fMCevent->GetTrack(itr);
c15ec07e 1137 if (!part->Charge()) continue;
1138 //
1139 //---------------------------------------- CHECK ------------------------------>>>
1140 /*
1141 Float_t dz = part->Zv() - vtxMC[2];
1142 if (TMath::Abs(dz)>1e-6) continue; // reject
1143 */
1144 //---------------------------------------- CHECK ------------------------------<<<
1145 //
1146 Float_t theta = part->Theta();
1147 if (theta<1e-6 || theta>TMath::Pi()-1e-6) continue;
1148 Float_t eta = part->Eta();
1149 if (eta<fEtaMin || eta>fEtaMax) continue;
1150 hprimEtaZ->Fill(eta, fESDVtx[2]);
1151 nprim++;
1152 }
1153 //
1154 ((TH1*)fHistosCustom->UncheckedAt(kHNPrimMeanMC))->Fill(fCurrCentBin,nprim);
1155 if (fNPart>0) {
1156 ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2PartMC))->Fill(fCurrCentBin,nprim/fNPart);
1157 ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2PartNpMC))->Fill(fNPart,nprim/fNPart);
1158 ((TH2*)fHistosCustom->UncheckedAt(kHNPartMC))->Fill(fNPart,fCurrCentBin);
1159 ((TH1*)fHistosCustom->UncheckedAt(kHNPartMeanMC))->Fill(fCurrCentBin,fNPart);
1160 }
1161 if (fNBColl>0) {
1162 ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2BCollMC))->Fill(fCurrCentBin,nprim/fNBColl);
1163 ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2BCollNpMC))->Fill(fNPart,nprim/fNBColl);
1164 ((TH2*)fHistosCustom->UncheckedAt(kHNBCollMC))->Fill(fNBColl,fCurrCentBin);
1165 ((TH1*)fHistosCustom->UncheckedAt(kHNBCollMeanMC))->Fill(fCurrCentBin,fNBColl);
1166 }
1167 //
1168}
1169
1170//_________________________________________________________________________
c79a7c1c 1171 void AliTrackletTaskMultipp::FillHistosSet(TObjArray* histos, double eta,
c15ec07e 1172 //double /*phi*/,double /*theta*/,
1173 double dphi,double dtheta,double dThetaX,
1174 double dist)
1175{
1176 // fill standard set of histos
1177 if (dist>fNStdDev) return;
1178 //
1179 double dphiS = TMath::Abs(dphi) - fDPhiShift;
1180 if (dphi<0) dphiS = -dphiS;
1181 //
1182 int offs = fCurrCentBin*kNStandardH;
1183 //
1184 if (histos->UncheckedAt(offs+kHDPhiSDThetaX))
1185 ((TH2*)histos->UncheckedAt(offs+kHDPhiSDThetaX))->Fill(dphiS,dThetaX);
1186 //
1187 if (histos->UncheckedAt(offs+kHDPhiDTheta))
1188 ((TH2*)histos->UncheckedAt(offs+kHDPhiDTheta))->Fill(dphi,dtheta);
1189 //
1190 if (histos->UncheckedAt(kHWDist))
1191 ((TH2*)histos->UncheckedAt(offs+kHWDist))->Fill(dist);
1192 //
1193 if (dist<fNStdCut && dphiS<fDPhiSCut && histos->UncheckedAt(offs+kHEtaZvCut))
1194 ((TH2*)histos->UncheckedAt(offs+kHEtaZvCut))->Fill(eta,fESDVtx[2]);
1195 //
1196}
1197
1198//__________________________________________________________________
c79a7c1c 1199void AliTrackletTaskMultipp::FillSpecies(Int_t primsec, Int_t id)
c15ec07e 1200{
1201 // fill PDGcode
1202 TH1 *hPart=0,*hParent=0;
1203 if (primsec==0) {
1204 hPart = (TH1*)fHistosCustom->UncheckedAt(kHPrimPDG);
1205 hParent = (TH1*)fHistosCustom->UncheckedAt(kHPrimParPDG);
1206 }
1207 else if (primsec==1) {
1208 hPart = (TH1*)fHistosCustom->UncheckedAt(kHSecPDG);
1209 hParent = (TH1*)fHistosCustom->UncheckedAt(kHSecParPDG);
1210 }
1211 else return;
1212 int ntr = fStack->GetNtrack();
1213 TParticle* part = fStack->Particle(id);
1214 int pdgCode = TMath::Abs(part->GetPdgCode());
1215 int pdgBin = GetPdgBin(pdgCode);
1216 int parID = part->GetFirstMother();
1217 int pdgCodePar = -1;
1218 int pdgBinPar = -1;
1219 while (parID>=0 && parID<ntr) {
1220 part = fStack->Particle(parID);
1221 pdgCodePar = TMath::Abs(part->GetPdgCode());
1222 parID = part->GetFirstMother();
1223 }
1224 if (pdgCodePar>0) pdgBinPar = GetPdgBin(pdgCodePar);
1225 //
1226 hPart->Fill(pdgBin,fCurrCentBin);
1227 hParent->Fill(pdgBinPar,fCurrCentBin);
1228 //
1229}
1230
1231//_________________________________________________________________________
c79a7c1c 1232Int_t AliTrackletTaskMultipp::GetCentralityBin(Float_t perc) const
c15ec07e 1233{
1234 // calculate centrality bin
1235 for (int i=0;i<fNCentBins;i++) if (perc>=fgkCentPerc[i] && perc<=fgkCentPerc[i+1]) return i;
1236 return -1;
1237}
1238
1239//_________________________________________________________________________
c79a7c1c 1240Int_t AliTrackletTaskMultipp::GetPdgBin(Int_t pdgCode)
c15ec07e 1241{
1242 // return my pdg bin
1243 int ncodes = sizeof(fgkPDGCodes)/sizeof(int);
1244 int pdgBin=0;
1245 for (pdgBin=0;pdgBin<ncodes;pdgBin++) if (pdgCode==fgkPDGCodes[pdgBin]) break;
1246 if (pdgBin>=ncodes) {
1247 if (float(pdgCode)>1e9) pdgBin = ncodes; // nuclei
1248 else pdgBin = ncodes+1; // unknown
1249 }
1250 return pdgBin;
1251}
1252
1253//_________________________________________________________________________
c79a7c1c 1254Bool_t AliTrackletTaskMultipp::HaveCommonParent(const float* clLabs0,const float* clLabs1)
c15ec07e 1255{
1256 // do 2 clusters have common parrent
1257 const int kMaxPar = 50;
1258 static int pars[2][50];
1259 int npars[2]={0,0};
1260 const float *labs[2] = {clLabs0,clLabs1};
1261 int ntr = fStack->GetNtrack();
3c78f321 1262 // printf("\nHaveCommonParent \n");
c15ec07e 1263 for (int il=0;il<2;il++) {
3c78f321 1264
c15ec07e 1265 for (int ilb=0;ilb<3;ilb++) {
1266 int lbl = (int)labs[il][ilb];
3c78f321 1267 if (lbl<2 || lbl>=ntr) continue;
c15ec07e 1268 //
1269 while (npars[il]<kMaxPar-1) {
c15ec07e 1270 TParticle* part = fStack->Particle(lbl);
1271 if (!part) break;
3c78f321 1272 int code = TMath::Abs(part->GetPdgCode());
1273 int q = (int)TMath::Abs(part->GetPDG()->Charge());
1274 if (code==21 || code<10 || q==1 || q==2 || q==4 ) break;
1275 //printf("%d/%d/%d: %4d (%d)%s|",il,ilb,npars[il],lbl,part->GetStatusCode(),part->GetName());
1276 pars[il][ npars[il]++ ] = lbl;
c15ec07e 1277 lbl = part->GetFirstMother();
1278 if (lbl<1 || lbl>=ntr) break;
1279 }
3c78f321 1280 // printf("\n");
c15ec07e 1281 }
3c78f321 1282 }
c15ec07e 1283 // compare array of parents
1284 for (int i0=npars[0];i0--;) for (int i1=npars[1];i1--;) if (pars[0][i0]==pars[1][i1]) return kTRUE;
1285 return kFALSE;
1286}
1287
1288
1289//_________________________________________________________________________
c79a7c1c 1290void AliTrackletTaskMultipp::CheckReconstructables()
c15ec07e 1291{
1292 // fill reconstructable tracklets hitsos
1293 static TArrayI trInd;
1294 static TBits isPrimArr;
1295 //
1296 if (!fMultReco || !fMultReco->IsRecoDone()) {AliInfo("To check reconstructables the reco had to be requested"); return;}
1297 if (!fStack) {AliInfo("MC Stack is not availalble"); return;}
1298 const double kPtMin = 0.05;
1299 //
1300 TClonesArray *clArr[2];
1301 for (int ilr=0;ilr<2;ilr++) {
1302 clArr[ilr] = fMultReco->GetClustersOfLayer(ilr);
1303 if (!clArr[ilr]) {AliInfo("Clusters are not available"); return;}
1304 }
1305 //
1306 int ntr = fStack->GetNtrack();
1307 if (!ntr) return;
1308 trInd.Reset();
1309 if (trInd.GetSize()<ntr) trInd.Set(ntr);
1310 isPrimArr.ResetAllBits();
1311 // count track wich may be reconstructable
1312 //
1313 int ntrStore=0,ntrStorePrim=0;
1314 Int_t *trIndArr = trInd.GetArray();
1315 for (Int_t it=0; it<ntr; it++) {
1316 TParticle* part = fStack->Particle(it);
1317 if (TMath::Abs(part->Eta())>2.2) continue;
1318 if (TMath::Abs(part->Pt())<kPtMin) continue;
1319 if (fStack->IsPhysicalPrimary(it)) {
1320 isPrimArr.SetBitNumber(it);
1321 ntrStorePrim++;
1322 }
1323 else { // check if secondary is worth cheking
1324 TParticlePDG* pdgPart = part->GetPDG();
1325 if (TMath::Abs(pdgPart->Charge())!=3) continue;
1326 if (part->R()>5.) continue;
1327 }
1328 trIndArr[it] = ++ntrStore;
1329 }
1330 //
1331 AliInfo(Form("Selected %d MC particles (%d prim) out of %d in the stack\n",ntrStore,ntrStorePrim,ntr));
1332 //
1333 const int kMaxCl=3;
1334 AliITSRecPoint **clIndL[2];
1335 clIndL[0] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer
1336 clIndL[1] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer
1337 memset(clIndL[0],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*));
1338 memset(clIndL[1],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*));
1339 //
1340 for (int ilr=0;ilr<2;ilr++) {
1341 TClonesArray *clusters = clArr[ilr];
1342 int ncl = clusters->GetEntriesFast();
1343 for (int icl=ncl;icl--;) {
1344 AliITSRecPoint *cl = (AliITSRecPoint*)clusters->UncheckedAt(icl);
1345 for (int ilb=3;ilb--;) {
1346 int lbl = cl->GetLabel(ilb); if (lbl<0 || lbl>=ntr) continue;
1347 int lblI = trIndArr[lbl];
1348 if (--lblI<0) continue; // not kept
1349 for (int icc=0;icc<kMaxCl;icc++) if (!clIndL[ilr][lblI+icc*ntrStore]) {clIndL[ilr][lblI+ntrStore*icc] = cl; break;} // first empty one
1350 }
1351 }
1352 }
1353 //
1354 Float_t clusterLay[2][AliITSMultReconstructor::kClNPar];
1355 double trComp[6][kMaxCl*kMaxCl];
1356 int indQual[kMaxCl*kMaxCl];
1357 //
1358 for (int itr=ntr;itr--;) {
1359 int lblI = trIndArr[itr];
1360 if (--lblI<0) continue; // discarded
1361 int ntrCand = 0; // number of tracklet candidates (including overlaps)
1362 for (int icl0=0;icl0<kMaxCl;icl0++) {
1363 AliITSRecPoint *cl0 = clIndL[0][lblI+icl0*ntrStore];
1364 if (!cl0 || !clIndL[1][lblI]) break;
1365 cl0->GetGlobalXYZ( clusterLay[0] );
1366 fMultReco->ClusterPos2Angles(clusterLay[0], fESDVtx);
1367 for (int icl1=0;icl1<kMaxCl;icl1++) {
1368 AliITSRecPoint *cl1 = clIndL[1][lblI+icl1*ntrStore];
1369 if (!cl1) break;
1370 cl1->GetGlobalXYZ( clusterLay[1] );
1371 fMultReco->ClusterPos2Angles(clusterLay[1], fESDVtx);
1372 trComp[AliITSMultReconstructor::kTrPhi][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClPh];
1373 trComp[AliITSMultReconstructor::kTrTheta][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClTh];
1374 trComp[AliITSMultReconstructor::kTrDTheta][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClTh] - clusterLay[1][AliITSMultReconstructor::kClTh];
1375 trComp[AliITSMultReconstructor::kTrDPhi][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClPh] - clusterLay[1][AliITSMultReconstructor::kClPh];
1376 trComp[AliITSMultReconstructor::kTrLab1][ntrCand] = icl1*10 + icl0;
1377 double &dphi = trComp[ntrCand][3];
1378 if (dphi>TMath::Pi()) dphi=2.*TMath::Pi()-dphi; // take into account boundary condition
1379 trComp[5][ntrCand] = fMultReco->CalcDist(trComp[AliITSMultReconstructor::kTrDPhi][ntrCand],
1380 trComp[AliITSMultReconstructor::kTrDTheta][ntrCand],
1381 trComp[AliITSMultReconstructor::kTrTheta][ntrCand]);
1382 ntrCand++;
1383 }
1384 }
1385 if (!ntrCand) continue; // no tracklets
1386 if (ntrCand>1) TMath::Sort(ntrCand,trComp[5],indQual,kFALSE); else indQual[0] = 0; // sort in weighted distance
1387 if (fRemoveOverlaps) ntrCand = 1; // select the best
1388 //
1389 // disable worst tracklet with shared cluster
1390 for (int itc=0;itc<ntrCand;itc++) {
1391 int ind = indQual[itc];
1392 if (trComp[AliITSMultReconstructor::kTrLab1][ind]<0) continue; // already disabled
1393 for (int jtc=itc+1;jtc<ntrCand;jtc++) {
1394 int jnd = indQual[jtc];
1395 if (trComp[AliITSMultReconstructor::kTrLab1][jnd]<0) continue; // already disabled
1396 if ( int(trComp[AliITSMultReconstructor::kTrLab1][ind])/10 == int(trComp[AliITSMultReconstructor::kTrLab1][jnd])/10 ||
1397 int(trComp[AliITSMultReconstructor::kTrLab1][ind])%10 == int(trComp[AliITSMultReconstructor::kTrLab1][jnd])%10) trComp[AliITSMultReconstructor::kTrLab1][jnd] = -1;
1398 }
1399 }
1400 //
1401 // store, but forbid cluster reusing
1402 TObjArray* histos = isPrimArr.TestBitNumber(itr) ? fHistosTrRcblPrim : fHistosTrRcblSec;
1403 for (int itc=0;itc<ntrCand;itc++) {
1404 int ind = indQual[itc];
1405 if (trComp[4][ind]<0) continue; // discarded
1406 double eta = -TMath::Log(TMath::Tan(trComp[AliITSMultReconstructor::kTrTheta][ind]/2));
1407 if (eta<fEtaMin || eta>fEtaMax) continue;
1408 double dThetaX = trComp[AliITSMultReconstructor::kTrTheta][ind];
1409 if (fScaleDTBySin2T) {
1410 double sint = TMath::Sin(trComp[AliITSMultReconstructor::kTrTheta][ind]);
1411 dThetaX /= (sint*sint);
1412 }
1413 FillHistosSet(histos,eta,
1414 //trComp[AliITSMultReconstructor::kTrPhi][ind],trComp[AliITSMultReconstructor::kTrTheta][ind],
1415 trComp[AliITSMultReconstructor::kTrDPhi][ind],trComp[AliITSMultReconstructor::kTrDTheta][ind],
1416 dThetaX,trComp[5][ind]);
1417 }
1418 }
1419 //
1420 delete[] clIndL[0];
1421 delete[] clIndL[1];
1422}
1423
1424//_________________________________________________________________________
c79a7c1c 1425void AliTrackletTaskMultipp::FillClusterInfo()
c15ec07e 1426{
1427 // fill info on clusters associated to good tracklets
1428 if (!fMultReco) return;
1429 int ntr = fMultReco->GetNTracklets();
1430 int clID[2];
1431 TH2F *hclU[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL1)};
1432 TH2F *hclA[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL1)};
1433 for (int itr=ntr;itr--;) {
1434 Float_t *trc = fMultReco->GetTracklet(itr);
1435 if (TMath::Abs(TMath::Abs(trc[AliITSMultReconstructor::kTrDPhi])-fDPhiShift)>fDPhiSCut) continue;
1436 if (fMultReco->CalcDist(trc[AliITSMultReconstructor::kTrDPhi],
1437 trc[AliITSMultReconstructor::kTrDTheta],
1438 trc[AliITSMultReconstructor::kTrTheta]) > fNStdCut) continue;
1439 clID[0] = (int)trc[AliITSMultReconstructor::kClID1];
1440 clID[1] = (int)trc[AliITSMultReconstructor::kClID2];
1441 for (int il=0;il<2;il++) {
1442 Float_t *clinf = fMultReco->GetClusterOfLayer(il,clID[il]);
1443 hclU[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]);
1444 }
1445 }
1446 //
1447 for (int il=0;il<2;il++) for (int ic=fMultReco->GetNClustersLayer(il);ic--;) {
1448 Float_t *clinf = fMultReco->GetClusterOfLayer(il,ic);
1449 hclA[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]);
1450 }
1451 //
1452}
1453
1454//_________________________________________________________________________
c79a7c1c 1455void AliTrackletTaskMultipp::FillClusterInfoFromMult(const AliMultiplicity* mlt, double zVertex)
c15ec07e 1456{
1457 // fill info on clusters taking them from Multiplicity object
1458 const double kRSPD1 = 3.9;
1459 TH2F *hclU = (TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0);
1460 TH2F *hclA = (TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0);
1461 int ntr = mlt->GetNumberOfTracklets();
1462 for (int itr=ntr;itr--;) {
1463 Bool_t goodTracklet = kTRUE;
1464 if (TMath::Abs( mlt->GetDeltaPhi(itr)-fDPhiShift)>fDPhiSCut) goodTracklet = kFALSE;
1465 if (mlt->CalcDist(itr) > fNStdCut) goodTracklet = kFALSE;
1466 double phi = mlt->GetPhi(itr);
1467 double z = kRSPD1/TMath::Tan(mlt->GetTheta(itr)) + zVertex;
1468 if (goodTracklet) hclU->Fill(z,phi);
1469 hclA->Fill(z,phi);
1470 }
1471 //
1472 int ncl = mlt->GetNumberOfSingleClusters();
1473 for (int icl=ncl;icl--;) {
1474 double phi = mlt->GetPhiSingle(icl);
1475 double z = kRSPD1/TMath::Tan(mlt->GetThetaSingle(icl)) + zVertex;
1476 hclA->Fill(z,phi);
1477 }
1478 //
1479}
1480
1481//_________________________________________________________________________
c79a7c1c 1482void AliTrackletTaskMultipp::FillClusterAutoCorrelationFromMult(const AliMultiplicity* mlt, double zVertex)
c15ec07e 1483{
1484 // fill mutual distance between SPD1 clusters
1485 const double kRSPD1 = 3.9;
1486 enum {kX=0,kY,kZ,kPhi,kNV};
1487 //
1488 TH2F *hclDstZAll = (TH2F*)fHistosCustom->UncheckedAt(kHclDstZAll);
1489 TH2F *hclDstPhiAll = (TH2F*)fHistosCustom->UncheckedAt(kHclDstPhiAll);
1490 TH2F *hclDstZUsed = (TH2F*)fHistosCustom->UncheckedAt(kHclDstZUsed);
1491 TH2F *hclDstPhiUsed = (TH2F*)fHistosCustom->UncheckedAt(kHclDstPhiUsed);
1492 if (!hclDstZAll && !hclDstPhiAll && !hclDstZUsed && !hclDstPhiUsed) return; // nothing to fill
1493
1494 int nCl = mlt->GetNumberOfTracklets() + mlt->GetNumberOfSingleClusters();
1495 if (nCl<2) return;
1496 double *clXYZPhi = new Double_t[kNV*nCl];
1497 int cnt = 0;
1498 for (int itr=mlt->GetNumberOfTracklets();itr--;) { // extract coordinates of used SPD1 clusters
1499 double phi = mlt->GetPhi(itr); if (phi<0) phi += 2*TMath::Pi();
1500 int offs = cnt*kNV;
1501 clXYZPhi[offs + kX] = TMath::Cos(phi)*kRSPD1;
1502 clXYZPhi[offs + kY] = TMath::Sin(phi)*kRSPD1;
1503 clXYZPhi[offs + kZ] = kRSPD1/TMath::Tan(mlt->GetTheta(itr)) + zVertex;
1504 clXYZPhi[offs + kPhi] = phi;
1505 cnt++;
1506 }
1507 //
1508 for (int icl=mlt->GetNumberOfSingleClusters();icl--;) { // extract coordinates of unused SPD1 clusters
1509 double phi = mlt->GetPhiSingle(icl); if (phi<0) phi += 2*TMath::Pi();
1510 int offs = cnt*kNV;
1511 clXYZPhi[offs + kX] = TMath::Cos(phi)*kRSPD1;
1512 clXYZPhi[offs + kY] = TMath::Sin(phi)*kRSPD1;
1513 clXYZPhi[offs + kZ] = kRSPD1/TMath::Tan(mlt->GetThetaSingle(icl)) + zVertex;
1514 clXYZPhi[offs + kPhi] = phi;
1515 cnt++;
1516 }
1517 //
1518 for (int icl=0;icl<nCl;icl++) {
1519 int offs = icl*kNV;
1520 if (hclDstZAll) hclDstZAll->Fill(clXYZPhi[offs+kZ],-1,2); // fill in underflow, to count the clusters
1521 if (hclDstPhiAll) hclDstPhiAll->Fill(clXYZPhi[offs+kPhi],-1,2); // fill in underflow, to count the clusters
1522 if (icl<mlt->GetNumberOfTracklets()) {
1523 if (hclDstZUsed) hclDstZUsed->Fill(clXYZPhi[offs+kZ],-1); // fill in underflow, to count the clusters
1524 if (hclDstPhiUsed) hclDstPhiUsed->Fill(clXYZPhi[offs+kPhi],-1); // fill in underflow, to count the clusters
1525 }
1526 for (int icl1=icl+1;icl1<nCl;icl1++) {
1527 int offs1 =icl1*kNV;
1528 double dx = clXYZPhi[offs+kX]-clXYZPhi[offs1+kX];
1529 double dy = clXYZPhi[offs+kY]-clXYZPhi[offs1+kY];
1530 double dz = clXYZPhi[offs+kZ]-clXYZPhi[offs1+kZ];
1531 double dst = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
1532 //
1533 if (hclDstZAll) {
1534 hclDstZAll->Fill(clXYZPhi[offs+kZ],dst);
1535 hclDstZAll->Fill(clXYZPhi[offs1+kZ],dst);
1536 }
1537 if (hclDstPhiAll) {
1538 hclDstPhiAll->Fill(clXYZPhi[offs+kPhi],dst);
1539 hclDstPhiAll->Fill(clXYZPhi[offs1+kPhi],dst);
1540 }
1541 //
1542 // tracklets with singles only
1543 if (icl>=mlt->GetNumberOfTracklets() || icl1<mlt->GetNumberOfTracklets()) continue;
1544 if (hclDstZUsed) hclDstZUsed->Fill(clXYZPhi[offs+kZ],dst);
1545 if (hclDstPhiUsed) hclDstPhiUsed->Fill(clXYZPhi[offs+kPhi],dst);
1546 //
1547 }
1548 }
1549
1550}