1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //_________________________________________________________________________
19 // Algorythm class to analyze PHOSv1 events:
20 // Construct histograms and displays them.
21 // Use the macro EditorBar.C for best access to the functionnalities
23 //*-- Author: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH)
24 //////////////////////////////////////////////////////////////////////////////
26 // --- ROOT system ---
32 #include "TParticle.h"
33 #include "TClonesArray.h"
39 // --- Standard library ---
44 // --- AliRoot header files ---
47 #include "AliPHOSAnalyze.h"
48 #include "AliPHOSClusterizerv1.h"
49 #include "AliPHOSTrackSegmentMakerv1.h"
50 #include "AliPHOSPIDv1.h"
51 #include "AliPHOSReconstructioner.h"
52 #include "AliPHOSDigit.h"
53 #include "AliPHOSTrackSegment.h"
54 #include "AliPHOSRecParticle.h"
55 #include "AliPHOSIndexToObject.h"
56 #include "AliPHOSHit.h"
57 #include "AliPHOSCpvRecPoint.h"
59 ClassImp(AliPHOSAnalyze)
61 //____________________________________________________________________________
62 AliPHOSAnalyze::AliPHOSAnalyze()
64 // default ctor (useless)
69 //____________________________________________________________________________
70 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
72 // ctor: analyze events from root file "name"
74 Bool_t ok = OpenRootFile(name) ;
76 cout << " AliPHOSAnalyze > Error opening " << name << endl ;
79 //========== Get AliRun object from file
80 gAlice = (AliRun*) fRootFile->Get("gAlice") ;
82 //=========== Get the PHOS object and associated geometry from the file
83 fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ;
84 fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
86 //========== Initializes the Index to Object converter
87 fObjGetter = AliPHOSIndexToObject::GetInstance(fPHOS) ;
88 //========== Current event number
100 //____________________________________________________________________________
101 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
104 ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
107 //____________________________________________________________________________
108 void AliPHOSAnalyze::Copy(TObject & obj)
110 // copy an analysis into an other one
112 // I do nothing more because the copy is silly but the Code checkers requires one
115 //____________________________________________________________________________
116 AliPHOSAnalyze::~AliPHOSAnalyze()
120 if(fRootFile->IsOpen()) fRootFile->Close() ;
121 if(fRootFile) {delete fRootFile ; fRootFile=0 ;}
122 if(fPHOS) {delete fPHOS ; fPHOS =0 ;}
123 if(fClu) {delete fClu ; fClu =0 ;}
124 if(fPID) {delete fPID ; fPID =0 ;}
125 if(fRec) {delete fRec ; fRec =0 ;}
126 if(fTrs) {delete fTrs ; fTrs =0 ;}
129 //____________________________________________________________________________
130 void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){
131 //Draws pimary particles and reconstructed
132 //digits, RecPoints, RecPartices etc
133 //for event Nevent in the module Nmod.
135 TH2F * digitOccupancy = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.);
136 TH2F * sdigitOccupancy = new TH2F("sdigitOccupancy","EMC sdigits", 64,-71.,71.,64,-71.,71.);
137 TH2F * emcOccupancy = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.);
138 TH2F * ppsdUp = new TH2F("ppsdUp","PPSD Up digits", 128,-71.,71.,128,-71.,71.) ;
139 TH2F * ppsdUpCl = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ;
140 TH2F * ppsdLow = new TH2F("ppsdLow","PPSD Low digits", 128,-71.,71.,128,-71.,71.) ;
141 TH2F * ppsdLowCl = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ;
142 TH2F * nbar = new TH2F("nbar","Primary nbar", 64,-71.,71.,64,-71.,71.);
143 TH2F * phot = new TH2F("phot","Primary Photon", 64,-71.,71.,64,-71.,71.);
144 TH2F * charg = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.);
145 TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.);
146 TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", 64,-71.,71.,64,-71.,71.);
148 //========== Create the Clusterizer
149 fClu = new AliPHOSClusterizerv1() ;
151 gAlice->GetEvent(Nevent);
153 TParticle * primary ;
155 for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
157 primary = gAlice->Particle(iPrimary) ;
158 Int_t primaryType = primary->GetPdgCode() ;
159 if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) {
161 Double_t primX, primZ ;
162 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
163 if(moduleNumber==Nmod)
164 charg->Fill(primZ,primX,primary->Energy()) ;
166 if( primaryType == 22 ) {
168 Double_t primX, primZ ;
169 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
170 if(moduleNumber==Nmod)
171 phot->Fill(primZ,primX,primary->Energy()) ;
174 if( primaryType == -2112 ) {
176 Double_t primX, primZ ;
177 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
178 if(moduleNumber==Nmod)
179 nbar->Fill(primZ,primX,primary->Energy()) ;
184 // fPHOS->SetTreeAddress() ;
186 gAlice->TreeS()->GetEvent(0) ;
189 AliPHOSDigit * sdigit ;
191 if(fPHOS->SDigits()){
192 for(iSDigit = 0; iSDigit < fPHOS->SDigits()->GetEntries(); iSDigit++)
194 sdigit = (AliPHOSDigit *) fPHOS->SDigits()->At(iSDigit) ;
196 fGeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
198 fGeom->RelPosInModule(relid,x,z) ;
199 Float_t e = fPHOS->Calibrate(sdigit->GetAmp()) ;
201 if(relid[1]==0) //EMC
202 sdigitOccupancy->Fill(x,z,e) ;
203 if((relid[1]>0)&&(relid[1]<17))
204 ppsdUp->Fill(x,z,e) ;
206 ppsdLow->Fill(x,z,e) ;
211 cout << "No SDigits read " << endl ;
214 gAlice->TreeD()->GetEvent(0) ;
218 AliPHOSDigit * digit ;
219 for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++)
221 digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ;
223 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
225 fGeom->RelPosInModule(relid,x,z) ;
226 Float_t e = fClu->Calibrate(digit->GetAmp()) ;
228 if(relid[1]==0) //EMC
229 digitOccupancy->Fill(x,z,e) ;
230 if((relid[1]>0)&&(relid[1]<17))
231 ppsdUp->Fill(x,z,e) ;
233 ppsdLow->Fill(x,z,e) ;
238 cout << "No Digits read " << endl ;
241 gAlice->TreeR()->GetEvent(0) ;
243 TObjArray * emcRecPoints = fPHOS->EmcRecPoints() ;
244 TObjArray * ppsdRecPoints = fPHOS->PpsdRecPoints() ;
245 TClonesArray * recParticleList = fPHOS->RecParticles() ;
252 for(irecp = 0; irecp < emcRecPoints->GetEntries() ; irecp ++){
253 AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)emcRecPoints->At(irecp) ;
254 if(emc->GetPHOSMod()==Nmod){
255 emc->GetLocalPosition(pos) ;
256 emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy());
261 cout << "No EMC rec points read " << endl ;
265 for(irecp = 0; irecp < ppsdRecPoints->GetEntries() ; irecp ++){
266 AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)ppsdRecPoints->At(irecp) ;
267 if(ppsd->GetPHOSMod()==Nmod){
268 ppsd->GetLocalPosition(pos) ;
270 ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
272 ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy());
277 cout << "No PPSD/CPV rec points read " << endl ;
280 AliPHOSRecParticle * recParticle ;
282 if(recParticleList ){
283 for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
285 recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
287 Int_t moduleNumberRec ;
288 Double_t recX, recZ ;
289 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
290 if(moduleNumberRec == Nmod){
292 Double_t minDistance = 5. ;
293 Int_t closestPrimary = -1 ;
295 Int_t numberofprimaries ;
296 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
298 TParticle * primary ;
299 Double_t distance = minDistance ;
301 for ( index = 0 ; index < numberofprimaries ; index++){
302 primary = gAlice->Particle(listofprimaries[index]) ;
304 Double_t primX, primZ ;
305 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
306 if(moduleNumberRec == moduleNumber)
307 distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
308 if(minDistance > distance)
310 minDistance = distance ;
311 closestPrimary = listofprimaries[index] ;
315 if(closestPrimary >=0 ){
317 Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
320 recPhot->Fill(recZ,recX,recParticle->Energy()) ;
322 if(primaryType==-2112)
323 recNbar->Fill(recZ,recX,recParticle->Energy()) ;
329 cout << "Not Rec Prticles read " << endl ;
332 digitOccupancy->Draw("box") ;
333 sdigitOccupancy->SetLineColor(5) ;
334 sdigitOccupancy->Draw("box") ;
335 emcOccupancy->SetLineColor(2) ;
336 emcOccupancy->Draw("boxsame") ;
337 ppsdUp->SetLineColor(3) ;
338 ppsdUp->Draw("boxsame") ;
339 ppsdLow->SetLineColor(4) ;
340 ppsdLow->Draw("boxsame") ;
341 phot->SetLineColor(8) ;
342 phot->Draw("boxsame") ;
343 nbar->SetLineColor(6) ;
344 nbar->Draw("boxsame") ;
347 //____________________________________________________________________________
348 void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent )
351 // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT)
352 // for events from FirstEvent to Nevents
355 for ( ievent=firstEvent; ievent<nevents; ievent++) {
356 if (ievent==firstEvent) {
357 cout << "Analyze > Starting Reconstructing " << endl ;
358 //========== Create the Clusterizer
359 fClu = new AliPHOSClusterizerv1() ;
361 //========== Creates the track segment maker
362 fTrs = new AliPHOSTrackSegmentMakerv1() ;
363 // fTrs->UnsetUnfoldFlag() ;
365 //========== Creates the particle identifier
366 fPID = new AliPHOSPIDv1() ;
367 fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ;
369 //========== Creates the Reconstructioner
370 fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
371 if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE);
374 if (fDebugLevel != 0 ||
375 (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
376 cout << "======= Analyze ======> Event " << ievent+1 << endl ;
379 gAlice->GetEvent(ievent) ;
380 gAlice->SetEvent(ievent) ;
382 if(gAlice->TreeS() == 0) gAlice->MakeTree("S");
383 fPHOS->MakeBranch("S") ;
385 fPHOS->Hits2SDigits() ;
387 if(gAlice->TreeD() == 0) gAlice->MakeTree("D");
388 fPHOS->MakeBranch("D") ;
390 fPHOS->SDigits2Digits() ;
392 if(gAlice->TreeR() == 0) gAlice->MakeTree("R");
394 fPHOS->Reconstruction(fRec);
396 gAlice->TreeS()->Fill() ;
397 gAlice->TreeS()->Write(0,TObject::kOverwrite);
399 gAlice->TreeD()->Fill() ;
400 gAlice->TreeD()->Write(0,TObject::kOverwrite);
404 if(fClu) {delete fClu ; fClu =0 ;}
405 if(fPID) {delete fPID ; fPID =0 ;}
406 if(fRec) {delete fRec ; fRec =0 ;}
407 if(fTrs) {delete fTrs ; fTrs =0 ;}
411 //-------------------------------------------------------------------------------------
412 void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast)
415 // // Read and print generated and reconstructed hits in CPV
416 // // for events from EvFirst to Nevent.
417 // // If only EvFirst is defined, print only this one event.
418 // // Author: Yuri Kharlov
419 // // 12 October 2000
422 // if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
423 // for ( Int_t ievent=EvFirst; ievent<=EvLast; ievent++) {
425 // //========== Event Number>
426 // cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ;
428 // //=========== Connects the various Tree's for evt
429 // Int_t ntracks = gAlice->GetEvent(ievent);
431 // //========== Creating branches ===================================
432 // AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
433 // gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
435 // AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
436 // gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
438 // // Read and print CPV hits
440 // AliPHOSCPVModule cpvModule;
441 // TClonesArray *cpvHits;
443 // AliPHOSCPVHit *cpvHit;
445 // Float_t xgen, zgen;
447 // Int_t nGenHits = 0;
448 // for (Int_t itrack=0; itrack<ntracks; itrack++) {
449 // //=========== Get the Hits Tree for the Primary track itrack
450 // gAlice->ResetHits();
451 // gAlice->TreeH()->GetEvent(itrack);
452 // Int_t iModule = 0 ;
453 // for (iModule=0; iModule < fGeom->GetNCPVModules(); iModule++) {
454 // cpvModule = fPHOS->GetCPVModule(iModule);
455 // cpvHits = cpvModule.Hits();
456 // nCPVhits = cpvHits->GetEntriesFast();
457 // for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
459 // cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
460 // p = cpvHit->GetMomentum();
461 // xgen = cpvHit->X();
462 // zgen = cpvHit->Y();
463 // ipart = cpvHit->GetIpart();
464 // printf("CPV hit in module %d: ",iModule+1);
465 // printf(" p = (%f, %f, %f, %f) GeV,\n",
466 // p.Px(),p.Py(),p.Pz(),p.Energy());
467 // printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d\n",
473 // // Read and print CPV reconstructed points
475 // //=========== Gets the Reconstruction TTree
476 // gAlice->TreeR()->GetEvent(0) ;
477 // printf("Recpoints: %d\n",(*fPHOS->CpvRecPoints())->GetEntries());
478 // TIter nextRP(*fPHOS->CpvRecPoints() ) ;
479 // AliPHOSCpvRecPoint *cpvRecPoint ;
480 // Int_t nRecPoints = 0;
481 // while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
484 // cpvRecPoint->GetLocalPosition(locpos);
485 // Int_t phosModule = cpvRecPoint->GetPHOSMod();
486 // printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cm\n",
487 // phosModule,locpos.X(),locpos.Z());
489 // printf("This event has %d generated hits and %d reconstructed points\n",
490 // nGenHits,nRecPoints);
494 //____________________________________________________________________________
495 void AliPHOSAnalyze::AnalyzeCPV(Int_t Nevents)
498 // // Analyzes CPV characteristics
499 // // Author: Yuri Kharlov
503 // // Book histograms
505 // TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.);
506 // TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.);
507 // TH1F *hDr = new TH1F("hDr" ,"CPV r-resolution@reconstruction",100, 0. , 5.);
508 // TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5);
509 // TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" , 21,-0.5,20.5);
510 // TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" , 21,-0.5,20.5);
512 // cout << "Start CPV Analysis"<< endl ;
513 // for ( Int_t ievent=0; ievent<Nevents; ievent++) {
515 // //========== Event Number>
516 // // if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
517 // cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ;
519 // //=========== Connects the various Tree's for evt
520 // Int_t ntracks = gAlice->GetEvent(ievent);
522 // //========== Creating branches ===================================
523 // AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ;
524 // gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ;
526 // AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ;
527 // gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ;
529 // // Create and fill arrays of hits for each CPV module
531 // Int_t nOfModules = fGeom->GetNModules();
532 // TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
533 // Int_t iModule = 0;
534 // for (iModule=0; iModule < nOfModules; iModule++)
535 // hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
537 // AliPHOSCPVModule cpvModule;
538 // TClonesArray *cpvHits;
540 // AliPHOSCPVHit *cpvHit;
545 // // First go through all primary tracks and fill the arrays
546 // // of hits per each CPV module
548 // for (Int_t itrack=0; itrack<ntracks; itrack++) {
549 // // Get the Hits Tree for the Primary track itrack
550 // gAlice->ResetHits();
551 // gAlice->TreeH()->GetEvent(itrack);
552 // for (Int_t iModule=0; iModule < nOfModules; iModule++) {
553 // cpvModule = fPHOS->GetCPVModule(iModule);
554 // cpvHits = cpvModule.Hits();
555 // nCPVhits = cpvHits->GetEntriesFast();
556 // for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
557 // cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
558 // p = cpvHit->GetMomentum();
559 // xzgen[0] = cpvHit->X();
560 // xzgen[1] = cpvHit->Y();
561 // ipart = cpvHit->GetIpart();
562 // TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
563 // new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit);
565 // cpvModule.Clear();
568 // for (iModule=0; iModule < nOfModules; iModule++) {
569 // Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
570 // printf("Module %d has %d hits\n",iModule,nsum);
573 // // Then go through reconstructed points and for each find
574 // // the closeset hit
575 // // The distance from the rec.point to the closest hit
576 // // gives the coordinate resolution of the CPV
578 // // Get the Reconstruction Tree
579 // gAlice->TreeR()->GetEvent(0) ;
580 // TIter nextRP(*fPHOS->PpsdRecPoints() ) ;
581 // AliPHOSCpvRecPoint *cpvRecPoint ;
582 // Float_t xgen, zgen;
583 // while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
585 // cpvRecPoint->GetLocalPosition(locpos);
586 // Int_t phosModule = cpvRecPoint->GetPHOSMod();
587 // Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity();
588 // Int_t rpMultX, rpMultZ;
589 // cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
590 // Float_t xrec = locpos.X();
591 // Float_t zrec = locpos.Z();
592 // Float_t dxmin = 1.e+10;
593 // Float_t dzmin = 1.e+10;
594 // Float_t r2min = 1.e+10;
597 // cpvHits = hitsPerModule[phosModule-1];
598 // Int_t nCPVhits = cpvHits->GetEntriesFast();
599 // for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
600 // cpvHit = (AliPHOSCPVHit*)cpvHits->UncheckedAt(ihit);
601 // xgen = cpvHit->X();
602 // zgen = cpvHit->Y();
603 // r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2);
604 // if ( r2 < r2min ) {
606 // dxmin = xgen - xrec;
607 // dzmin = zgen - zrec;
610 // hDx ->Fill(dxmin);
611 // hDz ->Fill(dzmin);
612 // hDr ->Fill(TMath::Sqrt(r2min));
613 // hNrp ->Fill(rpMult);
614 // hNrpX->Fill(rpMultX);
615 // hNrpZ->Fill(rpMultZ);
617 // delete [] hitsPerModule;
619 // // Save histograms
621 // Text_t outputname[80] ;
622 // sprintf(outputname,"%s.analyzed",fRootFile->GetName());
623 // TFile output(outputname,"RECREATE");
633 // // Plot histograms
635 // TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400);
636 // gStyle->SetOptStat(111111);
637 // gStyle->SetOptFit(1);
638 // gStyle->SetOptDate(1);
639 // cpvCanvas->Divide(3,2);
642 // gPad->SetFillColor(10);
643 // hNrp->SetFillColor(16);
647 // gPad->SetFillColor(10);
648 // hNrpX->SetFillColor(16);
652 // gPad->SetFillColor(10);
653 // hNrpZ->SetFillColor(16);
657 // gPad->SetFillColor(10);
658 // hDx->SetFillColor(16);
663 // gPad->SetFillColor(10);
664 // hDz->SetFillColor(16);
669 // gPad->SetFillColor(10);
670 // hDr->SetFillColor(16);
673 // cpvCanvas->Print("CPV.ps");
677 //____________________________________________________________________________
678 void AliPHOSAnalyze::InvariantMass(Int_t Nevents )
680 // Calculates Real and Mixed invariant mass distributions
682 const Int_t knMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution
683 Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/knMixedEvents) ;
685 //========== Booking Histograms
686 TH2D * hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ;
687 TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
688 TH2D * hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ;
689 TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
692 Int_t eventInMixedLoop ;
694 Int_t nRecParticles[4];//knMixedEvents] ;
696 AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", knMixedEvents*1000) ;
698 for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++ ){
701 for ( ievent=0; ievent < knMixedEvents; ievent++){
703 Int_t absEventNumber = eventInMixedLoop*knMixedEvents + ievent ;
705 //=========== Connects the various Tree's for evt
706 gAlice->GetEvent(absEventNumber);
708 //========== Creating branches ===================================
709 fPHOS->SetTreeAddress() ;
711 gAlice->TreeD()->GetEvent(0) ;
712 gAlice->TreeR()->GetEvent(0) ;
714 TClonesArray * recParticleList = fPHOS->RecParticles() ;
717 AliPHOSRecParticle * recParticle ;
719 for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
721 recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
722 if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)||
723 (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){
724 new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ;
729 nRecParticles[ievent] = iRecPhot-1 ;
732 //Now calculate invariant mass:
734 Int_t nCurEvent = 0 ;
736 for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
737 AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
739 for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
740 AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
743 invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
744 (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
745 (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
746 (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
749 invMass = TMath::Sqrt(invMass);
752 pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) +(rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py()));
754 if(irp1 > nRecParticles[nCurEvent])
757 if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
758 hRealEM->Fill(invMass,pt);
759 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
760 hRealPhot->Fill(invMass,pt);
763 hMixedEM->Fill(invMass,pt);
764 if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA))
765 hMixedPhot->Fill(invMass,pt);
768 } //loop over second rp
769 }//loop over first rp
770 allRecParticleList->Delete() ;
773 delete allRecParticleList ;
776 TFile output("invmass.root","RECREATE");
782 hMixedPhot->Write() ;
789 //____________________________________________________________________________
790 void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast)
793 // // Read and print generated and reconstructed hits in EMC
794 // // for events from EvFirst to Nevent.
795 // // If only EvFirst is defined, print only this one event.
796 // // Author: Yuri Kharlov
797 // // 24 November 2000
800 // if (EvFirst!=0 && EvLast==0) EvLast=EvFirst;
802 // for (ievent=EvFirst; ievent<=EvLast; ievent++) {
804 // //========== Event Number>
805 // cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ;
807 // //=========== Connects the various Tree's for evt
808 // Int_t ntracks = gAlice->GetEvent(ievent);
809 // fPHOS->SetTreeAddress() ;
811 // gAlice->TreeD()->GetEvent(0) ;
812 // gAlice->TreeR()->GetEvent(0) ;
814 // // Loop over reconstructed particles
816 // TClonesArray ** recParticleList = fPHOS->RecParticles() ;
817 // AliPHOSRecParticle * recParticle ;
818 // Int_t iRecParticle ;
821 // for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) {
822 // recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ;
823 // Float_t recE = recParticle->Energy();
824 // primList = recParticle->GetPrimaries(nPrimary);
825 // Int_t moduleNumberRec ;
826 // Double_t recX, recZ ;
827 // fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
828 // printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n",
829 // moduleNumberRec,recX,recZ,recE,*primList);
832 // // Read and print EMC hits from EMCn branches
834 // AliPHOSCPVModule emcModule;
835 // TClonesArray *emcHits;
837 // AliPHOSCPVHit *emcHit;
839 // Float_t xgen, zgen;
840 // Int_t ipart, primary;
841 // Int_t nGenHits = 0;
842 // for (Int_t itrack=0; itrack<ntracks; itrack++) {
843 // //=========== Get the Hits Tree for the Primary track itrack
844 // gAlice->ResetHits();
845 // gAlice->TreeH()->GetEvent(itrack);
846 // Int_t iModule = 0 ;
847 // for (iModule=0; iModule < fGeom->GetNModules(); iModule++) {
848 // emcModule = fPHOS->GetEMCModule(iModule);
849 // emcHits = emcModule.Hits();
850 // nEMChits = emcHits->GetEntriesFast();
851 // for (Int_t ihit=0; ihit<nEMChits; ihit++) {
853 // emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
854 // p = emcHit->GetMomentum();
855 // xgen = emcHit->X();
856 // zgen = emcHit->Y();
857 // ipart = emcHit->GetIpart();
858 // primary= emcHit->GetTrack();
859 // printf("EMC hit A: module %d, ",iModule+1);
860 // printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n",
861 // p.Px(),p.Py(),p.Pz(),p.Energy());
862 // printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
863 // xgen,zgen,ipart,primary);
868 // // // Read and print EMC hits from PHOS branch
870 // // for (Int_t itrack=0; itrack<ntracks; itrack++) {
871 // // //=========== Get the Hits Tree for the Primary track itrack
872 // // gAlice->ResetHits();
873 // // gAlice->TreeH()->GetEvent(itrack);
874 // // TClonesArray *hits = fPHOS->Hits();
875 // // AliPHOSHit *hit ;
877 // // for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) {
878 // // hit = (AliPHOSHit*)hits->At(ihit) ;
879 // // Float_t hitXYZ[3];
880 // // hitXYZ[0] = hit->X();
881 // // hitXYZ[1] = hit->Y();
882 // // hitXYZ[2] = hit->Z();
883 // // ipart = hit->GetPid();
884 // // primary = hit->GetPrimary();
885 // // Int_t absId = hit->GetId();
886 // // Int_t relId[4];
887 // // fGeom->AbsToRelNumbering(absId, relId) ;
888 // // Int_t module = relId[0];
889 // // if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0))
890 // // printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
891 // // module,hitXYZ[0],hitXYZ[2],ipart,primary);
898 //____________________________________________________________________________
899 void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents)
902 // // Read generated and reconstructed hits in EMC for Nevents events.
903 // // Plots the coordinate and energy resolution histograms.
904 // // Coordinate resolution is a difference between the reconstructed
905 // // coordinate and the exact coordinate on the face of the PHOS
906 // // Author: Yuri Kharlov
907 // // 27 November 2000
910 // // Book histograms
912 // TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.);
913 // TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.);
914 // TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.);
916 // TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.);
917 // TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.);
918 // TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.);
920 // cout << "Start EMC Analysis"<< endl ;
921 // for (Int_t ievent=0; ievent<Nevents; ievent++) {
923 // //========== Event Number>
924 // if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0)
925 // cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ;
927 // //=========== Connects the various Tree's for evt
928 // Int_t ntracks = gAlice->GetEvent(ievent);
930 // fPHOS->SetTreeAddress() ;
932 // gAlice->TreeD()->GetEvent(0) ;
933 // gAlice->TreeR()->GetEvent(0) ;
935 // // Create and fill arrays of hits for each EMC module
937 // Int_t nOfModules = fGeom->GetNModules();
938 // TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
940 // for (iModule=0; iModule < nOfModules; iModule++)
941 // hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100);
943 // AliPHOSCPVModule emcModule;
944 // TClonesArray *emcHits;
946 // AliPHOSCPVHit *emcHit;
948 // // First go through all primary tracks and fill the arrays
949 // // of hits per each EMC module
951 // for (Int_t itrack=0; itrack<ntracks; itrack++) {
952 // // Get the Hits Tree for the Primary track itrack
953 // gAlice->ResetHits();
954 // gAlice->TreeH()->GetEvent(itrack);
955 // for (Int_t iModule=0; iModule < nOfModules; iModule++) {
956 // emcModule = fPHOS->GetEMCModule(iModule);
957 // emcHits = emcModule.Hits();
958 // nEMChits = emcHits->GetEntriesFast();
959 // for (Int_t ihit=0; ihit<nEMChits; ihit++) {
960 // emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
961 // TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
962 // new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit);
964 // emcModule.Clear();
968 // // Loop over reconstructed particles
970 // TClonesArray ** recParticleList = fPHOS->RecParticles() ;
971 // AliPHOSRecParticle * recParticle ;
972 // Int_t nEMCrecs = (*recParticleList)->GetEntries();
973 // if (nEMCrecs == 1) {
974 // recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ;
975 // Float_t recE = recParticle->Energy();
977 // Double_t recX, recZ ;
978 // fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ;
980 // // for this rec.point take the hit list in the same PHOS module
982 // emcHits = hitsPerModule[phosModule-1];
983 // Int_t nEMChits = emcHits->GetEntriesFast();
984 // if (nEMChits == 1) {
985 // Float_t genX, genZ, genE;
986 // for (Int_t ihit=0; ihit<nEMChits; ihit++) {
987 // emcHit = (AliPHOSCPVHit*)emcHits->UncheckedAt(ihit);
988 // genX = emcHit->X();
989 // genZ = emcHit->Y();
990 // genE = emcHit->GetMomentum().E();
992 // Float_t dx = recX - genX;
993 // Float_t dz = recZ - genZ;
994 // Float_t de = recE - genE;
998 // hDx2 ->Fill(genE,dx);
999 // hDz2 ->Fill(genE,dz);
1000 // hDE2 ->Fill(genE,recE);
1003 // delete [] hitsPerModule;
1005 // // Save histograms
1007 // Text_t outputname[80] ;
1008 // sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1009 // TFile output(outputname,"RECREATE");
1019 // // Plot histograms
1021 // TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300);
1022 // gStyle->SetOptStat(111111);
1023 // gStyle->SetOptFit(1);
1024 // gStyle->SetOptDate(1);
1025 // emcCanvas->Divide(3,1);
1027 // emcCanvas->cd(1);
1028 // gPad->SetFillColor(10);
1029 // hDx1->SetFillColor(16);
1032 // emcCanvas->cd(2);
1033 // gPad->SetFillColor(10);
1034 // hDz1->SetFillColor(16);
1037 // emcCanvas->cd(3);
1038 // gPad->SetFillColor(10);
1039 // hDE1->SetFillColor(16);
1042 // emcCanvas->Print("EMC.ps");
1046 //____________________________________________________________________________
1047 void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents )
1049 // analyzes Nevents events and calculate Energy and Position resolution as well as
1050 // probaility of correct indentifiing of the incident particle
1052 //========== Booking Histograms
1053 cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ;
1054 BookResolutionHistograms();
1056 Int_t counter[9][5] ;
1057 Int_t i1,i2,totalInd = 0 ;
1058 for(i1 = 0; i1<9; i1++)
1059 for(i2 = 0; i2<5; i2++)
1060 counter[i1][i2] = 0 ;
1062 Int_t totalPrimary = 0 ;
1063 Int_t totalRecPart = 0 ;
1064 Int_t totalRPwithPrim = 0 ;
1067 cout << "Start Analysing"<< endl ;
1068 for ( ievent=0; ievent<Nevents; ievent++)
1071 //========== Event Number>
1072 // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. )
1073 cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ;
1075 //=========== Connects the various Tree's for evt
1076 gAlice->GetEvent(ievent);
1078 //=========== Gets the Kine TTree
1079 gAlice->TreeK()->GetEvent(0) ;
1081 //=========== Gets the list of Primari Particles
1083 TParticle * primary ;
1085 for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++)
1087 primary = gAlice->Particle(iPrimary) ;
1088 Int_t primaryType = primary->GetPdgCode() ;
1089 if( primaryType == 22 ) {
1090 Int_t moduleNumber ;
1091 Double_t primX, primZ ;
1092 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1094 fhPrimary->Fill(primary->Energy()) ;
1095 if(primary->Energy() > 0.3)
1101 fPHOS->SetTreeAddress() ;
1103 gAlice->TreeD()->GetEvent(0) ;
1104 gAlice->TreeR()->GetEvent(0) ;
1106 TClonesArray * recParticleList = fPHOS->RecParticles() ;
1108 AliPHOSRecParticle * recParticle ;
1109 Int_t iRecParticle ;
1110 for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ )
1112 recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ;
1113 fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ;
1115 Int_t moduleNumberRec ;
1116 Double_t recX, recZ ;
1117 fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
1119 Double_t minDistance = 100. ;
1120 Int_t closestPrimary = -1 ;
1122 Int_t numberofprimaries ;
1123 Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ;
1125 TParticle * primary ;
1126 Double_t distance = minDistance ;
1128 Double_t dXmin = 0.;
1129 Double_t dZmin = 0. ;
1130 for ( index = 0 ; index < numberofprimaries ; index++){
1131 primary = gAlice->Particle(listofprimaries[index]) ;
1132 Int_t moduleNumber ;
1133 Double_t primX, primZ ;
1134 fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
1135 if(moduleNumberRec == moduleNumber) {
1138 distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
1139 if(minDistance > distance) {
1140 minDistance = distance ;
1143 closestPrimary = listofprimaries[index] ;
1149 if(closestPrimary >=0 ){
1152 Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ;
1153 // TParticlePDG* pDGparticle = gAlice->ParticleAt(closestPrimary)->GetPDG();
1154 // Double_t charge = PDGparticle->Charge() ;
1156 // cout <<"Primary " <<primaryType << " E " << ((TParticle *)primaryList->At(closestPrimary))->Energy() << endl ;
1161 primaryCode = 0; //Photon
1162 fhAllEnergy ->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy()) ;
1163 fhAllPosition ->Fill(gAlice->Particle(closestPrimary)->Energy(), minDistance) ;
1164 fhAllPositionX->Fill(dXmin);
1165 fhAllPositionZ->Fill(dZmin);
1168 primaryCode = 1; //Electron
1171 primaryCode = 1; //positron
1174 primaryCode = 4; //K+
1177 primaryCode = 4; //K-
1180 primaryCode = 4; //K0s
1183 primaryCode = 4; //K0l
1186 primaryCode = 2; //K0l
1189 primaryCode = 2; //K0l
1192 primaryCode = 2; //K0l
1195 primaryCode = 2; //K0l
1198 primaryCode = 3; //ELSE
1202 switch(recParticle->GetType())
1204 case AliPHOSFastRecParticle::kGAMMA:
1205 if(primaryType == 22){
1206 fhPhotEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1207 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1208 fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1210 fhPhotPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1211 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1212 fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1214 fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1215 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1216 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1218 fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1220 if(primaryType == 2112){ //neutron
1221 fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1222 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1223 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1226 if(primaryType == -2112){ //neutron ~
1227 fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1228 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1229 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1232 if(primaryCode == 2){
1233 fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1234 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1235 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1238 fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1239 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1240 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1241 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1242 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1243 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1244 counter[0][primaryCode]++;
1246 case AliPHOSFastRecParticle::kELECTRON:
1247 if(primaryType == 22){
1248 fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1249 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1250 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1251 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1252 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1254 if(primaryType == 2112){ //neutron
1255 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1256 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1259 if(primaryType == -2112){ //neutron ~
1260 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1261 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1264 if(primaryCode == 2){
1265 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1266 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1269 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1270 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1271 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1272 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1273 counter[1][primaryCode]++;
1275 case AliPHOSFastRecParticle::kNEUTRALHA:
1276 if(primaryType == 22)
1277 fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1279 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1280 counter[2][primaryCode]++;
1282 case AliPHOSFastRecParticle::kNEUTRALEM:
1283 if(primaryType == 22){
1284 fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(),recParticle->Energy() ) ;
1285 fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance ) ;
1287 fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1288 fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1290 if(primaryType == 2112) //neutron
1291 fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1293 if(primaryType == -2112) //neutron ~
1294 fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1296 if(primaryCode == 2)
1297 fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1299 fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1300 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1301 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1303 counter[3][primaryCode]++;
1305 case AliPHOSFastRecParticle::kCHARGEDHA:
1306 if(primaryType == 22) //photon
1307 fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1309 counter[4][primaryCode]++ ;
1311 case AliPHOSFastRecParticle::kGAMMAHA:
1312 if(primaryType == 22){ //photon
1313 fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1314 fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ;
1315 fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ;
1316 fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1318 if(primaryType == 2112){ //neutron
1319 fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1322 if(primaryType == -2112){ //neutron ~
1323 fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1325 if(primaryCode == 2){
1326 fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1329 fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1330 fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1331 fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1332 counter[5][primaryCode]++ ;
1334 case AliPHOSFastRecParticle::kABSURDEM:
1335 counter[6][primaryCode]++ ;
1336 fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ;
1338 case AliPHOSFastRecParticle::kABSURDHA:
1339 counter[7][primaryCode]++ ;
1342 counter[8][primaryCode]++ ;
1349 cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ;
1350 cout << "Resolutions: Total primary " << totalPrimary << endl ;
1351 cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ;
1352 cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ;
1353 cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ;
1354 cout << " Detected as photon " << counter[0][0] << " " << counter[0][1] << " " << counter[0][2] << " " <<counter[0][3] << " " << counter[0][4] << endl ;
1355 cout << " Detected as electron " << counter[1][0] << " " << counter[1][1] << " " << counter[1][2] << " " <<counter[1][3] << " " << counter[1][4] << endl ;
1356 cout << " Detected as neutral hadron " << counter[2][0] << " " << counter[2][1] << " " << counter[2][2] << " " <<counter[2][3] << " " << counter[2][4] << endl ;
1357 cout << " Detected as neutral EM " << counter[3][0] << " " << counter[3][1] << " " << counter[3][2] << " " <<counter[3][3] << " " << counter[3][4] << endl ;
1358 cout << " Detected as charged hadron " << counter[4][0] << " " << counter[4][1] << " " << counter[4][2] << " " <<counter[4][3] << " " << counter[4][4] << endl ;
1359 cout << " Detected as gamma-hadron " << counter[5][0] << " " << counter[5][1] << " " << counter[5][2] << " " <<counter[5][3] << " " << counter[5][4] << endl ;
1360 cout << " Detected as Absurd EM " << counter[6][0] << " " << counter[6][1] << " " << counter[6][2] << " " <<counter[6][3] << " " << counter[6][4] << endl ;
1361 cout << " Detected as absurd hadron " << counter[7][0] << " " << counter[7][1] << " " << counter[7][2] << " " <<counter[7][3] << " " << counter[7][4] << endl ;
1362 cout << " Detected as undefined " << counter[8][0] << " " << counter[8][1] << " " << counter[8][2] << " " <<counter[8][3] << " " << counter[8][4] << endl ;
1364 for(i1 = 0; i1<9; i1++)
1365 for(i2 = 0; i2<5; i2++)
1366 totalInd+=counter[i1][i2] ;
1367 cout << "Indentified particles " << totalInd << endl ;
1372 //____________________________________________________________________________
1373 void AliPHOSAnalyze::BookingHistograms()
1375 // Books the histograms where the results of the analysis are stored (to be changed)
1378 delete fhVetoDigit ;
1379 delete fhConvertorDigit ;
1380 delete fhEmcCluster ;
1381 delete fhVetoCluster ;
1382 delete fhConvertorCluster ;
1383 delete fhConvertorEmc ;
1385 fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.);
1386 fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5);
1387 fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5);
1388 fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.);
1389 fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5);
1390 fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5);
1391 fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5);
1394 //____________________________________________________________________________
1395 void AliPHOSAnalyze::BookResolutionHistograms()
1397 // Books the histograms where the results of the Resolution analysis are stored
1400 // delete fhAllEnergy ;
1402 // delete fhPhotEnergy ;
1404 // delete fhEMEnergy ;
1406 // delete fhPPSDEnergy ;
1409 fhAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1410 fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1411 fhEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1412 fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1414 // if(fhAllPosition)
1415 // delete fhAllPosition ;
1416 // if(fhPhotPosition)
1417 // delete fhPhotPosition ;
1419 // delete fhEMPosition ;
1420 // if(fhPPSDPosition)
1421 // delete fhPPSDPosition ;
1424 fhAllPosition = new TH2F("hAllPosition", "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
1425 fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
1426 fhEMPosition = new TH2F("hEMPosition", "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.);
1427 fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon", 100, 0., 5., 100, 0., 5.);
1429 fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.);
1430 fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.);
1433 // delete fhAllReg ;
1435 // delete fhPhotReg ;
1439 // delete fhNBarReg ;
1441 // delete fhChargedReg ;
1443 fhAllReg = new TH1F("hAllReg", "All primaries registered as photon", 100, 0., 5.);
1444 fhPhotReg = new TH1F("hPhotReg", "Photon registered as photon", 100, 0., 5.);
1445 fhNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.);
1446 fhNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.);
1447 fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
1452 // delete fhPhotEM ;
1456 // delete fhNBarEM ;
1458 // delete fhChargedEM ;
1460 fhAllEM = new TH1F("hAllEM", "All primary registered as EM",100, 0., 5.);
1461 fhPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.);
1462 fhNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.);
1463 fhNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.);
1464 fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
1467 // delete fhAllPPSD ;
1469 // delete fhPhotPPSD ;
1473 // delete fhNBarPPSD ;
1474 // if(fhChargedPPSD)
1475 // delete fhChargedPPSD ;
1477 fhAllPPSD = new TH1F("hAllPPSD", "All primary registered as PPSD",100, 0., 5.);
1478 fhPhotPPSD = new TH1F("hPhotPPSD", "Photon registered as PPSD", 100, 0., 5.);
1479 fhNPPSD = new TH1F("hNPPSD", "N registered as PPSD", 100, 0., 5.);
1480 fhNBarPPSD = new TH1F("hNBarPPSD", "NBar registered as PPSD", 100, 0., 5.);
1481 fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.);
1484 // delete fhPrimary ;
1485 fhPrimary= new TH1F("hPrimary", "hPrimary", 100, 0., 5.);
1496 fhAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
1497 fhVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.);
1498 fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.);
1499 fhPPSD = new TH1F("hPPSD", "All PPSD photon particles", 100, 0., 5.);
1503 // delete fhPhotPhot ;
1505 // delete fhPhotElec ;
1507 // delete fhPhotNeuH ;
1509 // delete fhPhotNuEM ;
1511 // delete fhPhotChHa ;
1513 // delete fhPhotGaHa ;
1515 fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.); //Photon registered as photon
1516 fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.); //Photon registered as Electron
1517 fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.); //Photon registered as Neutral Hadron
1518 fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.); //Photon registered as Neutral EM
1519 fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.); //Photon registered as Charged Hadron
1520 fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.); //Photon registered as Gamma-Hadron
1523 //____________________________________________________________________________
1524 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
1526 // Open the root file named "name"
1528 fRootFile = new TFile(name, "update") ;
1529 return fRootFile->IsOpen() ;
1532 //____________________________________________________________________________
1533 void AliPHOSAnalyze::SaveHistograms()
1535 // Saves the histograms in a root file named "name.analyzed"
1537 Text_t outputname[80] ;
1538 sprintf(outputname,"%s.analyzed",fRootFile->GetName());
1539 TFile output(outputname,"RECREATE");
1543 fhAllEnergy->Write() ;
1545 fhPhotEnergy->Write() ;
1547 fhEMEnergy->Write() ;
1549 fhPPSDEnergy->Write() ;
1551 fhAllPosition->Write() ;
1553 fhAllPositionX->Write() ;
1555 fhAllPositionZ->Write() ;
1557 fhPhotPosition->Write() ;
1559 fhEMPosition->Write() ;
1561 fhPPSDPosition->Write() ;
1565 fhPhotReg->Write() ;
1569 fhNBarReg->Write() ;
1571 fhChargedReg->Write() ;
1581 fhChargedEM->Write() ;
1583 fhAllPPSD->Write() ;
1585 fhPhotPPSD->Write() ;
1589 fhNBarPPSD->Write() ;
1591 fhChargedPPSD->Write() ;
1593 fhPrimary->Write() ;
1603 fhPhotPhot->Write() ;
1605 fhPhotElec->Write() ;
1607 fhPhotNeuH->Write() ;
1609 fhPhotNuEM->Write() ;
1611 fhPhotNuEM->Write() ;
1613 fhPhotChHa->Write() ;
1615 fhPhotGaHa->Write() ;
1616 if(fhEnergyCorrelations)
1617 fhEnergyCorrelations->Write() ;
1622 //____________________________________________________________________________
1623 Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart)
1625 return ERecPart/0.8783 ;
1628 //____________________________________________________________________________
1629 void AliPHOSAnalyze::ResetHistograms()
1631 fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2)
1633 fhEmcDigit = 0 ; // Histo of digit energies in the Emc
1634 fhVetoDigit = 0 ; // Histo of digit energies in the Veto
1635 fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor
1636 fhEmcCluster = 0 ; // Histo of Cluster energies in Emc
1637 fhVetoCluster = 0 ; // Histo of Cluster energies in Veto
1638 fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor
1639 fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies
1642 fhPhotEnergy = 0 ; // Total spectrum of detected photons
1643 fhEMEnergy = 0 ; // Spectrum of detected electrons with electron primary
1646 fhAllPositionX = 0 ;
1647 fhAllPositionZ = 0 ;
1648 fhPhotPosition = 0 ;
1650 fhPPSDPosition = 0 ;