X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSAnalyze.cxx;h=0a0ec831e11e709f4802e67542a2ee81073fee07;hb=adcca1e6768f2a684cb2908f220b2cddc77a7caa;hp=cc3d6c6baf9a455f5954295c2c6de00ab111311b;hpb=055a87bcb268af8cedd043c9558c25a7ecf125f0;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSAnalyze.cxx b/PHOS/AliPHOSAnalyze.cxx index cc3d6c6baf9..0a0ec831e11 100644 --- a/PHOS/AliPHOSAnalyze.cxx +++ b/PHOS/AliPHOSAnalyze.cxx @@ -28,8 +28,7 @@ // Contamination(...) - calculates contamination of the photon spectrum and // pobability of reconstruction of several primaries as // kGAMMA,kELECTRON etc. -// -// User Case: +//// User Case: // root [0] AliPHOSAnalyze * a = new AliPHOSAnalyze("galice.root") // // set the file you want to analyse // root [1] a->DrawRecon(1,3) @@ -70,27 +69,24 @@ #include "TH2.h" #include "TParticle.h" #include "TClonesArray.h" -#include "TTree.h" #include "TMath.h" #include "TROOT.h" -#include "TFolder.h" // --- Standard library --- -#include - // --- AliRoot header files --- -#include "AliRun.h" -#include "AliPHOSv1.h" +//#include "AliRun.h" +#include "AliStack.h" +#include "AliPHOSGeometry.h" #include "AliPHOSAnalyze.h" #include "AliPHOSDigit.h" #include "AliPHOSSDigitizer.h" +#include "AliPHOSEmcRecPoint.h" +#include "AliPHOSCpvRecPoint.h" #include "AliPHOSTrackSegment.h" #include "AliPHOSRecParticle.h" -#include "AliPHOSCpvRecPoint.h" -#include "AliPHOSPpsdRecPoint.h" -#include "AliPHOSGetter.h" +#include "AliPHOSLoader.h" ClassImp(AliPHOSAnalyze) @@ -100,18 +96,25 @@ ClassImp(AliPHOSAnalyze) { // default ctor (useless) fCorrection = 1.2 ; //Value calculated for default parameters of reconstruction + fRunLoader = 0x0; } //____________________________________________________________________________ AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName) { // ctor: analyze events from root file "name" - ffileName = fileName ; + ffileName = fileName; fCorrection = 1.05 ; //Value calculated for default parameters of reconstruction + fRunLoader = AliRunLoader::Open(fileName,"AliPHOSAnalyze"); + if (fRunLoader == 0x0) + { + Error("AliPHOSAnalyze","Error Loading session"); + } } //____________________________________________________________________________ AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana) + : TObject(ana) { // copy ctor ( (AliPHOSAnalyze &)ana ).Copy(*this) ; @@ -124,43 +127,96 @@ AliPHOSAnalyze::~AliPHOSAnalyze() } //____________________________________________________________________________ -void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,const char* branchTitle){ +void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){ //Draws pimary particles and reconstructed //digits, RecPoints, RecPartices etc //for event Nevent in the module Nmod. + + //========== Create ObjectLoader + if (fRunLoader == 0x0) + { + Error("DrawRecon","Error Loading session"); + return; + } - TH2F * digitOccupancy = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.); - TH2F * sdigitOccupancy = new TH2F("sdigitOccupancy","EMC sdigits", 64,-71.,71.,64,-71.,71.); - TH2F * emcOccupancy = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.); - TH2F * ppsdUp = new TH2F("ppsdUp","PPSD Up digits", 128,-71.,71.,128,-71.,71.) ; - TH2F * ppsdUpCl = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ; - TH2F * ppsdLow = new TH2F("ppsdLow","PPSD Low digits", 128,-71.,71.,128,-71.,71.) ; - TH2F * ppsdLowCl = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ; - TH2F * nbar = new TH2F("nbar","Primary nbar", 64,-71.,71.,64,-71.,71.); - TH2F * phot = new TH2F("phot","Primary Photon", 64,-71.,71.,64,-71.,71.); - TH2F * charg = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.); - TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.); - TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", 64,-71.,71.,64,-71.,71.); + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + Error("DrawRecon","Could not obtain the Loader object !"); + return ; + } - //========== Create ObjectGetter - AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ; + + if(Nevent >= fRunLoader->GetNumberOfEvents() ) { + Error("DrawRecon", "There is no event %d only %d events available", Nevent, fRunLoader->GetNumberOfEvents() ) ; + return ; + } const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; - gime->Event(Nevent); + fRunLoader->GetEvent(Nevent); + + Int_t nx = phosgeom->GetNPhi() ; + Int_t nz = phosgeom->GetNZ() ; + Float_t * cri= phosgeom->GetEMCAGeometry()->GetCrystalHalfSize() ; + Float_t x = nx*cri[0] ; + Float_t z = nz*cri[2] ; + Int_t nxCPV = (Int_t) (nx*phosgeom->GetPadSizePhi()/(2.*cri[0])) ; + Int_t nzCPV = (Int_t) (nz*phosgeom->GetPadSizeZ()/(2.*cri[2])) ; + + TH2F * emcDigits = (TH2F*) gROOT->FindObject("emcDigits") ; + if(emcDigits) + emcDigits->Delete() ; + emcDigits = new TH2F("emcDigits","EMC digits", nx,-x,x,nz,-z,z); + TH2F * emcSdigits =(TH2F*) gROOT->FindObject("emcSdigits") ; + if(emcSdigits) + emcSdigits->Delete() ; + emcSdigits = new TH2F("emcSdigits","EMC sdigits", nx,-x,x,nz,-z,z); + TH2F * emcRecPoints = (TH2F*)gROOT->FindObject("emcRecPoints") ; + if(emcRecPoints) + emcRecPoints->Delete() ; + emcRecPoints = new TH2F("emcRecPoints","EMC RecPoints",nx,-x,x,nz,-z,z); + TH2F * cpvSdigits =(TH2F*) gROOT->FindObject("cpvSdigits") ; + if(cpvSdigits) + cpvSdigits->Delete() ; + cpvSdigits = new TH2F("cpvSdigits","CPV sdigits", nx,-x,x,nz,-z,z); + TH2F * cpvDigits = (TH2F*)gROOT->FindObject("cpvDigits") ; + if(cpvDigits) + cpvDigits->Delete() ; + cpvDigits = new TH2F("cpvDigits","CPV digits", nxCPV,-x,x,nzCPV,-z,z) ; + TH2F * cpvRecPoints= (TH2F*)gROOT->FindObject("cpvRecPoints") ; + if(cpvRecPoints) + cpvRecPoints->Delete() ; + cpvRecPoints = new TH2F("cpvRecPoints","CPV RecPoints", nxCPV,-x,x,nzCPV,-z,z) ; + + TH2F * phot = (TH2F*)gROOT->FindObject("phot") ; + if(phot) + phot->Delete() ; + phot = new TH2F("phot","Primary Photon", nx,-x,x,nz,-z,z); + TH2F * recPhot = (TH2F*)gROOT->FindObject("recPhot") ; + if(recPhot) + recPhot->Delete() ; + recPhot = new TH2F("recPhot","RecParticles with primary Photon",nx,-x,x,nz,-z,z); + //Plot Primary Particles + + if (fRunLoader->Stack() == 0x0) fRunLoader->LoadKinematics("READ"); + + const TParticle * primary ; Int_t iPrimary ; - for ( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++) + for ( iPrimary = 0 ; iPrimary < fRunLoader->Stack()->GetNprimary() ; iPrimary++) { - primary = gime->Primary(iPrimary) ; - Int_t primaryType = primary->GetPdgCode() ; - if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) { - Int_t moduleNumber ; - Double_t primX, primZ ; - phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; - if(moduleNumber==Nmod) - charg->Fill(primZ,primX,primary->Energy()) ; - } + primary = fRunLoader->Stack()->Particle(iPrimary); + + Int_t primaryType = primary->GetPdgCode(); +// if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) +// ||(primaryType == 11)||(primaryType == -11) ) { +// Int_t moduleNumber ; +// Double_t primX, primZ ; +// phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; +// if(moduleNumber==Nmod) +// charg->Fill(primZ,primX,primary->Energy()) ; +// } if( primaryType == 22 ) { Int_t moduleNumber ; Double_t primX, primZ ; @@ -168,213 +224,271 @@ void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,c if(moduleNumber==Nmod) phot->Fill(primZ,primX,primary->Energy()) ; } - else{ - if( primaryType == -2112 ) { - Int_t moduleNumber ; - Double_t primX, primZ ; - phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; - if(moduleNumber==Nmod) - nbar->Fill(primZ,primX,primary->Energy()) ; - } - } +// else{ +// if( primaryType == -2112 ) { +// Int_t moduleNumber ; +// Double_t primX, primZ ; +// phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; +// if(moduleNumber==Nmod) +// nbar->Fill(primZ,primX,primary->Energy()) ; +// } +// } } + Int_t iSDigit ; - const AliPHOSDigit * sdigit ; - - for(iSDigit = 0; iSDigit < gime->NSDigits(); iSDigit++) + AliPHOSDigit * sdigit ; + const TClonesArray * sdigits = gime->SDigits() ; + Int_t nsdig[5] = {0,0,0,0,0} ; + if(sdigits){ + for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++) { - sdigit = gime->SDigit(iSDigit) ; - Int_t relid[4]; - phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ; - Float_t x,z ; - phosgeom->RelPosInModule(relid,x,z) ; - Float_t e = gime->SDigitizer()->Calibrate(sdigit->GetAmp()) ; - if(relid[0]==Nmod){ - if(relid[1]==0) //EMC - sdigitOccupancy->Fill(x,z,e) ; - if((relid[1]>0)&&(relid[1]<17)) - ppsdUp->Fill(x,z,e) ; - if(relid[1]>16) - ppsdLow->Fill(x,z,e) ; - } + sdigit = (AliPHOSDigit *) sdigits->At(iSDigit) ; + Int_t relid[4]; + phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ; + Float_t x,z ; + phosgeom->RelPosInModule(relid,x,z); + AliPHOSSDigitizer* sd = dynamic_cast(gime->SDigitizer()); + Float_t e = sd->Calibrate(sdigit->GetAmp()) ; + nsdig[relid[0]-1]++ ; + if(relid[0]==Nmod){ + if(relid[1]==0) //EMC + emcSdigits->Fill(x,z,e) ; + if( relid[1]!=0 ) + cpvSdigits->Fill(x,z,e) ; + } } - + } + TString message ; + message = "Number of EMC + CPV SDigits per module: \n" ; + message += "%d %d %d %d %d\n"; + Info("DrawRecon", message.Data(), nsdig[0], nsdig[1], nsdig[2], nsdig[3], nsdig[4] ) ; + //Plot digits Int_t iDigit ; - const AliPHOSDigit * digit ; - for(iDigit = 0; iDigit < gime->NDigits(); iDigit++) - { - digit = gime->Digit(iDigit) ; - Int_t relid[4]; - phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; - Float_t x,z ; - phosgeom->RelPosInModule(relid,x,z) ; - Float_t e = gime->SDigitizer()->Calibrate(digit->GetAmp()) ; - if(relid[0]==Nmod){ - if(relid[1]==0) //EMC - digitOccupancy->Fill(x,z,e) ; - if((relid[1]>0)&&(relid[1]<17)) - ppsdUp->Fill(x,z,e) ; - if(relid[1]>16) - ppsdLow->Fill(x,z,e) ; + AliPHOSDigit * digit ; + const TClonesArray * digits = gime->Digits(); + if(digits) { + for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++) + { + digit = (AliPHOSDigit *) digits->At(iDigit) ; + Int_t relid[4]; + phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; + Float_t x,z ; + phosgeom->RelPosInModule(relid,x,z) ; + AliPHOSSDigitizer* sd = dynamic_cast(gime->SDigitizer()); + Float_t e = sd->Calibrate(digit->GetAmp()) ; + if(relid[0]==Nmod){ + if(relid[1]==0) //EMC + emcDigits->Fill(x,z,e) ; + if( relid[1]!=0 ) + cpvDigits->Fill(x,z,e) ; + } } - } - - + } + + //Plot RecPoints Int_t irecp ; TVector3 pos ; - - for(irecp = 0; irecp < gime->NEmcRecPoints() ; irecp ++){ - const AliPHOSEmcRecPoint * emc= gime->EmcRecPoint(irecp) ; - if(emc->GetPHOSMod()==Nmod){ - emc->GetLocalPosition(pos) ; - emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy()); - } - } - - - for(irecp = 0; irecp < gime->NCpvRecPoints() ; irecp ++){ - const AliPHOSRecPoint * cpv = gime->CpvRecPoint(irecp) ; - if((strcmp(cpv->ClassName(),"AliPHOSPpsdRecPoint" )) == 0){ // PPSD Rec Point - AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) cpv ; - ppsd->GetLocalPosition(pos) ; - if(ppsd->GetPHOSMod()==Nmod){ - ppsd->GetLocalPosition(pos) ; - if(ppsd->GetUp()) - ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy()); - else - ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy()); + TObjArray * emcrp = gime->EmcRecPoints() ; + if(emcrp) { + for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){ + AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ; + if(emc->GetPHOSMod()==Nmod){ + emc->GetLocalPosition(pos) ; + emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy()); } } - else{ - AliPHOSCpvRecPoint * cpv1 = (AliPHOSCpvRecPoint*) cpv ; - if(cpv1->GetPHOSMod()==Nmod){ - cpv1->GetLocalPosition(pos) ; - ppsdUpCl->Fill(pos.X(),pos.Z(),cpv1->GetEnergy()); + } + + TObjArray * cpvrp = gime->CpvRecPoints() ; + if(cpvrp) { + for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){ + AliPHOSRecPoint * cpv = (AliPHOSCpvRecPoint *) cpvrp->At(irecp) ; + if(cpv->GetPHOSMod()==Nmod){ + cpv->GetLocalPosition(pos) ; + cpvRecPoints->Fill(pos.X(),pos.Z(),cpv->GetEnergy()); } } } - - + //Plot RecParticles - const AliPHOSRecParticle * recParticle ; + AliPHOSRecParticle * recParticle ; Int_t iRecParticle ; - for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ ) - { - recParticle = gime->RecParticle(iRecParticle) ; - Int_t moduleNumberRec ; - Double_t recX, recZ ; - phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; - if(moduleNumberRec == Nmod){ - - Double_t minDistance = 5. ; - Int_t closestPrimary = -1 ; - - - //extract list of primaries: it is stored at EMC RecPoints - Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ; - Int_t numberofprimaries ; - Int_t * listofprimaries = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries) ; - Int_t index ; - const TParticle * primary ; - Double_t distance = minDistance ; - - for ( index = 0 ; index < numberofprimaries ; index++){ - primary = gime->Primary(listofprimaries[index]) ; - Int_t moduleNumber ; - Double_t primX, primZ ; - phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; - if(moduleNumberRec == moduleNumber) - distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ; - if(minDistance > distance) - { - minDistance = distance ; - closestPrimary = listofprimaries[index] ; - } - } - - if(closestPrimary >=0 ){ - - Int_t primaryType = gime->Primary(closestPrimary)->GetPdgCode() ; - - if(primaryType==22) - recPhot->Fill(recZ,recX,recParticle->Energy()) ; - else - if(primaryType==-2112) - recNbar->Fill(recZ,recX,recParticle->Energy()) ; - } + TClonesArray * rp = gime->RecParticles() ; + TClonesArray * ts = gime->TrackSegments() ; + if(rp && ts && emcrp) { + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ ) + { + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; + Int_t moduleNumberRec ; + Double_t recX, recZ ; + phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; + if(moduleNumberRec == Nmod){ + + Double_t minDistance = 5. ; + Int_t closestPrimary = -1 ; + + //extract list of primaries: it is stored at EMC RecPoints + Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; + Int_t numberofprimaries ; + Int_t * listofprimaries = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; + Int_t index ; + const TParticle * primary ; + Double_t distance = minDistance ; + + for ( index = 0 ; index < numberofprimaries ; index++){ + primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ; + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumberRec == moduleNumber) + distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ; + if(minDistance > distance) + { + minDistance = distance ; + closestPrimary = listofprimaries[index] ; + } + } + + if(closestPrimary >=0 ){ + + Int_t primaryType = fRunLoader->Stack()->Particle(closestPrimary)->GetPdgCode() ; + + if(primaryType==22) + recPhot->Fill(recZ,recX,recParticle->Energy()) ; +// else +// if(primaryType==-2112) +// recNbar->Fill(recZ,recX,recParticle->Energy()) ; + } + } } - } + } //Plot made histograms - digitOccupancy->Draw("box") ; - sdigitOccupancy->SetLineColor(5) ; - sdigitOccupancy->Draw("box") ; - emcOccupancy->SetLineColor(2) ; - emcOccupancy->Draw("boxsame") ; - ppsdUp->SetLineColor(3) ; - ppsdUp->Draw("boxsame") ; - ppsdLow->SetLineColor(4) ; - ppsdLow->Draw("boxsame") ; - phot->SetLineColor(8) ; - phot->Draw("boxsame") ; - nbar->SetLineColor(6) ; - nbar->Draw("boxsame") ; + emcSdigits->Draw("box") ; + emcDigits->SetLineColor(5) ; + emcDigits->Draw("boxsame") ; + emcRecPoints->SetLineColor(2) ; + emcRecPoints->Draw("boxsame") ; + cpvSdigits->SetLineColor(1) ; + cpvSdigits->Draw("boxsame") ; } //____________________________________________________________________________ void AliPHOSAnalyze::Ls(){ //lists branches and titles of PHOS-related branches of TreeR, TreeD, TreeS - AliPHOSGetter::GetInstance(ffileName.Data()) ; + if (fRunLoader == 0x0) + { + Error("Ls","Error Loading session"); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + Error("Ls","Could not obtain the Loader object !"); + return ; + } + Int_t ibranch; TObjArray * branches; - branches = gAlice->TreeS()->GetListOfBranches() ; + if (gime->TreeS() == 0x0) + { + if (gime->LoadSDigits("READ")) + { + Error("Ls","Problems with loading summable digits"); + return; + } + } + branches = gime->TreeS()->GetListOfBranches() ; - cout << "TreeS: " << endl ; + TString message ; + message = "TreeS:\n" ; for(ibranch = 0;ibranch GetEntries();ibranch++){ TBranch * branch=(TBranch *) branches->At(ibranch) ; - if(strstr(branch->GetName(),"PHOS") ) - cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ; + if(strstr(branch->GetName(),"PHOS") ){ + message += " " ; + message += branch->GetName() ; + message += " " ; + message += branch->GetTitle() ; + message += "\n" ; + } } - cout << endl ; - - branches = gAlice->TreeD()->GetListOfBranches() ; + if (gime->TreeD() == 0x0) + { + if (gime->LoadDigits("READ")) + { + Error("Ls","Problems with loading digits"); + return; + } + } + + branches = gime->TreeD()->GetListOfBranches() ; - cout << "TreeD: " << endl ; + message += "TreeD:\n" ; for(ibranch = 0;ibranch GetEntries();ibranch++){ TBranch * branch=(TBranch *) branches->At(ibranch) ; - if(strstr(branch->GetName(),"PHOS") ) - cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ; + if(strstr(branch->GetName(),"PHOS") ) { + message += " "; + message += branch->GetName() ; + message += " " ; + message += branch->GetTitle() ; + message +="\n" ; + } } - cout << endl ; + if (gime->TreeR() == 0x0) + { + if (gime->LoadRecPoints("READ")) + { + Error("Ls","Problems with loading rec points"); + return; + } + } - branches = gAlice->TreeR()->GetListOfBranches() ; + branches = gime->TreeR()->GetListOfBranches() ; - cout << "TreeR: " << endl ; + message += "TreeR: \n" ; for(ibranch = 0;ibranch GetEntries();ibranch++){ TBranch * branch=(TBranch *) branches->At(ibranch) ; - if(strstr(branch->GetName(),"PHOS") ) - cout << " " << branch->GetName() << " " << branch->GetTitle() << endl ; + if(strstr(branch->GetName(),"PHOS") ) { + message += " " ; + message += branch->GetName() ; + message += " " ; + message += branch->GetTitle() ; + message += "\n" ; + } } - cout << endl ; - - + Info("LS", message.Data()) ; } //____________________________________________________________________________ - void AliPHOSAnalyze::InvariantMass(const char* branchTitle) + void AliPHOSAnalyze::InvariantMass() { // Calculates Real and Mixed invariant mass distributions - AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ; - + if (fRunLoader == 0x0) + { + Error("DrawRecon","Error Loading session"); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + Error("DrawRecon","Could not obtain the Loader object !"); + return ; + } + + gime->LoadRecParticles("READ"); + Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution - + + //opening file TFile * mfile = new TFile("invmass.root","update"); @@ -410,76 +524,90 @@ void AliPHOSAnalyze::Ls(){ //scan over all events Int_t event ; - Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; + + + if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader(); + + + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries(); // for(event = 0; event < gime->MaxEvent(); event++ ){ + + + for(event = 0; event < maxevent; event++ ){ - gime->Event(event); + fRunLoader->GetEvent(event); //will read only TreeR //copy EM RecParticles to the "total" list const AliPHOSRecParticle * recParticle ; Int_t iRecParticle ; - for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ ) + TClonesArray * rp = gime->RecParticles() ; + if(!rp){ + Error("InvariantMass", "Can't find RecParticles") ; + return ; + } + + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ) { - recParticle = gime->RecParticle(iRecParticle) ; - if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)|| - (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)) - new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ; + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; + if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)) + new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ; } Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle nRecParticles[mevent] = iRecPhot-1 ; - + //check, if it is time to calculate invariant mass? Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; if((mevent == 0) && (event +1 == maxevent)){ - - // if((mevent == 0) && (event +1 == gime->MaxEvent())){ + + // if((mevent == 0) && (event +1 == gime->MaxEvent())){ //calculate invariant mass: Int_t irp1,irp2 ; Int_t nCurEvent = 0 ; for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){ - AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ; - - for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){ - AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ; - - Double_t invMass ; - invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())- - (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())- - (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())- - (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ; - - if(invMass> 0) - invMass = TMath::Sqrt(invMass); - - Double_t pt ; - pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) + - (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) ); - - if(irp1 > nRecParticles[nCurEvent]) - nCurEvent++; - - if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event - hRealEM->Fill(invMass,pt); - if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&& - (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) ) - hRealPhot->Fill(invMass,pt); - } - else{ - hMixedEM->Fill(invMass,pt); - if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&& - (rp2->GetType() == AliPHOSFastRecParticle::kGAMMA) ) - hMixedPhot->Fill(invMass,pt); - } //real-mixed - + AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ; + + for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){ + AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ; + + Double_t invMass ; + invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())- + (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())- + (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())- + (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ; + + if(invMass> 0) + invMass = TMath::Sqrt(invMass); + + Double_t pt ; + pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) + + (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) ); + + if(irp1 > nRecParticles[nCurEvent]) + nCurEvent++; + + if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event + hRealEM->Fill(invMass,pt); + if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&& + (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) ) + hRealPhot->Fill(invMass,pt); + } + else{ + hMixedEM->Fill(invMass,pt); + if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&& + (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) ) + hMixedPhot->Fill(invMass,pt); + } //real-mixed + } //loop over second rp }//loop over first rp //Make some cleanings for(Int_t index = 0; index < nMixedEvents; index ++) - nRecParticles[index] = 0 ; + nRecParticles[index] = 0 ; iRecPhot = 0 ; allRecParticleList->Clear() ; @@ -503,7 +631,7 @@ void AliPHOSAnalyze::Ls(){ } //____________________________________________________________________________ - void AliPHOSAnalyze::EnergyResolution(const char * branchTitle) + void AliPHOSAnalyze::EnergyResolution() { //fills two dimentional histo: energy of primary vs. energy of reconstructed @@ -527,21 +655,53 @@ void AliPHOSAnalyze::Ls(){ hEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.); - AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ; - const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; + if (fRunLoader == 0x0) + { + Error("DrawRecon","Error Loading session"); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + Error("DrawRecon","Could not obtain the Loader object !"); + return ; + } + + + const AliPHOSGeometry * phosgeom = gime->PHOSGeometry(); Int_t ievent; - Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; - // for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){ + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries(); + + fRunLoader->LoadKinematics("READ"); + gime->LoadTracks("READ"); + for ( ievent=0; ievent < maxevent ; ievent++){ //read the current event - gime->Event(ievent) ; + fRunLoader->GetEvent(ievent) ; const AliPHOSRecParticle * recParticle ; Int_t iRecParticle ; - for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ ){ - recParticle = gime->RecParticle(iRecParticle) ; + TClonesArray * rp = gime->RecParticles() ; + if(!rp) { + Error("EnergyResolution", "Event %d, Can't find RecParticles ", ievent) ; + return ; + } + TClonesArray * ts = gime->TrackSegments() ; + if(!ts) { + Error("EnergyResolution", "Event %d, Can't find TrackSegments", ievent) ; + return ; + } + TObjArray * emcrp = gime->EmcRecPoints() ; + if(!emcrp){ + Error("EnergyResolution", "Event %d, Can't find EmcRecPoints") ; + return ; + } + + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){ + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; //find the closest primary Int_t moduleNumberRec ; @@ -552,10 +712,10 @@ void AliPHOSAnalyze::Ls(){ Int_t closestPrimary = -1 ; //extract list of primaries: it is stored at EMC RecPoints - Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ; + Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; Int_t numberofprimaries ; - Int_t * listofprimaries = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries) ; - + Int_t * listofprimaries = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; + Int_t index ; const TParticle * primary ; Double_t distance = minDistance ; @@ -563,36 +723,38 @@ void AliPHOSAnalyze::Ls(){ Double_t dXmin = 0.; Double_t dZmin = 0. ; for ( index = 0 ; index < numberofprimaries ; index++){ - primary = gime->Primary(listofprimaries[index]) ; - Int_t moduleNumber ; - Double_t primX, primZ ; - phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; - if(moduleNumberRec == moduleNumber) { - dX = recX - primX; - dZ = recZ - primZ; - distance = TMath::Sqrt(dX*dX + dZ*dZ) ; - if(minDistance > distance) { - minDistance = distance ; - dXmin = dX; - dZmin = dZ; - closestPrimary = listofprimaries[index] ; - } - } + + primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ; + + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumberRec == moduleNumber) { + dX = recX - primX; + dZ = recZ - primZ; + distance = TMath::Sqrt(dX*dX + dZ*dZ) ; + if(minDistance > distance) { + minDistance = distance ; + dXmin = dX; + dZmin = dZ; + closestPrimary = listofprimaries[index] ; + } + } } //if found primary, fill histograms if(closestPrimary >=0 ){ - const TParticle * primary = gime->Primary(closestPrimary) ; - if(primary->GetPdgCode() == 22){ - hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ; - if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){ - hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; - hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; - } - else - if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) - hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; - } + const TParticle * primary = fRunLoader->Stack()->Particle(closestPrimary) ; + if(primary->GetPdgCode() == 22){ + hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ; + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){ + hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; + hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; + } + else + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) + hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; + } } } } @@ -608,7 +770,7 @@ void AliPHOSAnalyze::Ls(){ } //____________________________________________________________________________ -void AliPHOSAnalyze::PositionResolution(const char * branchTitle) +void AliPHOSAnalyze::PositionResolution() { //fills two dimentional histo: energy vs. primary - reconstructed distance @@ -628,40 +790,68 @@ void AliPHOSAnalyze::PositionResolution(const char * branchTitle) hAllPosition = (TH2F*)pfile->Get("hAllPosition"); if(hAllPosition == 0) hAllPosition = new TH2F("hAllPosition", - "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.); + "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.); hPhotPosition= (TH2F*)pfile->Get("hPhotPosition"); if(hPhotPosition == 0) hPhotPosition = new TH2F("hPhotPosition", - "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.); + "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.); hEMPosition= (TH2F*)pfile->Get("hEMPosition") ; if(hEMPosition == 0) hEMPosition = new TH2F("hEMPosition", - "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.); - hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ; + "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.); + hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ; if(hAllPositionX == 0) hAllPositionX = new TH1F("hAllPositionX", - "Delta X of any RP with primary photon",100, -2., 2.); + "Delta X of any RP with primary photon",100, -2., 2.); hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ; if(hAllPositionZ == 0) hAllPositionZ = new TH1F("hAllPositionZ", - "Delta X of any RP with primary photon",100, -2., 2.); + "Delta X of any RP with primary photon",100, -2., 2.); + if (fRunLoader == 0x0) + { + Error("DrawRecon","Error Loading session"); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + Error("DrawRecon","Could not obtain the Loader object !"); + return ; + } + + if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader(); - AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ; const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; Int_t ievent; - Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; - // for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){ + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries() ; for ( ievent=0; ievent < maxevent ; ievent++){ - + //read the current event - gime->Event(ievent) ; + fRunLoader->GetEvent(ievent) ; + TClonesArray * rp = gime->RecParticles() ; + if(!rp) { + Error("PositionResolution", "Event %d, Can't find RecParticles", ievent) ; + return ; + } + TClonesArray * ts = gime->TrackSegments() ; + if(!ts) { + Error("PositionResolution", "Event %d, Can't find TrackSegments", ievent) ; + return ; + } + TObjArray * emcrp = gime->EmcRecPoints() ; + if(!emcrp){ + Error("PositionResolution", "Event %d, Can't find EmcRecPoints", ievent) ; + return ; + } + const AliPHOSRecParticle * recParticle ; Int_t iRecParticle ; - for(iRecParticle = 0; iRecParticle < gime->NRecParticles() ;iRecParticle++ ){ - recParticle = gime->RecParticle(iRecParticle) ; + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){ + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; //find the closest primary Int_t moduleNumberRec ; @@ -672,9 +862,9 @@ void AliPHOSAnalyze::PositionResolution(const char * branchTitle) Int_t closestPrimary = -1 ; //extract list of primaries: it is stored at EMC RecPoints - Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ; + Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; Int_t numberofprimaries ; - Int_t * listofprimaries = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries) ; + Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; Int_t index ; const TParticle * primary ; @@ -684,38 +874,38 @@ void AliPHOSAnalyze::PositionResolution(const char * branchTitle) Double_t dXmin = 0.; Double_t dZmin = 0. ; for ( index = 0 ; index < numberofprimaries ; index++){ - primary = gime->Primary(listofprimaries[index]) ; - Int_t moduleNumber ; - Double_t primX, primZ ; - phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; - if(moduleNumberRec == moduleNumber) { - dX = recX - primX; - dZ = recZ - primZ; - distance = TMath::Sqrt(dX*dX + dZ*dZ) ; - if(minDistance > distance) { - minDistance = distance ; - dXmin = dX; - dZmin = dZ; - closestPrimary = listofprimaries[index] ; - } - } + primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ; + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumberRec == moduleNumber) { + dX = recX - primX; + dZ = recZ - primZ; + distance = TMath::Sqrt(dX*dX + dZ*dZ) ; + if(minDistance > distance) { + minDistance = distance ; + dXmin = dX; + dZmin = dZ; + closestPrimary = listofprimaries[index] ; + } + } } //if found primary, fill histograms if(closestPrimary >=0 ){ - const TParticle * primary = gime->Primary(closestPrimary) ; - if(primary->GetPdgCode() == 22){ - hAllPosition->Fill(primary->Energy(), minDistance) ; - hAllPositionX->Fill(primary->Energy(), dX) ; - hAllPositionZ->Fill(primary->Energy(), dZ) ; - if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){ - hPhotPosition->Fill(primary->Energy(), minDistance ) ; - hEMPosition->Fill(primary->Energy(), minDistance ) ; - } - else - if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) - hEMPosition->Fill(primary->Energy(), minDistance ) ; - } + const TParticle * primary = fRunLoader->Stack()->Particle(closestPrimary) ; + if(primary->GetPdgCode() == 22){ + hAllPosition->Fill(primary->Energy(), minDistance) ; + hAllPositionX->Fill(primary->Energy(), dX) ; + hAllPositionZ->Fill(primary->Energy(), dZ) ; + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){ + hPhotPosition->Fill(primary->Energy(), minDistance ) ; + hEMPosition->Fill(primary->Energy(), minDistance ) ; + } + else + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) + hEMPosition->Fill(primary->Energy(), minDistance ) ; + } } } } @@ -734,7 +924,7 @@ void AliPHOSAnalyze::PositionResolution(const char * branchTitle) } //____________________________________________________________________________ -void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ +void AliPHOSAnalyze::Contamination(){ // fills spectra of primary photons and several kinds of // reconstructed particles, so that analyzing them one can // estimate conatmination, efficiency of registration etc. @@ -743,7 +933,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ TH1F * hPrimary = 0; //spectrum (P_t distribution) of primary photons TH1F * hAllRP = 0; //spectrum of all RecParticles in PHOS TH1F * hPhot = 0; //spectrum of kGAMMA RecParticles - TH1F * hPPSD = 0; //spectrum of all RecParticles with PPSD signal TH1F * hShape = 0; //spectrum of all EM RecParticles TH1F * hVeto = 0; //spectrum of all neutral RecParticles @@ -751,27 +940,22 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ //primary - photon TH1F * hPhotReg = 0; //Registeres as photon TH1F * hPhotEM = 0; //Registered as EM - TH1F * hPhotPPSD= 0; //Registered as RecParticle with PPSD signal //primary - n TH1F * hNReg = 0; //Registeres as photon TH1F * hNEM = 0; //Registered as EM - TH1F * hNPPSD= 0; //Registered as RecParticle with PPSD signal //primary - nBar TH1F * hNBarReg = 0; //Registeres as photon TH1F * hNBarEM = 0; //Registered as EM - TH1F * hNBarPPSD= 0; //Registered as RecParticle with PPSD signal //primary - charged hadron (pBar excluded) TH1F * hChargedReg = 0; //Registeres as photon TH1F * hChargedEM = 0; //Registered as EM - TH1F * hChargedPPSD= 0; //Registered as RecParticle with PPSD signal //primary - pBar TH1F * hPbarReg = 0; //Registeres as photon TH1F * hPbarEM = 0; //Registered as EM - TH1F * hPbarPPSD= 0; //Registered as RecParticle with PPSD signal //Reading histograms from the file @@ -787,9 +971,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ hPhot = (TH1F*)cfile->Get("hPhot") ; if(hPhot == 0) hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.); - hPPSD = (TH1F*)cfile->Get("hPPSD") ; - if(hPPSD == 0) - hPPSD = new TH1F("hPPSD", "All PPSD Recparticles", 100, 0., 5.); hShape = (TH1F*) cfile->Get("hShape") ; if(hShape == 0) hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.); @@ -805,9 +986,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ hPhotEM =(TH1F*)cfile->Get("hPhotEM"); if(hPhotEM== 0) hPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.); - hPhotPPSD= (TH1F*)cfile->Get("hPhotPPSD"); - if(hPhotPPSD== 0) - hPhotPPSD = new TH1F("hPhotPPSD","Photon registered as PPSD", 100, 0., 5.); //primary - n hNReg = (TH1F*)cfile->Get("hNReg"); @@ -816,9 +994,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ hNEM = (TH1F*)cfile->Get("hNEM"); if(hNEM== 0) hNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.); - hNPPSD=(TH1F*)cfile->Get("hNPPSD"); - if(hNPPSD== 0) - hNPPSD = new TH1F("hNPPSD","N registered as PPSD", 100, 0., 5.); //primary - nBar hNBarReg =(TH1F*)cfile->Get("hNBarReg"); @@ -827,9 +1002,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ hNBarEM =(TH1F*)cfile->Get("hNBarEM"); if(hNBarEM== 0) hNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.); - hNBarPPSD=(TH1F*)cfile->Get("hNBarPPSD"); - if(hNBarPPSD== 0) - hNBarPPSD = new TH1F("hNBarPPSD","NBar registered as PPSD", 100, 0., 5.); //primary - charged hadron (pBar excluded) hChargedReg = (TH1F*)cfile->Get("hChargedReg"); @@ -838,10 +1010,7 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ hChargedEM = (TH1F*)cfile->Get("hChargedEM"); if(hChargedEM== 0) hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.); - hChargedPPSD= (TH1F*)cfile->Get("hChargedPPSD"); - if(hChargedPPSD== 0) - hChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.); - + //primary - pBar hPbarReg = (TH1F*)cfile->Get("hPbarReg"); if(hPbarReg== 0) @@ -849,9 +1018,6 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ hPbarEM = (TH1F*)cfile->Get("hPbarEM"); if(hPbarEM== 0) hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.); - hPbarPPSD= (TH1F*)cfile->Get("hPbarPPSD"); - if(hPbarPPSD== 0) - hPbarPPSD= new TH1F("hPbarPPSD","Pbar as PPSD",100, 0., 5.); //Now make some initializations @@ -864,42 +1030,72 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ - AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),RecPointsTitle) ; + if (fRunLoader == 0x0) + { + Error("DrawRecon","Error Loading session"); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + Error("DrawRecon","Could not obtain the Loader object !"); + return ; + } + + if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader(); const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; Int_t ievent; - Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; - // for ( ievent=0; ievent < gime->MaxEvent() ; ievent++){ + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries() ; for ( ievent=0; ievent < maxevent ; ievent++){ - gime->Event(ievent) ; + fRunLoader->GetEvent(ievent) ; + + TClonesArray * rp = gime->RecParticles() ; + if(!rp) { + Error("Contamination", "Event %d, Can't find RecParticles", ievent) ; + return ; + } + TClonesArray * ts = gime->TrackSegments() ; + if(!ts) { + Error("Contamination", "Event %d, Can't find TrackSegments", ievent) ; + return ; + } + TObjArray * emcrp = gime->EmcRecPoints() ; + if(!emcrp){ + Error("Contamination", "Event %d, Can't find EmcRecPoints", ievent) ; + return ; + } + //=========== Make spectrum of the primary photons const TParticle * primary ; Int_t iPrimary ; - for( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++){ - primary = gime->Primary(iPrimary) ; + for( iPrimary = 0 ; iPrimary < fRunLoader->Stack()->GetNprimary() ; iPrimary++){ + primary = fRunLoader->Stack()->Particle(iPrimary) ; Int_t primaryType = primary->GetPdgCode() ; if( primaryType == 22 ) { - //check, if photons folls onto PHOS - Int_t moduleNumber ; - Double_t primX, primZ ; - phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; - if(moduleNumber) - hPrimary->Fill(primary->Energy()) ; - + //check, if photons folls onto PHOS + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumber) + hPrimary->Fill(primary->Energy()) ; + } } + //========== Now scan over RecParticles const AliPHOSRecParticle * recParticle ; Int_t iRecParticle ; - for(iRecParticle = 0; iRecParticle < gime->NRecParticles(); iRecParticle++ ){ - recParticle = gime->RecParticle(iRecParticle) ; + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){ + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; //fill histo spectrum of all RecParticles hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ; - //==========find the closest primary + //==========find the closest primary Int_t moduleNumberRec ; Double_t recX, recZ ; phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; @@ -908,9 +1104,9 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ Int_t closestPrimary = -1 ; //extract list of primaries: it is stored at EMC RecPoints - Int_t emcIndex = gime->TrackSegment(recParticle->GetPHOSTSIndex())->GetEmcIndex() ; + Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; Int_t numberofprimaries ; - Int_t * listofprimaries = gime->EmcRecPoint(emcIndex)->GetPrimaries(numberofprimaries) ; + Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; Int_t index ; const TParticle * primary ; Double_t distance = minDistance ; @@ -918,128 +1114,105 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ Double_t dXmin = 0.; Double_t dZmin = 0. ; for ( index = 0 ; index < numberofprimaries ; index++){ - primary = gime->Primary(listofprimaries[index]) ; - Int_t moduleNumber ; - Double_t primX, primZ ; - phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; - if(moduleNumberRec == moduleNumber) { - dX = recX - primX; - dZ = recZ - primZ; - distance = TMath::Sqrt(dX*dX + dZ*dZ) ; - if(minDistance > distance) { - minDistance = distance ; - dXmin = dX; - dZmin = dZ; - closestPrimary = listofprimaries[index] ; - } - } + primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ; + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumberRec == moduleNumber) { + dX = recX - primX; + dZ = recZ - primZ; + distance = TMath::Sqrt(dX*dX + dZ*dZ) ; + if(minDistance > distance) { + minDistance = distance ; + dXmin = dX; + dZmin = dZ; + closestPrimary = listofprimaries[index] ; + } + } } //===========define the "type" of closest primary if(closestPrimary >=0 ){ - Int_t primaryCode = -1; - const TParticle * primary = gime->Primary(closestPrimary) ; - Int_t primaryType = primary->GetPdgCode() ; - if(primaryType == 22) // photon ? - primaryCode = 0 ; - else - if(primaryType == 2112) // neutron - primaryCode = 1 ; - else - if(primaryType == -2112) // Anti neutron - primaryCode = 2 ; - else - if(primaryType == -2122) //Anti proton - primaryCode = 4 ; - else { - TParticle tempo(*primary) ; - if(tempo.GetPDG()->Charge()) - primaryCode = 3 ; - } - - //==========Now look at the type of RecParticle - Float_t energy = CorrectedEnergy(recParticle->Energy()) ; - if(recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA){ - hPhot->Fill(energy ) ; - switch(primaryCode){ - case 0: - hPhotReg->Fill(energy ) ; - break ; - case 1: - hNReg->Fill(energy ) ; - break ; - case 2: - hNBarReg->Fill(energy ) ; - break ; - case 3: - hChargedReg->Fill(energy ) ; - break ; - case 4: - hPbarReg->Fill(energy ) ; - break ; - default: - break ; - } - } - if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)|| - (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA)){ //with PPSD signal - hPPSD->Fill(energy ) ; - switch(primaryCode){ - case 0: - hPhotPPSD->Fill(energy ) ; - break ; - case 1: - hNPPSD->Fill(energy ) ; - break ; - case 2: - hNBarPPSD->Fill(energy ) ; - break ; - case 3: - hChargedPPSD->Fill(energy ) ; - break ; - case 4: - hPbarPPSD->Fill(energy ) ; - break ; - default: - break ; - } - } - if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)|| - (recParticle->GetType() == AliPHOSFastRecParticle::kELECTRON)|| - (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)|| - (recParticle->GetType() == AliPHOSFastRecParticle::kABSURDEM) ){ //with EM shower - hShape->Fill(energy ) ; - switch(primaryCode){ - case 0: - hPhotEM->Fill(energy ) ; - break ; - case 1: - hNEM->Fill(energy ) ; - break ; - case 2: - hNBarEM->Fill(energy ) ; - break ; - case 3: - hChargedEM->Fill(energy ) ; - break ; - case 4: - hPbarEM->Fill(energy ) ; - break ; - default: - break ; - } - } - - if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)|| - (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHA) || - (recParticle->GetType() == AliPHOSFastRecParticle::kGAMMAHA) || - (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM) ) //nuetral - hVeto->Fill(energy ) ; - - //fill number of primaries identified as ... - if(primaryCode >= 0) // Primary code defined - counter[recParticle->GetType()][primaryCode]++ ; - + Int_t primaryCode = -1; + const TParticle * primary = fRunLoader->Stack()->Particle(closestPrimary) ; + Int_t primaryType = primary->GetPdgCode() ; + if(primaryType == 22) // photon ? + primaryCode = 0 ; + else + if(primaryType == 2112) // neutron + primaryCode = 1 ; + else + if(primaryType == -2112) // Anti neutron + primaryCode = 2 ; + else + if(primaryType == -2122) //Anti proton + primaryCode = 4 ; + else { + TParticle tempo(*primary) ; + if(tempo.GetPDG()->Charge()) + primaryCode = 3 ; + } + + //==========Now look at the type of RecParticle + Float_t energy = CorrectedEnergy(recParticle->Energy()) ; + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){ + hPhot->Fill(energy ) ; + switch(primaryCode){ + case 0: + hPhotReg->Fill(energy ) ; + break ; + case 1: + hNReg->Fill(energy ) ; + break ; + case 2: + hNBarReg->Fill(energy ) ; + break ; + case 3: + hChargedReg->Fill(energy ) ; + break ; + case 4: + hPbarReg->Fill(energy ) ; + break ; + default: + break ; + } + } + if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower + hShape->Fill(energy ) ; + switch(primaryCode){ + case 0: + hPhotEM->Fill(energy ) ; + break ; + case 1: + hNEM->Fill(energy ) ; + break ; + case 2: + hNBarEM->Fill(energy ) ; + break ; + case 3: + hChargedEM->Fill(energy ) ; + break ; + case 4: + hPbarEM->Fill(energy ) ; + break ; + default: + break ; + } + } + + if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) || + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) || + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral + hVeto->Fill(energy ) ; + + //fill number of primaries identified as ... + if(primaryCode >= 0) // Primary code defined + counter[recParticle->GetType()][primaryCode]++ ; + } } // no closest primary found @@ -1051,24 +1224,18 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ hPrimary->Write(0,kOverwrite); hAllRP->Write(0,kOverwrite); hPhot->Write(0,kOverwrite); - hPPSD->Write(0,kOverwrite); hShape->Write(0,kOverwrite); hVeto->Write(0,kOverwrite); hPhotReg->Write(0,kOverwrite); hPhotEM->Write(0,kOverwrite); - hPhotPPSD->Write(0,kOverwrite); hNReg ->Write(0,kOverwrite); hNEM ->Write(0,kOverwrite); - hNPPSD->Write(0,kOverwrite); hNBarReg ->Write(0,kOverwrite); hNBarEM ->Write(0,kOverwrite); - hNBarPPSD->Write(0,kOverwrite); hChargedReg ->Write(0,kOverwrite); hChargedEM ->Write(0,kOverwrite); - hChargedPPSD->Write(0,kOverwrite); hPbarReg ->Write(0,kOverwrite); hPbarEM ->Write(0,kOverwrite); - hPbarPPSD->Write(0,kOverwrite); cfile->Write(0,kOverwrite); cfile->Close(); @@ -1077,46 +1244,47 @@ void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){ //print Final Table maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; - // cout << "Resolutions: Analyzed " << gime->MaxEvent() << " event(s)" << endl ; - cout << "Resolutions: Analyzed " << maxevent << " event(s)" << endl ; - cout << endl ; - - cout << " Primary: Photon Neutron Antineutron Charged hadron AntiProton" << endl ; - cout << "--------------------------------------------------------------------------------" << endl ; - cout << " kGAMMA: " - << setw(8) << counter[2][0] << setw(9) << counter[2][1] << setw(13) << counter[2][2] - << setw(15)<< counter[2][3] << setw(13) << counter[2][4] << endl ; - cout << " kGAMMAHA: " - << setw(8) << counter[3][0] << setw(9) << counter[3][1] << setw(13) << counter[3][2] - << setw(15)<< counter[3][3] << setw(13) << counter[3][4] << endl ; - cout << " kNEUTRALEM: " - << setw(8) << counter[0][0] << setw(9) << counter[0][1] << setw(13) << counter[0][2] - << setw(15)<< counter[0][3] << setw(13) << counter[0][4] << endl ; - cout << " kNEUTRALHA: " - << setw(8) << counter[1][0] << setw(9) << counter[1][1] << setw(13) << counter[1][2] - << setw(15)<< counter[1][3] << setw(13) << counter[1][4] << endl ; - cout << " kABSURDEM: " - << setw(8) << counter[4][0] << setw(9) << counter[4][1] << setw(13) << counter[4][2] - << setw(15)<< counter[4][3] << setw(13) << counter[4][4] << endl ; - cout << " kABSURDHA: " - << setw(8) << counter[5][0] << setw(9) << counter[5][1] << setw(13) << counter[5][2] - << setw(15)<< counter[5][3] << setw(13) << counter[5][4] << endl ; - cout << " kELECTRON: " - << setw(8) << counter[6][0] << setw(9) << counter[6][1] << setw(13) << counter[6][2] - << setw(15)<< counter[6][3] << setw(13) << counter[6][4] << endl ; - cout << " kCHARGEDHA: " - << setw(8) << counter[7][0] << setw(9) << counter[7][1] << setw(13) << counter[7][2] - << setw(15)<< counter[7][3] << setw(13) << counter[7][4] << endl ; - cout << "--------------------------------------------------------------------------------" << endl ; - + TString message ; + message = "Resolutions: Analyzed %d event(s)\n" ; + + message += " Primary: Photon Neutron Antineutron Charged hadron AntiProton\n" ; + message += "--------------------------------------------------------------------------------\n" ; + message += " kGAMMA: " ; + message += "%d %d %d %d %d\n" ; + message += " kGAMMAHA: " ; + message += "%d %d %d %d %d\n" ; + message += " kNEUTRALEM: " ; + message += "%d %d %d %d %d\n" ; + message += " kNEUTRALHA: " ; + message += "%d %d %d %d %d\n" ; + message += " kABSURDEM: "; + message += "%d %d %d %d %d\n" ; + message += " kABSURDHA: " ; + message += "%d %d %d %d %d\n" ; + message += " kELECTRON: " ; + message += "%d %d %d %d %d\n" ; + message += " kCHARGEDHA: " ; + message += "%d %d %d %d %d\n" ; + + message += "--------------------------------------------------------------------------------" ; + + Int_t totalInd = 0 ; for(i1 = 0; i1<8; i1++) for(i2 = 0; i2<5; i2++) totalInd+=counter[i1][i2] ; - cout << "Indentified particles: " << totalInd << endl ; + message += "Indentified particles: %d" ; -} - - + Info("Contamination", message.Data(), maxevent, + counter[2][0], counter[2][1], counter[2][2], counter[2][3], counter[2][4], + counter[3][0], counter[3][1], counter[3][2], counter[3][3], counter[3][4], + counter[0][0], counter[0][1], counter[0][2], counter[0][3], counter[0][4], + counter[1][0], counter[1][1], counter[1][2], counter[1][3], counter[1][4], + counter[4][0], counter[4][1], counter[4][2], counter[4][3], counter[4][4], + counter[5][0], counter[5][1], counter[5][2], counter[5][3], counter[5][4], + counter[6][0], counter[6][1], counter[6][2], counter[6][3], counter[6][4], + counter[7][0], counter[7][1], counter[7][2], counter[7][3], counter[7][4], + totalInd ) ; +}