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