bug fix
[u/mrichter/AliRoot.git] / PWGPP / forward / SPDClustTask / AliSPDClustTask.cxx
CommitLineData
d854076c 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 AliSPDClustTask //
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 "TList.h"
35#include "TObjArray.h"
36#include "TGeoGlobalMagField.h"
37
38#include "AliAnalysisManager.h"
39
40#include "AliMultiplicity.h"
41#include "AliESDEvent.h"
42#include "AliESDInputHandler.h"
43#include "AliESDInputHandlerRP.h"
44#include "AliCDBPath.h"
45#include "AliCDBManager.h"
46#include "AliCDBEntry.h"
47#include "AliCDBStorage.h"
48#include "AliGeomManager.h"
49#include "AliMagF.h"
50#include "AliESDVZERO.h"
51#include "AliESDZDC.h"
52#include "AliRunLoader.h"
53#include "AliMCEventHandler.h"
54#include "AliMCEvent.h"
55#include "AliMCParticle.h"
56#include "AliStack.h"
57#include "AliGenEventHeader.h"
58#include "AliCentrality.h"
59#include "../ITS/AliITSRecPoint.h"
60#include "../ITS/AliITSgeomTGeo.h"
61#include "../ITS/AliITSMultReconstructor.h"
62
63#include "AliLog.h"
64
65#include "AliPhysicsSelection.h"
66#include "AliSPDClustTask.h"
67#include "AliITSMultReconstructor.h"
68#include "AliGenEventHeader.h"
69#include "AliGenHijingEventHeader.h"
70#include "AliGenDPMjetEventHeader.h"
71#include "AliESDtrackCuts.h"
72
73ClassImp(AliSPDClustTask)
74
75
76//________________________________________________________________________
77/*//Default constructor
78AliSPDClustTask::AliSPDClustTask(const char *name)
79 : AliAnalysisTaskSE(name),
80*/
81//________________________________________________________________________
82AliSPDClustTask::AliSPDClustTask(const char *name)
83 : AliAnalysisTaskSE(name),
84//
85 fOutput(0),
86 fUseMC(kFALSE),
87 fHistos(0),
88//
89 fEtaMin(-3.0),
90 fEtaMax(3.0),
91 fZVertexMin(-20),
92 fZVertexMax( 20),
93//
94 fScaleDTBySin2T(kFALSE),
95 fCutOnDThetaX(kFALSE),
96 fNStdDev(1.),
97 fDPhiWindow(0.08),
98 fDThetaWindow(0.025),
99 fDPhiShift(0.0045),
100 fPhiOverlapCut(0.005),
101 fZetaOverlap(0.05),
102 fRemoveOverlaps(kFALSE),
103//
104 fDPhiSCut(0.06),
105 fNStdCut(1.),
106//
107 fMultReco(0),
108 fRPTree(0),
109 fStack(0),
110 fMCEvent(0),
111 //
112 fDontMerge(kFALSE),
488209ed 113 fhPtPionIn(0),
114 fhPtKaonIn(0),
115 fhPtProtonIn(0)
d854076c 116{
117 // Constructor
118
119 DefineOutput(1, TList::Class());
120 //
121 SetScaleDThetaBySin2T();
122 SetNStdDev();
123 SetPhiWindow();
124 SetThetaWindow();
125 SetPhiShift();
126 SetPhiOverlapCut();
127 SetZetaOverlapCut();
128 SetRemoveOverlaps();
129 //
130}
131
132//________________________________________________________________________
133AliSPDClustTask::~AliSPDClustTask()
134{
135 // Destructor
136 // histograms are in the output list and deleted when the output
137 // list is deleted by the TSelector dtor
138 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { //RRR
139 printf("Deleteing output\n");
140 delete fOutput;
141 fOutput = 0;
142 }
143 //
144 delete fMultReco;
145 delete fHistos;
488209ed 146 if (fhPtPionIn) delete fhPtPionIn;
147 if (fhPtKaonIn) delete fhPtKaonIn;
148 if (fhPtProtonIn) delete fhPtProtonIn;
d854076c 149 //
150}
151
152//________________________________________________________________________
153void AliSPDClustTask::UserCreateOutputObjects()
154{
155 //
156 // OpenFile(1); fDontMerge = kTRUE;
157 //
158
159 fOutput = new TList();
160 fOutput->SetOwner();
161 //
162 AliCDBManager *man = AliCDBManager::Instance();
163 if (fUseMC) {
164 Bool_t newGeom = kTRUE;
165 man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
166 if (newGeom) {
167 // new geom
168 AliCDBEntry* obj = man->Get("GRP/Geometry/Data",130844,8);
169 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
170 if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,6,-1)) AliFatal("Failed to misalign geometry");
171 }
172 else {
173 // old geom
174 AliCDBEntry* obj = man->Get("GRP/Geometry/Data",130845,7);
175 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
176 if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130845,5,-1)) AliFatal("Failed to misalign geometry");
177 }
178 }
179 else {
180 man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB"); //man->SetRun(137366);
181 AliCDBEntry* obj = man->Get("GRP/Geometry/Data",137366);
182 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
183 if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137366,-1,-1)) AliFatal("Failed to misalign geometry");
184 }
185 //
186 // Create histograms
187 fHistos = BookHistos();
188 for(Int_t i=0; i<fHistos->GetEntriesFast(); i++) {
189 if(fHistos->At(i)){
190 fOutput->AddLast(fHistos->At(i));
191 }
192 }
193
194 PostData(1, fOutput);
195 //
196}
197
198//________________________________________________________________________
199void AliSPDClustTask::UserExec(Option_t *)
200{
201 // Main loop
202 //
203 //
204 AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
205 fRPTree = 0;
206 AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler();
207 if (!handRP) { AliError("No RP handler"); return; }
208 //
209 fRPTree = handRP->GetTreeR("ITS");
210 if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
211 //
212 AliESDEvent *esd = handRP->GetEvent();
213 if (!esd) { AliError("No AliESDEvent"); return; }
214 //
215 AliMCEventHandler* eventHandler = 0;
216 fMCEvent = 0;
217 fStack = 0;
218 //
219 if (fUseMC) {
220 eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
221 if (!eventHandler) { AliError("Could not retrieve MC event handler"); return; }
222 fMCEvent = eventHandler->MCEvent();
223 if (!fMCEvent) { AliError("Could not retrieve MC event"); return; }
224 fStack = fMCEvent->Stack();
225 if (!fStack) { AliError("Stack not available"); return; }
226 }
227 //
228 const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
229 if (vtxESD->GetNContributors()<1) return;
230 //
231 // ------------------RS: Put here additional event selection if needed ------->>
232 // at the moment, I am just selecting on vertex
233 TString vtxTyp = vtxESD->GetTitle();
234 if (vtxTyp.Contains("vertexer: Z")) {
235 if (vtxESD->GetDispersion()>0.04) return;
236 if (vtxESD->GetZRes()>0.25) return;
237 }
238 //
239 if (vtxESD->GetZ()<fZVertexMin || vtxESD->GetZ()>fZVertexMax) return;
240 //
241 // pile-up rejection
242 if (esd->IsPileupFromSPD(3,0.8)) return;
243 //
244 // ------------------RS: Put here additional event selection if needed -------<<
245 //
246 fVtx[0] = vtxESD->GetX();
247 fVtx[1] = vtxESD->GetY();
248 fVtx[2] = vtxESD->GetZ();
249 //
250 // do we need to initialize the field?
251 AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
252 if (!field && !esd->InitMagneticField()) {printf("Failed to initialize the B field\n");return;}
253 //
254 AliStack *stack = NULL;
255 if (fUseMC) {
256 if (!fMCEvent) { AliError("Could not retrieve MC event"); return; }
257 stack = fMCEvent->Stack();
258 if (!stack) { AliError("ERROR: Could not get stack"); return; }
259 for (Int_t iPart = 0; iPart < stack->GetNtrack(); ++iPart) {
260 if (!stack->IsPhysicalPrimary(iPart)) continue;
261 TParticle* p = stack->Particle(iPart);
262 if (TMath::Abs(p->Y()) > 0.5) continue;
263 if (TMath::Abs(p->GetPdgCode()) != 211) continue;
264 TH1D *hPt = (TH1D*)fHistos->At(kHPt);
265 hPt->Fill(p->Pt());
266 }
267 }
268 //
269 InitMultReco();
270 fMultReco->Reconstruct(fRPTree, fVtx);
271 // AliMultiplicity* mlt = fMultReco->GetMultiplicity();
272 FillHistos(stack);
273 //
274 delete fMultReco;
275 fMultReco = 0;
276 //
277}
278
279//________________________________________________________________________
280void AliSPDClustTask::Terminate(Option_t *)
281{
282 Printf("Terminating...");
283}
284
285//_________________________________________________________________________
286void AliSPDClustTask::InitMultReco()
287{
288 // create mult reconstructor
289 if (fMultReco) delete fMultReco;
290 fMultReco = new AliITSMultReconstructor();
291 fMultReco->SetCreateClustersCopy(kTRUE);
292 fMultReco->SetScaleDThetaBySin2T(fScaleDTBySin2T);
293 fMultReco->SetNStdDev(fNStdDev);
294 fMultReco->SetPhiWindow( fDPhiWindow );
295 fMultReco->SetThetaWindow( fDThetaWindow );
296 fMultReco->SetPhiShift( fDPhiShift );
297 fMultReco->SetRemoveClustersFromOverlaps(fRemoveOverlaps);
298 fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
299 fMultReco->SetZetaOverlapCut(fZetaOverlap);
300 fMultReco->SetHistOn(kFALSE);
301 //
302}
303
304//_________________________________________________________________________
305TObjArray* AliSPDClustTask::BookHistos()
306{
307 // book set of histos
308 //
309 TObjArray* histos = new TObjArray();
310 histos->SetOwner(kFALSE);
311 // this array contains histos stored at slots:
312 // 0-99 : tracklet related histos
313 // 100-199 : SPD1 clusters not used by tracklets
314 // 200-299 : SPD2 clusters not used by tracklets
315 // 300-399 : SPD1 clusters used by tracklets
316 // 400-499 : SPD2 clusters used by tracklets
317 //
318 int nEtaBins = int((fEtaMax-fEtaMin)/0.1);
319 if (nEtaBins<1) nEtaBins = 1;
320 //
321 // just an example: cluster type vs eta
322 for (int iused=0;iused<2;iused++) {
323 for (int spd=0;spd<2;spd++) {
324 TH2F* h = new TH2F(Form("clType_SPD%d_%s",spd,iused ? "used":"free"),
325 Form("clType SPD%d %s",spd,iused ? "used":"free"),
326 nEtaBins, fEtaMin,fEtaMax,
327 15, -0.5, 14.5);
328 histos->AddAtAndExpand(h, kHClusters+iused*200+spd*100 + kClTypevsEta);
329 TH1F* hZ = new TH1F(Form("clZ_SPD%d_%s",spd,iused ? "used":"free"),
330 Form("clZ SPD%d %s",spd,iused ? "used":"free"),
331 200, -15, 15);
332 histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZ);
333 TH1F* hEta = new TH1F(Form("clEta_SPD%d_%s",spd,iused ? "used":"free"),
334 Form("clEta SPD%d %s",spd,iused ? "used":"free"),
335 nEtaBins, fEtaMin,fEtaMax);
336 histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEta);
337
338 hZ = new TH1F(Form("clZ_pions_weighted_SPD%d_%s",spd,iused ? "used":"free"),
339 Form("clZ pions weighted SPD%d %s",spd,iused ? "used":"free"),
340 200, -15, 15);
341 histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZPionsW);
342 hEta = new TH1F(Form("clEta_pions_weighted_SPD%d_%s",spd,iused ? "used":"free"),
343 Form("clEta pions weighted SPD%d %s",spd,iused ? "used":"free"),
344 nEtaBins, fEtaMin,fEtaMax);
345 histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEtaPionsW);
346
347 hZ = new TH1F(Form("clZ_pions_SPD%d_%s",spd,iused ? "used":"free"),
348 Form("clZ pions SPD%d %s",spd,iused ? "used":"free"),
349 200, -15, 15);
350 histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZPions);
351 hEta = new TH1F(Form("clEta_pions_SPD%d_%s",spd,iused ? "used":"free"),
352 Form("clEta pions SPD%d %s",spd,iused ? "used":"free"),
353 nEtaBins, fEtaMin,fEtaMax);
354 histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEtaPions);
355 }
356 }
357 if (fUseMC) {
488209ed 358 TH1D *hPt = new TH1D("hPt","Pions Pt spectra (MC)",fhPtPionIn->GetXaxis()->GetNbins(),fhPtPionIn->GetXaxis()->GetXbins()->GetArray());
d854076c 359 histos->AddAtAndExpand(hPt, kHPt);
360 }
361 //
362 return histos;
363}
364
365//_________________________________________________________________________
366void AliSPDClustTask::FillHistos(AliStack *stack)
367{
368 // fill info on clusters, separately on associated to tracklets and singles
369 const int kUsedBit = BIT(15);
370 //
371 // 0) --------- just in case: reset kUsedBit of clusters --->>>
372 for (int ilr=0;ilr<2;ilr++) {
373 for (int icl=fMultReco->GetNClustersLayer(ilr);icl--;) { // loop on clusters of layer
374 AliITSRecPoint *clus = fMultReco->GetRecPoint(ilr,icl);
375 if (!clus) {
376 AliError(Form("Failed to extract cluster %d of layer %d",icl,ilr));
377 return;
378 }
379 clus->ResetBit(kUsedBit);
380 } // loop on clusters of layer
381 } // loop on layers
382 //
383 // 1) --------- mark used clusters by kUsedBit ------------>>>
384 int ntr = fMultReco->GetNTracklets();
385 for (int itr=ntr;itr--;) {
386 Float_t *trc = fMultReco->GetTracklet(itr);
387 AliITSRecPoint *cl0 = fMultReco->GetRecPoint(0,(int)trc[AliITSMultReconstructor::kClID1]); // cluster on SPD1
388 AliITSRecPoint *cl1 = fMultReco->GetRecPoint(1,(int)trc[AliITSMultReconstructor::kClID2]); // cluster on SPD1
389 //
390 if (!cl0 || !cl1) {
391 AliError(Form("Failed to extract clusters associated with tracklet %d",itr));
392 continue;
393 }
394 cl0->SetBit(kUsedBit);
395 cl1->SetBit(kUsedBit);
396 }
397 //
398 // 2) --------- fill needed info for used and unuses clusters
399 TH2F* h2d = 0;
400 TH1F* h1d = 0;
401 //
402 for (int ilr=0;ilr<2;ilr++) { // loop on layers
403 int ncl = fMultReco->GetNClustersLayer(ilr);
404 for (int icl=ncl;icl--;) { // loop on clusters of layer
405 AliITSRecPoint *clus = fMultReco->GetRecPoint(ilr,icl); // extract the cluster to access clType
406 Float_t* clInfo = fMultReco->GetClusterOfLayer(ilr,icl); // we already have in the MultReco cluster params extracted and
407 // stored in the float array, see AliITSMultReconstructor
408 //
409 float phi = clInfo[AliITSMultReconstructor::kClPh];
410 float eta = -TMath::Log(TMath::Tan(clInfo[AliITSMultReconstructor::kClTh]/2));
411 float z = clInfo[AliITSMultReconstructor::kClZ];
412 //
413 int used = clus->TestBit(kUsedBit) ? 1:0;
414 //
415 // Annalisa, here you should fill the info on used/unused clusters
416 if (phi > 0 && phi < 1.6) {
417 h2d = (TH2F*)fHistos->At(kHClusters+used*200+ilr*100 + kClTypevsEta);
418 h2d->Fill(eta, clus->GetClusterType());
419 h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZ);
420 h1d->Fill(z);
421 h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEta);
422 h1d->Fill(eta);
488209ed 423 Double_t weight = 1.;
d854076c 424 if (stack) {
425 for(Int_t iLabel = 0; iLabel < 3; ++iLabel) {
426 Int_t clLabel = clus->GetLabel(iLabel);
427 if (clLabel <= 0) continue;
428 Int_t moLabel = FindMotherParticle(stack, clLabel);
429 if (moLabel <= 0) continue;
430 TParticle* p = stack->Particle(moLabel);
488209ed 431 if ((TMath::Abs(p->GetPdgCode()) != 211) &&
432 (TMath::Abs(p->GetPdgCode()) != 321) &&
433 (TMath::Abs(p->GetPdgCode()) != 2212)) continue;
434 weight = PtWeight(p->Pt(),p->GetPdgCode());
d854076c 435 h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZPionsW);
436 h1d->Fill(z,weight);
437 h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEtaPionsW);
438 h1d->Fill(eta,weight);
439 h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZPions);
440 h1d->Fill(z);
441 h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEtaPions);
442 h1d->Fill(eta);
443 break;
444 }
445 }
446 }
447 } // loop on clusters of layer
448 } // loop on layers
449 //
450}
451
452//________________________________________________________________________
453void AliSPDClustTask::SetInput(const char *filename)
454{
455 // Read pions Pt spectra
456 // from a local input file
457
458 AliInfo(Form("Reading input histograms from %s",filename));
459 if (fUseMC) {
488209ed 460 TFile *fPt = TFile::Open("FitOutput.root");
d854076c 461 if (!fPt) { AliError("Failed to open input file"); return; }
488209ed 462 fhPtPionIn = (TH1F*)fPt->Get("pionPtWeight")->Clone("fhPtPionIn");
463 fhPtPionIn->SetDirectory(0);
464 fhPtKaonIn = (TH1F*)fPt->Get("kaonPtWeight")->Clone("fhPtKaonIn");
465 fhPtKaonIn->SetDirectory(0);
466 fhPtProtonIn = (TH1F*)fPt->Get("protPtWeight")->Clone("fhPtProtonIn");
467 fhPtProtonIn->SetDirectory(0);
d854076c 468 fPt->Close();
469 }
470}
471
472//________________________________________________________________________
473Int_t AliSPDClustTask::FindMotherParticle(AliStack* stack, Int_t i)
474{
475 if(stack->IsPhysicalPrimary(i)) return i;
476 TParticle* p = stack->Particle(i);
477 Int_t imo = p->GetFirstMother();
478
479 // printf("imo = %d\n",imo);
480
481 if(imo<=0)
482 {
483 AliInfo("particle is not a PhysPrim and has no mother");
484 return imo;
485 }
486 return FindMotherParticle(stack, imo);
487
488}
489
490//________________________________________________________________________
488209ed 491Double_t AliSPDClustTask::PtWeight(Double_t pt, Int_t pdgCode)
d854076c 492{
ec9cb017 493 switch (TMath::Abs(pdgCode)) {
488209ed 494 case 211:
495 return fhPtPionIn->GetBinContent(fhPtPionIn->FindBin(pt));
496 case 321:
497 return fhPtKaonIn->GetBinContent(fhPtKaonIn->FindBin(pt));
498 case 2212:
499 return fhPtProtonIn->GetBinContent(fhPtProtonIn->FindBin(pt));
500 default:
d854076c 501 return 1.;
488209ed 502 }
d854076c 503}