]>
Commit | Line | Data |
---|---|---|
6ad0bfa0 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, 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 | ||
b2a60966 | 16 | /* $Id$ */ |
17 | ||
6ad0bfa0 | 18 | //_________________________________________________________________________ |
efad3788 | 19 | // Algorythm class to analyze PHOS events. In this class we demostrate, |
20 | // how to handle reconstructed objects with AliPHSOIndexToObject. | |
21 | // As an example we propose sulotions for four most frequently used tasks: | |
22 | // DrawRecon(...) - to draw reconstructed objects in the PHOS plane, | |
23 | // very usefull in the debuging | |
24 | // InvarianMass(...) - to calculate "REAL" and "MIXED" photon pairs | |
25 | // invariant mass distributions | |
26 | // EnergyResoluition(...) -\ Energy and position resolutions of the | |
27 | // PositionResolution(...)-/ reconstructed photons | |
28 | // Contamination(...) - calculates contamination of the photon spectrum and | |
29 | // pobability of reconstruction of several primaries as | |
30 | // kGAMMA,kELECTRON etc. | |
31 | // | |
32 | // User Case: | |
e4761a49 | 33 | // root [0] AliPHOSAnalyze * a = new AliPHOSAnalyze("galice.root") |
efad3788 | 34 | // // set the file you want to analyse |
35 | // root [1] a->DrawRecon(1,3) | |
36 | // // plot RecObjects, made in event 1, PHOS module 3 | |
37 | // root [2] a->DrawRecon(1,3,"PHOSRP","another PID") | |
38 | // // plot RecObjets made in the event 1, PHOS module 3, | |
39 | // // produced in the another reconstruction pass, | |
40 | // // which produced PHOS RecParticles ("PHOSRP") with | |
41 | // // title "another PID". | |
42 | // root [3] a->InvariantMass() | |
43 | // // Calculates "REAL" and "MIXED" invariant mass | |
44 | // // distributions of kGAMMA and (kGAMMA+kNEUTRALEM) | |
45 | // // and APPENDS this to the file "invmass.root" | |
46 | // root [4] a->PositionResolution() | |
47 | // // calculates two dimentional histos: energy of the primary | |
48 | // // photon vs distance betwin incedence point and reconstructed | |
49 | // // poisition. One can analyse the produced file position.root | |
50 | // // with macro PhotonPosition.C | |
51 | // root [5] a->EnergyResolution() | |
52 | // // calculates two dimentional histos: energy of the primary | |
53 | // // photon vs energy of the reconstructed particle. One can | |
54 | // // analyse the produced file energy.root | |
55 | // // with macro PhotonEnergy.C | |
56 | // root [6] a->Contamination() | |
57 | // // fills spectra of primary photons and several kinds of | |
58 | // // reconstructed particles, so that analyzing them one can | |
59 | // // estimate conatmination, efficiency of registration etc. | |
a3dfe79c | 60 | //*-- |
efad3788 | 61 | //*-- Author: Dmitri Peressounko (SUBATECH & RRC Kurchatov Institute) |
6ad0bfa0 | 62 | ////////////////////////////////////////////////////////////////////////////// |
63 | ||
baef0810 | 64 | |
6ad0bfa0 | 65 | // --- ROOT system --- |
66 | ||
92862013 | 67 | #include "TFile.h" |
68 | #include "TH1.h" | |
6ad0bfa0 | 69 | #include "TH2.h" |
83448140 | 70 | #include "TH2.h" |
6ad0bfa0 | 71 | #include "TParticle.h" |
72 | #include "TClonesArray.h" | |
73 | #include "TTree.h" | |
74 | #include "TMath.h" | |
a88eea6f | 75 | #include "TROOT.h" |
76 | #include "TFolder.h" | |
6ad0bfa0 | 77 | |
78 | // --- Standard library --- | |
79 | ||
efad3788 | 80 | #include <iomanip.h> |
92862013 | 81 | |
6ad0bfa0 | 82 | // --- AliRoot header files --- |
83 | ||
84 | #include "AliRun.h" | |
83448140 | 85 | #include "AliPHOSv1.h" |
6ad0bfa0 | 86 | #include "AliPHOSAnalyze.h" |
6ad0bfa0 | 87 | #include "AliPHOSDigit.h" |
efad3788 | 88 | #include "AliPHOSSDigitizer.h" |
6ad0bfa0 | 89 | #include "AliPHOSTrackSegment.h" |
90 | #include "AliPHOSRecParticle.h" | |
867ede0e | 91 | #include "AliPHOSCpvRecPoint.h" |
83448140 | 92 | #include "AliPHOSPpsdRecPoint.h" |
7b7c1533 | 93 | #include "AliPHOSGetter.h" |
6ad0bfa0 | 94 | |
baef0810 | 95 | |
6ad0bfa0 | 96 | ClassImp(AliPHOSAnalyze) |
97 | ||
6ad0bfa0 | 98 | //____________________________________________________________________________ |
99 | AliPHOSAnalyze::AliPHOSAnalyze() | |
100 | { | |
b2a60966 | 101 | // default ctor (useless) |
efad3788 | 102 | fCorrection = 1.2 ; //Value calculated for default parameters of reconstruction |
6ad0bfa0 | 103 | } |
104 | ||
105 | //____________________________________________________________________________ | |
efad3788 | 106 | AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName) |
6ad0bfa0 | 107 | { |
83974468 | 108 | // ctor: analyze events from root file "name" |
efad3788 | 109 | ffileName = fileName ; |
7b7c1533 | 110 | fCorrection = 1.05 ; //Value calculated for default parameters of reconstruction |
6ad0bfa0 | 111 | } |
112 | ||
88714635 | 113 | //____________________________________________________________________________ |
114 | AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana) | |
115 | { | |
116 | // copy ctor | |
117 | ( (AliPHOSAnalyze &)ana ).Copy(*this) ; | |
118 | } | |
119 | ||
6ad0bfa0 | 120 | //____________________________________________________________________________ |
121 | AliPHOSAnalyze::~AliPHOSAnalyze() | |
122 | { | |
123 | // dtor | |
124 | ||
6ad0bfa0 | 125 | } |
6ad0bfa0 | 126 | //____________________________________________________________________________ |
efad3788 | 127 | void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,const char* branchTitle){ |
128 | //Draws pimary particles and reconstructed | |
129 | //digits, RecPoints, RecPartices etc | |
130 | //for event Nevent in the module Nmod. | |
131 | ||
132 | TH2F * digitOccupancy = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.); | |
133 | TH2F * sdigitOccupancy = new TH2F("sdigitOccupancy","EMC sdigits", 64,-71.,71.,64,-71.,71.); | |
134 | TH2F * emcOccupancy = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.); | |
135 | TH2F * ppsdUp = new TH2F("ppsdUp","PPSD Up digits", 128,-71.,71.,128,-71.,71.) ; | |
136 | TH2F * ppsdUpCl = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ; | |
137 | TH2F * ppsdLow = new TH2F("ppsdLow","PPSD Low digits", 128,-71.,71.,128,-71.,71.) ; | |
138 | TH2F * ppsdLowCl = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ; | |
139 | TH2F * nbar = new TH2F("nbar","Primary nbar", 64,-71.,71.,64,-71.,71.); | |
140 | TH2F * phot = new TH2F("phot","Primary Photon", 64,-71.,71.,64,-71.,71.); | |
141 | TH2F * charg = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.); | |
142 | TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.); | |
143 | TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", 64,-71.,71.,64,-71.,71.); | |
83974468 | 144 | |
7b7c1533 | 145 | //========== Create ObjectGetter |
146 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ; | |
147 | const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; | |
148 | gime->Event(Nevent); | |
ed4205d8 | 149 | |
efad3788 | 150 | //Plot Primary Particles |
7b7c1533 | 151 | const TParticle * primary ; |
efad3788 | 152 | Int_t iPrimary ; |
7b7c1533 | 153 | for ( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++) |
efad3788 | 154 | { |
7b7c1533 | 155 | primary = gime->Primary(iPrimary) ; |
efad3788 | 156 | Int_t primaryType = primary->GetPdgCode() ; |
157 | if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) { | |
158 | Int_t moduleNumber ; | |
159 | Double_t primX, primZ ; | |
7b7c1533 | 160 | phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; |
efad3788 | 161 | if(moduleNumber==Nmod) |
162 | charg->Fill(primZ,primX,primary->Energy()) ; | |
163 | } | |
164 | if( primaryType == 22 ) { | |
165 | Int_t moduleNumber ; | |
166 | Double_t primX, primZ ; | |
7b7c1533 | 167 | phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; |
168 | if(moduleNumber==Nmod) | |
efad3788 | 169 | phot->Fill(primZ,primX,primary->Energy()) ; |
170 | } | |
171 | else{ | |
172 | if( primaryType == -2112 ) { | |
173 | Int_t moduleNumber ; | |
174 | Double_t primX, primZ ; | |
7b7c1533 | 175 | phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; |
efad3788 | 176 | if(moduleNumber==Nmod) |
177 | nbar->Fill(primZ,primX,primary->Energy()) ; | |
178 | } | |
179 | } | |
180 | } | |
ed4205d8 | 181 | |
efad3788 | 182 | Int_t iSDigit ; |
cce01dd5 | 183 | AliPHOSDigit * sdigit ; |
ed4205d8 | 184 | |
cce01dd5 | 185 | TClonesArray * sdigits = gime->SDigits(branchTitle,ffileName) ; |
186 | if(sdigits) | |
187 | for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast(); iSDigit++) | |
efad3788 | 188 | { |
cce01dd5 | 189 | sdigit = (AliPHOSDigit*)sdigits->At(iSDigit) ; |
efad3788 | 190 | Int_t relid[4]; |
7b7c1533 | 191 | phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ; |
efad3788 | 192 | Float_t x,z ; |
7b7c1533 | 193 | phosgeom->RelPosInModule(relid,x,z) ; |
cce01dd5 | 194 | Float_t e = gime->SDigitizer(branchTitle)->Calibrate(sdigit->GetAmp()) ; |
efad3788 | 195 | if(relid[0]==Nmod){ |
196 | if(relid[1]==0) //EMC | |
197 | sdigitOccupancy->Fill(x,z,e) ; | |
198 | if((relid[1]>0)&&(relid[1]<17)) | |
199 | ppsdUp->Fill(x,z,e) ; | |
200 | if(relid[1]>16) | |
201 | ppsdLow->Fill(x,z,e) ; | |
202 | } | |
203 | } | |
ed4205d8 | 204 | |
efad3788 | 205 | //Plot digits |
206 | Int_t iDigit ; | |
cce01dd5 | 207 | AliPHOSDigit * digit ; |
208 | TClonesArray * digits = gime->Digits(branchTitle) ; | |
209 | if(digits) | |
210 | for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++) | |
211 | { | |
212 | digit = (AliPHOSDigit*) digits->At(iDigit) ; | |
213 | Int_t relid[4]; | |
214 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
215 | Float_t x,z ; | |
216 | phosgeom->RelPosInModule(relid,x,z) ; | |
217 | Float_t e = gime->SDigitizer(branchTitle)->Calibrate(digit->GetAmp()) ; | |
218 | if(relid[0]==Nmod){ | |
219 | if(relid[1]==0) //EMC | |
220 | digitOccupancy->Fill(x,z,e) ; | |
221 | if((relid[1]>0)&&(relid[1]<17)) | |
222 | ppsdUp->Fill(x,z,e) ; | |
223 | if(relid[1]>16) | |
224 | ppsdLow->Fill(x,z,e) ; | |
225 | } | |
efad3788 | 226 | } |
cce01dd5 | 227 | |
228 | ||
efad3788 | 229 | //Plot RecPoints |
230 | Int_t irecp ; | |
231 | TVector3 pos ; | |
cce01dd5 | 232 | TObjArray * emcrp = gime->EmcRecPoints(branchTitle) ; |
233 | TObjArray * cpvrp = gime->CpvRecPoints(branchTitle) ; | |
234 | if(cpvrp) | |
235 | for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){ | |
236 | const AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint *) emcrp->At(irecp) ; | |
237 | if(emc->GetPHOSMod()==Nmod){ | |
238 | emc->GetLocalPosition(pos) ; | |
239 | emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy()); | |
240 | } | |
efad3788 | 241 | } |
cce01dd5 | 242 | |
243 | if(cpvrp) | |
244 | for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){ | |
245 | const AliPHOSRecPoint * cpv = (AliPHOSRecPoint *) cpvrp->At(irecp) ; | |
246 | if((strcmp(cpv->ClassName(),"AliPHOSPpsdRecPoint" )) == 0){ // PPSD Rec Point | |
247 | AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) cpv ; | |
efad3788 | 248 | ppsd->GetLocalPosition(pos) ; |
cce01dd5 | 249 | if(ppsd->GetPHOSMod()==Nmod){ |
250 | ppsd->GetLocalPosition(pos) ; | |
251 | if(ppsd->GetUp()) | |
252 | ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy()); | |
253 | else | |
254 | ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy()); | |
255 | } | |
efad3788 | 256 | } |
cce01dd5 | 257 | else{ |
258 | AliPHOSCpvRecPoint * cpv1 = (AliPHOSCpvRecPoint*) cpv ; | |
259 | if(cpv1->GetPHOSMod()==Nmod){ | |
260 | cpv1->GetLocalPosition(pos) ; | |
261 | ppsdUpCl->Fill(pos.X(),pos.Z(),cpv1->GetEnergy()); | |
262 | } | |
efad3788 | 263 | } |
264 | } | |
ed4205d8 | 265 | |
ed4205d8 | 266 | |
efad3788 | 267 | //Plot RecParticles |
7b7c1533 | 268 | const AliPHOSRecParticle * recParticle ; |
efad3788 | 269 | Int_t iRecParticle ; |
cce01dd5 | 270 | TClonesArray * rps = gime->RecParticles(branchTitle) ; |
271 | TClonesArray * tss = gime->TrackSegments(branchTitle) ; | |
272 | if(rps && tss && emcrp) | |
273 | for(iRecParticle = 0; iRecParticle < rps->GetEntriesFast() ;iRecParticle++ ) | |
274 | { | |
275 | recParticle = (AliPHOSRecParticle *) rps->At(iRecParticle) ; | |
276 | Int_t moduleNumberRec ; | |
277 | Double_t recX, recZ ; | |
278 | phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; | |
279 | if(moduleNumberRec == Nmod){ | |
ed4205d8 | 280 | |
cce01dd5 | 281 | Double_t minDistance = 5. ; |
282 | Int_t closestPrimary = -1 ; | |
ad8cfaf4 | 283 | |
cce01dd5 | 284 | |
285 | //extract list of primaries: it is stored at EMC RecPoints | |
286 | Int_t emcIndex = ((AliPHOSTrackSegment *) tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; | |
287 | Int_t numberofprimaries ; | |
288 | Int_t * listofprimaries = ((AliPHOSEmcRecPoint *)emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; | |
289 | Int_t index ; | |
290 | const TParticle * primary ; | |
291 | Double_t distance = minDistance ; | |
292 | ||
293 | for ( index = 0 ; index < numberofprimaries ; index++){ | |
294 | primary = gime->Primary(listofprimaries[index]) ; | |
295 | Int_t moduleNumber ; | |
296 | Double_t primX, primZ ; | |
297 | phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; | |
298 | if(moduleNumberRec == moduleNumber) | |
299 | distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ; | |
300 | if(minDistance > distance) | |
301 | { | |
302 | minDistance = distance ; | |
303 | closestPrimary = listofprimaries[index] ; | |
304 | } | |
305 | } | |
306 | ||
307 | if(closestPrimary >=0 ){ | |
308 | ||
309 | Int_t primaryType = gime->Primary(closestPrimary)->GetPdgCode() ; | |
310 | ||
311 | if(primaryType==22) | |
312 | recPhot->Fill(recZ,recX,recParticle->Energy()) ; | |
313 | else | |
314 | if(primaryType==-2112) | |
315 | recNbar->Fill(recZ,recX,recParticle->Energy()) ; | |
316 | } | |
efad3788 | 317 | } |
318 | } | |
cce01dd5 | 319 | |
ed4205d8 | 320 | |
efad3788 | 321 | //Plot made histograms |
322 | digitOccupancy->Draw("box") ; | |
323 | sdigitOccupancy->SetLineColor(5) ; | |
324 | sdigitOccupancy->Draw("box") ; | |
325 | emcOccupancy->SetLineColor(2) ; | |
326 | emcOccupancy->Draw("boxsame") ; | |
327 | ppsdUp->SetLineColor(3) ; | |
328 | ppsdUp->Draw("boxsame") ; | |
329 | ppsdLow->SetLineColor(4) ; | |
330 | ppsdLow->Draw("boxsame") ; | |
331 | phot->SetLineColor(8) ; | |
332 | phot->Draw("boxsame") ; | |
333 | nbar->SetLineColor(6) ; | |
334 | nbar->Draw("boxsame") ; | |
2bed9e3e | 335 | |
2bed9e3e | 336 | } |
e4761a49 | 337 | //____________________________________________________________________________ |
338 | void AliPHOSAnalyze::Ls(){ | |
339 | //lists branches and titles of PHOS-related branches of TreeR, TreeD, TreeS | |
340 | ||
7b7c1533 | 341 | AliPHOSGetter::GetInstance(ffileName.Data()) ; |
134ce69a | 342 | |
e4761a49 | 343 | Int_t ibranch; |
344 | TObjArray * branches; | |
345 | ||
346 | branches = gAlice->TreeS()->GetListOfBranches() ; | |
347 | ||
348 | cout << "TreeS: " << endl ; | |
349 | for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){ | |
350 | TBranch * branch=(TBranch *) branches->At(ibranch) ; | |
351 | if(strstr(branch->GetName(),"PHOS") ) | |
352 | cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ; | |
353 | } | |
354 | cout << endl ; | |
355 | ||
356 | branches = gAlice->TreeD()->GetListOfBranches() ; | |
357 | ||
358 | cout << "TreeD: " << endl ; | |
359 | for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){ | |
360 | TBranch * branch=(TBranch *) branches->At(ibranch) ; | |
361 | if(strstr(branch->GetName(),"PHOS") ) | |
362 | cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ; | |
363 | } | |
364 | cout << endl ; | |
365 | ||
366 | ||
367 | branches = gAlice->TreeR()->GetListOfBranches() ; | |
368 | ||
369 | cout << "TreeR: " << endl ; | |
370 | for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){ | |
371 | TBranch * branch=(TBranch *) branches->At(ibranch) ; | |
372 | if(strstr(branch->GetName(),"PHOS") ) | |
373 | cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ; | |
374 | } | |
375 | cout << endl ; | |
376 | ||
377 | ||
378 | } | |
2bed9e3e | 379 | //____________________________________________________________________________ |
efad3788 | 380 | void AliPHOSAnalyze::InvariantMass(const char* branchTitle) |
2bed9e3e | 381 | { |
efad3788 | 382 | // Calculates Real and Mixed invariant mass distributions |
7b7c1533 | 383 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ; |
8c0cd6e9 | 384 | |
efad3788 | 385 | Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution |
2bed9e3e | 386 | |
efad3788 | 387 | //opening file |
388 | TFile * mfile = new TFile("invmass.root","update"); | |
69183710 | 389 | |
efad3788 | 390 | //========== Reading /Booking Histograms |
391 | TH2D * hRealEM = 0 ; | |
392 | hRealEM = (TH2D*) mfile->Get("hRealEM") ; | |
393 | if(hRealEM == 0) | |
394 | hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ; | |
395 | TH2D * hRealPhot = 0 ; | |
396 | ||
397 | hRealPhot = (TH2D*)mfile->Get("hRealPhot"); | |
398 | if(hRealPhot == 0) | |
399 | hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ; | |
400 | ||
401 | TH2D * hMixedEM = 0 ; | |
402 | hMixedEM = (TH2D*) mfile->Get("hMixedEM") ; | |
403 | if(hMixedEM == 0) | |
404 | hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ; | |
405 | ||
406 | TH2D * hMixedPhot = 0 ; | |
407 | hMixedPhot = (TH2D*) mfile->Get("hMixedPhot") ; | |
408 | if(hMixedPhot == 0) | |
409 | hMixedPhot = new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ; | |
69183710 | 410 | |
efad3788 | 411 | |
412 | //reading event and copyng it to TConesArray of all photons | |
413 | ||
414 | TClonesArray * allRecParticleList = new TClonesArray("AliPHOSRecParticle", 1000) ; | |
fd0eb6da | 415 | Int_t * nRecParticles = new Int_t[nMixedEvents] ; // to mark boundaries of each event in the total list |
efad3788 | 416 | for(Int_t index = 0; index < nMixedEvents; index ++) |
417 | nRecParticles[index] = 0 ; | |
418 | Int_t iRecPhot = 0 ; // number of EM particles in total list | |
69183710 | 419 | |
efad3788 | 420 | //scan over all events |
421 | Int_t event ; | |
154fc86d | 422 | Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; |
7b7c1533 | 423 | // for(event = 0; event < gime->MaxEvent(); event++ ){ |
424 | for(event = 0; event < maxevent; event++ ){ | |
425 | gime->Event(event); | |
efad3788 | 426 | |
427 | //copy EM RecParticles to the "total" list | |
7b7c1533 | 428 | const AliPHOSRecParticle * recParticle ; |
efad3788 | 429 | Int_t iRecParticle ; |
cce01dd5 | 430 | TClonesArray * rp = gime->RecParticles(branchTitle) ; |
431 | if(rp) | |
432 | for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ) | |
433 | { | |
434 | recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; | |
435 | if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)|| | |
436 | (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)) | |
437 | new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ; | |
438 | } | |
efad3788 | 439 | |
440 | Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle | |
441 | nRecParticles[mevent] = iRecPhot-1 ; | |
ed4205d8 | 442 | |
efad3788 | 443 | //check, if it is time to calculate invariant mass? |
154fc86d | 444 | Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; |
7b7c1533 | 445 | if((mevent == 0) && (event +1 == maxevent)){ |
446 | ||
447 | // if((mevent == 0) && (event +1 == gime->MaxEvent())){ | |
69183710 | 448 | |
efad3788 | 449 | //calculate invariant mass: |
450 | Int_t irp1,irp2 ; | |
451 | Int_t nCurEvent = 0 ; | |
69183710 | 452 | |
efad3788 | 453 | for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){ |
454 | AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ; | |
455 | ||
456 | for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){ | |
457 | AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ; | |
458 | ||
459 | Double_t invMass ; | |
460 | invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())- | |
461 | (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())- | |
462 | (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())- | |
463 | (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ; | |
464 | ||
465 | if(invMass> 0) | |
466 | invMass = TMath::Sqrt(invMass); | |
467 | ||
468 | Double_t pt ; | |
469 | pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) + | |
470 | (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) ); | |
471 | ||
472 | if(irp1 > nRecParticles[nCurEvent]) | |
473 | nCurEvent++; | |
474 | ||
475 | if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event | |
476 | hRealEM->Fill(invMass,pt); | |
477 | if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&& | |
478 | (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) ) | |
479 | hRealPhot->Fill(invMass,pt); | |
480 | } | |
481 | else{ | |
482 | hMixedEM->Fill(invMass,pt); | |
483 | if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&& | |
484 | (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) ) | |
485 | hMixedPhot->Fill(invMass,pt); | |
486 | } //real-mixed | |
487 | ||
488 | } //loop over second rp | |
489 | }//loop over first rp | |
ed4205d8 | 490 | |
efad3788 | 491 | //Make some cleanings |
492 | for(Int_t index = 0; index < nMixedEvents; index ++) | |
493 | nRecParticles[index] = 0 ; | |
494 | iRecPhot = 0 ; | |
495 | allRecParticleList->Clear() ; | |
69183710 | 496 | |
efad3788 | 497 | } |
498 | } | |
499 | delete allRecParticleList ; | |
69183710 | 500 | |
efad3788 | 501 | //writing output |
502 | mfile->cd(); | |
69183710 | 503 | |
efad3788 | 504 | hRealEM->Write(0,kOverwrite) ; |
505 | hRealPhot->Write(0,kOverwrite) ; | |
506 | hMixedEM->Write(0,kOverwrite) ; | |
507 | hMixedPhot->Write(0,kOverwrite) ; | |
69183710 | 508 | |
efad3788 | 509 | mfile->Write(); |
510 | mfile->Close(); | |
511 | delete mfile ; | |
fd0eb6da | 512 | delete nRecParticles; |
c3b9b3f9 | 513 | |
514 | } | |
515 | ||
ed4205d8 | 516 | //____________________________________________________________________________ |
efad3788 | 517 | void AliPHOSAnalyze::EnergyResolution(const char * branchTitle) |
ed4205d8 | 518 | { |
efad3788 | 519 | //fills two dimentional histo: energy of primary vs. energy of reconstructed |
ed4205d8 | 520 | |
efad3788 | 521 | TH2F * hAllEnergy = 0 ; //all reconstructed with primary photon |
522 | TH2F * hPhotEnergy= 0 ; //kGamma with primary photon | |
523 | TH2F * hEMEnergy = 0 ; //electromagnetic with primary photon | |
ed4205d8 | 524 | |
efad3788 | 525 | //opening file and reading histograms if any |
526 | TFile * efile = new TFile("energy.root","update"); | |
ed4205d8 | 527 | |
efad3788 | 528 | hAllEnergy = (TH2F*)efile->Get("hAllEnergy") ; |
529 | if(hAllEnergy == 0) | |
530 | hAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.); | |
531 | ||
532 | hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ; | |
533 | if(hPhotEnergy == 0) | |
534 | hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.); | |
535 | ||
536 | hEMEnergy =(TH2F*) efile->Get("hEMEnergy"); | |
537 | if(hEMEnergy == 0) | |
538 | hEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.); | |
ed4205d8 | 539 | |
ed4205d8 | 540 | |
7b7c1533 | 541 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ; |
542 | const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; | |
efad3788 | 543 | |
544 | Int_t ievent; | |
154fc86d | 545 | Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; |
7b7c1533 | 546 | // for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){ |
547 | for ( ievent=0; ievent < maxevent ; ievent++){ | |
efad3788 | 548 | |
549 | //read the current event | |
7b7c1533 | 550 | gime->Event(ievent) ; |
efad3788 | 551 | |
7b7c1533 | 552 | const AliPHOSRecParticle * recParticle ; |
efad3788 | 553 | Int_t iRecParticle ; |
cce01dd5 | 554 | TClonesArray * rp = gime->RecParticles(branchTitle) ; |
555 | TClonesArray * tss = gime->TrackSegments(branchTitle) ; | |
556 | TObjArray * emcrp = gime->EmcRecPoints(branchTitle) ; | |
557 | if(rp && tss && emcrp) | |
558 | for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ ){ | |
559 | recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; | |
560 | ||
561 | //find the closest primary | |
562 | Int_t moduleNumberRec ; | |
563 | Double_t recX, recZ ; | |
564 | phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; | |
565 | ||
566 | Double_t minDistance = 100. ; | |
567 | Int_t closestPrimary = -1 ; | |
568 | ||
569 | //extract list of primaries: it is stored at EMC RecPoints | |
570 | Int_t emcIndex = ((AliPHOSTrackSegment *)tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; | |
571 | Int_t numberofprimaries ; | |
572 | Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; | |
573 | ||
574 | Int_t index ; | |
575 | const TParticle * primary ; | |
576 | Double_t distance = minDistance ; | |
577 | Double_t dX, dZ; | |
578 | Double_t dXmin = 0.; | |
579 | Double_t dZmin = 0. ; | |
580 | for ( index = 0 ; index < numberofprimaries ; index++){ | |
581 | primary = gime->Primary(listofprimaries[index]) ; | |
582 | Int_t moduleNumber ; | |
583 | Double_t primX, primZ ; | |
584 | phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; | |
585 | if(moduleNumberRec == moduleNumber) { | |
586 | dX = recX - primX; | |
587 | dZ = recZ - primZ; | |
588 | distance = TMath::Sqrt(dX*dX + dZ*dZ) ; | |
589 | if(minDistance > distance) { | |
590 | minDistance = distance ; | |
591 | dXmin = dX; | |
592 | dZmin = dZ; | |
593 | closestPrimary = listofprimaries[index] ; | |
594 | } | |
efad3788 | 595 | } |
596 | } | |
cce01dd5 | 597 | |
598 | //if found primary, fill histograms | |
599 | if(closestPrimary >=0 ){ | |
600 | const TParticle * primary = gime->Primary(closestPrimary) ; | |
601 | if(primary->GetPdgCode() == 22){ | |
602 | hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ; | |
603 | if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){ | |
604 | hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; | |
efad3788 | 605 | hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; |
cce01dd5 | 606 | } |
607 | else | |
608 | if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) | |
609 | hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; | |
610 | } | |
efad3788 | 611 | } |
612 | } | |
efad3788 | 613 | } |
cce01dd5 | 614 | |
efad3788 | 615 | //write filled histograms |
616 | efile->cd() ; | |
617 | hAllEnergy->Write(0,kOverwrite) ; | |
618 | hPhotEnergy->Write(0,kOverwrite) ; | |
619 | hEMEnergy->Write(0,kOverwrite) ; | |
620 | // efile->Write() ; | |
621 | efile->Close() ; | |
622 | delete efile ; | |
cce01dd5 | 623 | |
efad3788 | 624 | } |
fc879520 | 625 | //____________________________________________________________________________ |
efad3788 | 626 | void AliPHOSAnalyze::PositionResolution(const char * branchTitle) |
fc879520 | 627 | { |
efad3788 | 628 | //fills two dimentional histo: energy vs. primary - reconstructed distance |
629 | ||
630 | ||
631 | ||
632 | TH2F * hAllPosition = 0; // Position of any RP with primary photon | |
633 | TH2F * hPhotPosition = 0; // Position of kGAMMA with primary photon | |
634 | TH2F * hEMPosition = 0; // Position of EM with primary photon | |
635 | ||
636 | TH1F * hAllPositionX = 0; // X-Position Resolution of photons with photon primary | |
637 | TH1F * hAllPositionZ = 0; // Z-Position Resolution of photons with photon primary | |
638 | ||
639 | ||
640 | //opening file and reading histograms if any | |
641 | TFile * pfile = new TFile("position.root","update"); | |
642 | ||
643 | hAllPosition = (TH2F*)pfile->Get("hAllPosition"); | |
644 | if(hAllPosition == 0) | |
645 | hAllPosition = new TH2F("hAllPosition", | |
646 | "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.); | |
647 | hPhotPosition= (TH2F*)pfile->Get("hPhotPosition"); | |
648 | if(hPhotPosition == 0) | |
649 | hPhotPosition = new TH2F("hPhotPosition", | |
650 | "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.); | |
651 | hEMPosition= (TH2F*)pfile->Get("hEMPosition") ; | |
652 | if(hEMPosition == 0) | |
653 | hEMPosition = new TH2F("hEMPosition", | |
654 | "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.); | |
655 | hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ; | |
656 | if(hAllPositionX == 0) | |
657 | hAllPositionX = new TH1F("hAllPositionX", | |
658 | "Delta X of any RP with primary photon",100, -2., 2.); | |
659 | hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ; | |
660 | if(hAllPositionZ == 0) | |
661 | hAllPositionZ = new TH1F("hAllPositionZ", | |
662 | "Delta X of any RP with primary photon",100, -2., 2.); | |
663 | ||
664 | ||
7b7c1533 | 665 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ; |
666 | const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; | |
a88eea6f | 667 | |
efad3788 | 668 | Int_t ievent; |
154fc86d | 669 | Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; |
7b7c1533 | 670 | // for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){ |
671 | for ( ievent=0; ievent < maxevent ; ievent++){ | |
672 | ||
efad3788 | 673 | //read the current event |
7b7c1533 | 674 | gime->Event(ievent) ; |
efad3788 | 675 | |
7b7c1533 | 676 | const AliPHOSRecParticle * recParticle ; |
efad3788 | 677 | Int_t iRecParticle ; |
cce01dd5 | 678 | TClonesArray * rp = gime->RecParticles(branchTitle) ; |
679 | TClonesArray * tss= gime->TrackSegments(branchTitle) ; | |
680 | TObjArray * emcrp = gime->EmcRecPoints(branchTitle) ; | |
681 | if( rp && tss && emcrp ) | |
682 | for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){ | |
683 | recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; | |
684 | ||
685 | //find the closest primary | |
686 | Int_t moduleNumberRec ; | |
687 | Double_t recX, recZ ; | |
688 | phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; | |
689 | ||
690 | Double_t minDistance = 100. ; | |
691 | Int_t closestPrimary = -1 ; | |
692 | ||
693 | //extract list of primaries: it is stored at EMC RecPoints | |
694 | Int_t emcIndex = ((AliPHOSTrackSegment *) tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; | |
695 | Int_t numberofprimaries ; | |
696 | Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; | |
697 | ||
698 | Int_t index ; | |
699 | const TParticle * primary ; | |
700 | Double_t distance = minDistance ; | |
701 | Double_t dX = 1000; // incredible number | |
702 | Double_t dZ = 1000; // for the case if no primary will be found | |
703 | Double_t dXmin = 0.; | |
704 | Double_t dZmin = 0. ; | |
705 | for ( index = 0 ; index < numberofprimaries ; index++){ | |
706 | primary = gime->Primary(listofprimaries[index]) ; | |
707 | Int_t moduleNumber ; | |
708 | Double_t primX, primZ ; | |
709 | phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; | |
710 | if(moduleNumberRec == moduleNumber) { | |
711 | dX = recX - primX; | |
712 | dZ = recZ - primZ; | |
713 | distance = TMath::Sqrt(dX*dX + dZ*dZ) ; | |
714 | if(minDistance > distance) { | |
715 | minDistance = distance ; | |
716 | dXmin = dX; | |
717 | dZmin = dZ; | |
718 | closestPrimary = listofprimaries[index] ; | |
719 | } | |
efad3788 | 720 | } |
721 | } | |
cce01dd5 | 722 | |
723 | //if found primary, fill histograms | |
724 | if(closestPrimary >=0 ){ | |
725 | const TParticle * primary = gime->Primary(closestPrimary) ; | |
726 | if(primary->GetPdgCode() == 22){ | |
727 | hAllPosition->Fill(primary->Energy(), minDistance) ; | |
728 | hAllPositionX->Fill(primary->Energy(), dX) ; | |
729 | hAllPositionZ->Fill(primary->Energy(), dZ) ; | |
730 | if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){ | |
731 | hPhotPosition->Fill(primary->Energy(), minDistance ) ; | |
efad3788 | 732 | hEMPosition->Fill(primary->Energy(), minDistance ) ; |
cce01dd5 | 733 | } |
734 | else | |
735 | if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) | |
736 | hEMPosition->Fill(primary->Energy(), minDistance ) ; | |
737 | } | |
efad3788 | 738 | } |
739 | } | |
efad3788 | 740 | } |
741 | ||
742 | //Write output histgrams | |
743 | pfile->cd() ; | |
744 | hAllPosition->Write(0,kOverwrite) ; | |
745 | hAllPositionX->Write(0,kOverwrite) ; | |
746 | hAllPositionZ->Write(0,kOverwrite) ; | |
747 | hPhotPosition->Write(0,kOverwrite) ; | |
748 | hEMPosition->Write(0,kOverwrite) ; | |
749 | pfile->Write() ; | |
750 | pfile->Close() ; | |
751 | delete pfile ; | |
cce01dd5 | 752 | |
753 | ||
efad3788 | 754 | } |
92862013 | 755 | //____________________________________________________________________________ |
efad3788 | 756 | void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ |
757 | // fills spectra of primary photons and several kinds of | |
758 | // reconstructed particles, so that analyzing them one can | |
759 | // estimate conatmination, efficiency of registration etc. | |
760 | ||
761 | //define several general histograms | |
762 | TH1F * hPrimary = 0; //spectrum (P_t distribution) of primary photons | |
763 | TH1F * hAllRP = 0; //spectrum of all RecParticles in PHOS | |
764 | TH1F * hPhot = 0; //spectrum of kGAMMA RecParticles | |
765 | TH1F * hPPSD = 0; //spectrum of all RecParticles with PPSD signal | |
766 | TH1F * hShape = 0; //spectrum of all EM RecParticles | |
767 | TH1F * hVeto = 0; //spectrum of all neutral RecParticles | |
768 | ||
769 | //Now separate histograms in accoradance with primary | |
770 | //primary - photon | |
771 | TH1F * hPhotReg = 0; //Registeres as photon | |
772 | TH1F * hPhotEM = 0; //Registered as EM | |
773 | TH1F * hPhotPPSD= 0; //Registered as RecParticle with PPSD signal | |
774 | ||
775 | //primary - n | |
776 | TH1F * hNReg = 0; //Registeres as photon | |
777 | TH1F * hNEM = 0; //Registered as EM | |
778 | TH1F * hNPPSD= 0; //Registered as RecParticle with PPSD signal | |
779 | ||
780 | //primary - nBar | |
781 | TH1F * hNBarReg = 0; //Registeres as photon | |
782 | TH1F * hNBarEM = 0; //Registered as EM | |
783 | TH1F * hNBarPPSD= 0; //Registered as RecParticle with PPSD signal | |
784 | ||
785 | //primary - charged hadron (pBar excluded) | |
786 | TH1F * hChargedReg = 0; //Registeres as photon | |
787 | TH1F * hChargedEM = 0; //Registered as EM | |
788 | TH1F * hChargedPPSD= 0; //Registered as RecParticle with PPSD signal | |
789 | ||
790 | //primary - pBar | |
791 | TH1F * hPbarReg = 0; //Registeres as photon | |
792 | TH1F * hPbarEM = 0; //Registered as EM | |
793 | TH1F * hPbarPPSD= 0; //Registered as RecParticle with PPSD signal | |
794 | ||
795 | ||
796 | //Reading histograms from the file | |
797 | TFile * cfile = new TFile("contamination.root","update") ; | |
798 | ||
799 | //read general histograms | |
800 | hPrimary = (TH1F*) cfile->Get("hPrimary") ; | |
801 | if(hPrimary == 0) | |
802 | hPrimary= new TH1F("hPrimary", "Primary photon spectrum", 100, 0., 5.); | |
803 | hAllRP = (TH1F*)cfile->Get("hAllRP") ; | |
804 | if(hAllRP == 0) | |
805 | hAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.); | |
806 | hPhot = (TH1F*)cfile->Get("hPhot") ; | |
807 | if(hPhot == 0) | |
808 | hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.); | |
809 | hPPSD = (TH1F*)cfile->Get("hPPSD") ; | |
810 | if(hPPSD == 0) | |
811 | hPPSD = new TH1F("hPPSD", "All PPSD Recparticles", 100, 0., 5.); | |
812 | hShape = (TH1F*) cfile->Get("hShape") ; | |
813 | if(hShape == 0) | |
814 | hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.); | |
815 | hVeto= (TH1F*)cfile->Get("hVeto") ; | |
816 | if(hVeto == 0) | |
817 | hVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.); | |
818 | ||
819 | ||
820 | //primary - photon | |
821 | hPhotReg = (TH1F*)cfile->Get("hPhotReg"); | |
822 | if(hPhotReg == 0) | |
823 | hPhotReg = new TH1F("hPhotReg","Photon registered as photon",100, 0., 5.); | |
824 | hPhotEM =(TH1F*)cfile->Get("hPhotEM"); | |
825 | if(hPhotEM== 0) | |
826 | hPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.); | |
827 | hPhotPPSD= (TH1F*)cfile->Get("hPhotPPSD"); | |
828 | if(hPhotPPSD== 0) | |
829 | hPhotPPSD = new TH1F("hPhotPPSD","Photon registered as PPSD", 100, 0., 5.); | |
830 | ||
831 | //primary - n | |
832 | hNReg = (TH1F*)cfile->Get("hNReg"); | |
833 | if(hNReg== 0) | |
834 | hNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.); | |
835 | hNEM = (TH1F*)cfile->Get("hNEM"); | |
836 | if(hNEM== 0) | |
837 | hNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.); | |
838 | hNPPSD=(TH1F*)cfile->Get("hNPPSD"); | |
839 | if(hNPPSD== 0) | |
840 | hNPPSD = new TH1F("hNPPSD","N registered as PPSD", 100, 0., 5.); | |
841 | ||
842 | //primary - nBar | |
843 | hNBarReg =(TH1F*)cfile->Get("hNBarReg"); | |
844 | if(hNBarReg== 0) | |
845 | hNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.); | |
846 | hNBarEM =(TH1F*)cfile->Get("hNBarEM"); | |
847 | if(hNBarEM== 0) | |
848 | hNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.); | |
849 | hNBarPPSD=(TH1F*)cfile->Get("hNBarPPSD"); | |
850 | if(hNBarPPSD== 0) | |
851 | hNBarPPSD = new TH1F("hNBarPPSD","NBar registered as PPSD", 100, 0., 5.); | |
852 | ||
853 | //primary - charged hadron (pBar excluded) | |
854 | hChargedReg = (TH1F*)cfile->Get("hChargedReg"); | |
855 | if(hChargedReg== 0) | |
856 | hChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.); | |
857 | hChargedEM = (TH1F*)cfile->Get("hChargedEM"); | |
858 | if(hChargedEM== 0) | |
859 | hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.); | |
860 | hChargedPPSD= (TH1F*)cfile->Get("hChargedPPSD"); | |
861 | if(hChargedPPSD== 0) | |
862 | hChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.); | |
863 | ||
864 | //primary - pBar | |
865 | hPbarReg = (TH1F*)cfile->Get("hPbarReg"); | |
866 | if(hPbarReg== 0) | |
867 | hPbarReg= new TH1F("hPbarReg", "pBar registered as photon",100, 0., 5.); | |
868 | hPbarEM = (TH1F*)cfile->Get("hPbarEM"); | |
869 | if(hPbarEM== 0) | |
870 | hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.); | |
871 | hPbarPPSD= (TH1F*)cfile->Get("hPbarPPSD"); | |
872 | if(hPbarPPSD== 0) | |
873 | hPbarPPSD= new TH1F("hPbarPPSD","Pbar as PPSD",100, 0., 5.); | |
eecb6765 | 874 | |
92862013 | 875 | |
efad3788 | 876 | //Now make some initializations |
877 | ||
878 | Int_t counter[8][5] ; //# of registered particles | |
879 | Int_t i1,i2 ; | |
880 | for(i1 = 0; i1<8; i1++) | |
881 | for(i2 = 0; i2<5; i2++) | |
882 | counter[i1][i2] = 0 ; | |
883 | ||
884 | ||
885 | ||
7b7c1533 | 886 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),RecPointsTitle) ; |
887 | const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; | |
154fc86d | 888 | |
efad3788 | 889 | Int_t ievent; |
154fc86d | 890 | Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; |
7b7c1533 | 891 | // for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){ |
892 | for ( ievent=0; ievent < maxevent ; ievent++){ | |
efad3788 | 893 | |
7b7c1533 | 894 | gime->Event(ievent) ; |
efad3788 | 895 | |
896 | //=========== Make spectrum of the primary photons | |
7b7c1533 | 897 | const TParticle * primary ; |
efad3788 | 898 | Int_t iPrimary ; |
7b7c1533 | 899 | for( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++){ |
900 | primary = gime->Primary(iPrimary) ; | |
efad3788 | 901 | Int_t primaryType = primary->GetPdgCode() ; |
902 | if( primaryType == 22 ) { | |
903 | //check, if photons folls onto PHOS | |
904 | Int_t moduleNumber ; | |
905 | Double_t primX, primZ ; | |
7b7c1533 | 906 | phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; |
efad3788 | 907 | if(moduleNumber) |
908 | hPrimary->Fill(primary->Energy()) ; | |
909 | ||
910 | } | |
911 | ||
912 | } | |
913 | //========== Now scan over RecParticles | |
7b7c1533 | 914 | const AliPHOSRecParticle * recParticle ; |
efad3788 | 915 | Int_t iRecParticle ; |
cce01dd5 | 916 | TClonesArray * rp = gime->RecParticles(RecPointsTitle) ; |
917 | TClonesArray * tss= gime->TrackSegments(RecPointsTitle) ; | |
918 | TObjArray * emcrp = gime->EmcRecPoints(RecPointsTitle) ; | |
919 | if( rp && tss && emcrp ) | |
920 | for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){ | |
921 | recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; | |
922 | //fill histo spectrum of all RecParticles | |
923 | hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ; | |
924 | ||
925 | //==========find the closest primary | |
926 | Int_t moduleNumberRec ; | |
927 | Double_t recX, recZ ; | |
928 | phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; | |
929 | ||
930 | Double_t minDistance = 100. ; | |
931 | Int_t closestPrimary = -1 ; | |
932 | ||
933 | //extract list of primaries: it is stored at EMC RecPoints | |
934 | Int_t emcIndex = ((AliPHOSTrackSegment *) tss->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; | |
935 | Int_t numberofprimaries ; | |
936 | Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; | |
937 | Int_t index ; | |
938 | const TParticle * primary ; | |
939 | Double_t distance = minDistance ; | |
940 | Double_t dX, dZ; | |
941 | Double_t dXmin = 0.; | |
942 | Double_t dZmin = 0. ; | |
943 | for ( index = 0 ; index < numberofprimaries ; index++){ | |
944 | primary = gime->Primary(listofprimaries[index]) ; | |
945 | Int_t moduleNumber ; | |
946 | Double_t primX, primZ ; | |
947 | phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; | |
948 | if(moduleNumberRec == moduleNumber) { | |
949 | dX = recX - primX; | |
950 | dZ = recZ - primZ; | |
951 | distance = TMath::Sqrt(dX*dX + dZ*dZ) ; | |
952 | if(minDistance > distance) { | |
953 | minDistance = distance ; | |
954 | dXmin = dX; | |
955 | dZmin = dZ; | |
956 | closestPrimary = listofprimaries[index] ; | |
957 | } | |
efad3788 | 958 | } |
959 | } | |
cce01dd5 | 960 | |
961 | //===========define the "type" of closest primary | |
962 | if(closestPrimary >=0 ){ | |
963 | Int_t primaryCode = -1; | |
964 | const TParticle * primary = gime->Primary(closestPrimary) ; | |
965 | Int_t primaryType = primary->GetPdgCode() ; | |
966 | if(primaryType == 22) // photon ? | |
967 | primaryCode = 0 ; | |
efad3788 | 968 | else |
cce01dd5 | 969 | if(primaryType == 2112) // neutron |
970 | primaryCode = 1 ; | |
efad3788 | 971 | else |
cce01dd5 | 972 | if(primaryType == -2112) // Anti neutron |
973 | primaryCode = 2 ; | |
974 | else | |
975 | if(primaryType == -2122) //Anti proton | |
976 | primaryCode = 4 ; | |
977 | else { | |
978 | TParticle tempo(*primary) ; | |
979 | if(tempo.GetPDG()->Charge()) | |
980 | primaryCode = 3 ; | |
981 | } | |
982 | ||
983 | //==========Now look at the type of RecParticle | |
984 | Float_t energy = CorrectedEnergy(recParticle->Energy()) ; | |
985 | if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){ | |
986 | hPhot->Fill(energy ) ; | |
987 | switch(primaryCode){ | |
988 | case 0: | |
989 | hPhotReg->Fill(energy ) ; | |
990 | break ; | |
991 | case 1: | |
992 | hNReg->Fill(energy ) ; | |
993 | break ; | |
994 | case 2: | |
995 | hNBarReg->Fill(energy ) ; | |
996 | break ; | |
997 | case 3: | |
998 | hChargedReg->Fill(energy ) ; | |
999 | break ; | |
1000 | case 4: | |
1001 | hPbarReg->Fill(energy ) ; | |
1002 | break ; | |
1003 | default: | |
1004 | break ; | |
1005 | } | |
efad3788 | 1006 | } |
cce01dd5 | 1007 | if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)|| |
1008 | (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA)){ //with PPSD signal | |
1009 | hPPSD->Fill(energy ) ; | |
1010 | switch(primaryCode){ | |
1011 | case 0: | |
1012 | hPhotPPSD->Fill(energy ) ; | |
1013 | break ; | |
1014 | case 1: | |
1015 | hNPPSD->Fill(energy ) ; | |
1016 | break ; | |
1017 | case 2: | |
1018 | hNBarPPSD->Fill(energy ) ; | |
1019 | break ; | |
1020 | case 3: | |
1021 | hChargedPPSD->Fill(energy ) ; | |
1022 | break ; | |
1023 | case 4: | |
1024 | hPbarPPSD->Fill(energy ) ; | |
1025 | break ; | |
1026 | default: | |
1027 | break ; | |
1028 | } | |
efad3788 | 1029 | } |
cce01dd5 | 1030 | if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)|| |
1031 | (recParticle->GetType() == AliPHOSFastRecParticle::kELECTRON)|| | |
1032 | (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)|| | |
1033 | (recParticle->GetType() == AliPHOSFastRecParticle::kABSURDEM) ){ //with EM shower | |
1034 | hShape->Fill(energy ) ; | |
1035 | switch(primaryCode){ | |
1036 | case 0: | |
1037 | hPhotEM->Fill(energy ) ; | |
1038 | break ; | |
1039 | case 1: | |
1040 | hNEM->Fill(energy ) ; | |
1041 | break ; | |
1042 | case 2: | |
1043 | hNBarEM->Fill(energy ) ; | |
1044 | break ; | |
1045 | case 3: | |
1046 | hChargedEM->Fill(energy ) ; | |
1047 | break ; | |
1048 | case 4: | |
1049 | hPbarEM->Fill(energy ) ; | |
1050 | break ; | |
1051 | default: | |
1052 | break ; | |
1053 | } | |
efad3788 | 1054 | } |
cce01dd5 | 1055 | |
1056 | if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)|| | |
1057 | (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHA) || | |
1058 | (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA) || | |
1059 | (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) ) //nuetral | |
1060 | hVeto->Fill(energy ) ; | |
1061 | ||
1062 | //fill number of primaries identified as ... | |
1063 | if(primaryCode >= 0) // Primary code defined | |
1064 | counter[recParticle->GetType()][primaryCode]++ ; | |
1065 | ||
efad3788 | 1066 | } |
1067 | ||
cce01dd5 | 1068 | } // no closest primary found |
efad3788 | 1069 | } |
46b146ca | 1070 | |
46b146ca | 1071 | |
efad3788 | 1072 | //=================== SaveHistograms |
1073 | cfile->cd() ; | |
1074 | hPrimary->Write(0,kOverwrite); | |
1075 | hAllRP->Write(0,kOverwrite); | |
1076 | hPhot->Write(0,kOverwrite); | |
1077 | hPPSD->Write(0,kOverwrite); | |
1078 | hShape->Write(0,kOverwrite); | |
1079 | hVeto->Write(0,kOverwrite); | |
1080 | hPhotReg->Write(0,kOverwrite); | |
1081 | hPhotEM->Write(0,kOverwrite); | |
1082 | hPhotPPSD->Write(0,kOverwrite); | |
1083 | hNReg ->Write(0,kOverwrite); | |
1084 | hNEM ->Write(0,kOverwrite); | |
1085 | hNPPSD->Write(0,kOverwrite); | |
1086 | hNBarReg ->Write(0,kOverwrite); | |
1087 | hNBarEM ->Write(0,kOverwrite); | |
1088 | hNBarPPSD->Write(0,kOverwrite); | |
1089 | hChargedReg ->Write(0,kOverwrite); | |
1090 | hChargedEM ->Write(0,kOverwrite); | |
1091 | hChargedPPSD->Write(0,kOverwrite); | |
1092 | hPbarReg ->Write(0,kOverwrite); | |
1093 | hPbarEM ->Write(0,kOverwrite); | |
1094 | hPbarPPSD->Write(0,kOverwrite); | |
46b146ca | 1095 | |
efad3788 | 1096 | cfile->Write(0,kOverwrite); |
1097 | cfile->Close(); | |
1098 | delete cfile ; | |
1099 | ||
69183710 | 1100 | |
efad3788 | 1101 | //print Final Table |
154fc86d | 1102 | maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; |
7b7c1533 | 1103 | // cout << "Resolutions: Analyzed " << gime->MaxEvent() << " event(s)" << endl ; |
2aad621e | 1104 | |
7b7c1533 | 1105 | cout << "Resolutions: Analyzed " << maxevent << " event(s)" << endl ; |
efad3788 | 1106 | cout << endl ; |
b2a60966 | 1107 | |
efad3788 | 1108 | cout << " Primary: Photon Neutron Antineutron Charged hadron AntiProton" << endl ; |
1109 | cout << "--------------------------------------------------------------------------------" << endl ; | |
1110 | cout << " kGAMMA: " | |
1111 | << setw(8) << counter[2][0] << setw(9) << counter[2][1] << setw(13) << counter[2][2] | |
1112 | << setw(15)<< counter[2][3] << setw(13) << counter[2][4] << endl ; | |
1113 | cout << " kGAMMAHA: " | |
1114 | << setw(8) << counter[3][0] << setw(9) << counter[3][1] << setw(13) << counter[3][2] | |
1115 | << setw(15)<< counter[3][3] << setw(13) << counter[3][4] << endl ; | |
1116 | cout << " kNEUTRALEM: " | |
1117 | << setw(8) << counter[0][0] << setw(9) << counter[0][1] << setw(13) << counter[0][2] | |
1118 | << setw(15)<< counter[0][3] << setw(13) << counter[0][4] << endl ; | |
1119 | cout << " kNEUTRALHA: " | |
1120 | << setw(8) << counter[1][0] << setw(9) << counter[1][1] << setw(13) << counter[1][2] | |
1121 | << setw(15)<< counter[1][3] << setw(13) << counter[1][4] << endl ; | |
1122 | cout << " kABSURDEM: " | |
1123 | << setw(8) << counter[4][0] << setw(9) << counter[4][1] << setw(13) << counter[4][2] | |
1124 | << setw(15)<< counter[4][3] << setw(13) << counter[4][4] << endl ; | |
1125 | cout << " kABSURDHA: " | |
1126 | << setw(8) << counter[5][0] << setw(9) << counter[5][1] << setw(13) << counter[5][2] | |
1127 | << setw(15)<< counter[5][3] << setw(13) << counter[5][4] << endl ; | |
1128 | cout << " kELECTRON: " | |
1129 | << setw(8) << counter[6][0] << setw(9) << counter[6][1] << setw(13) << counter[6][2] | |
1130 | << setw(15)<< counter[6][3] << setw(13) << counter[6][4] << endl ; | |
1131 | cout << " kCHARGEDHA: " | |
1132 | << setw(8) << counter[7][0] << setw(9) << counter[7][1] << setw(13) << counter[7][2] | |
1133 | << setw(15)<< counter[7][3] << setw(13) << counter[7][4] << endl ; | |
1134 | cout << "--------------------------------------------------------------------------------" << endl ; | |
1135 | ||
1136 | Int_t totalInd = 0 ; | |
1137 | for(i1 = 0; i1<8; i1++) | |
1138 | for(i2 = 0; i2<5; i2++) | |
1139 | totalInd+=counter[i1][i2] ; | |
1140 | cout << "Indentified particles: " << totalInd << endl ; | |
46b146ca | 1141 | |
69183710 | 1142 | } |
1143 | ||
46b146ca | 1144 | |
efad3788 | 1145 |