int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
if( index >= 0)
+ {
+ Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
+ Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
+ short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
+ // timebinOffset is timebin value at maximum (maxrev)
+ short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
+ if( maxf < fAmpCut || ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
{
- Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
- Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
- short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
- // timebinOffset is timebin value at maximum (maxrev)
- short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
- if( maxf < fAmpCut || ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
- {
- return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
- }
- else if ( maxf >= fAmpCut )
- {
- int first = 0;
- int last = 0;
- SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut );
- int nsamples = last - first + 1;
-
- if( ( nsamples ) >= fNsampleCut )
+ return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
+ }
+ else if ( maxf >= fAmpCut )
+ {
+ int first = 0;
+ int last = 0;
+ SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut );
+ int nsamples = last - first + 1;
+
+ if( ( nsamples ) >= fNsampleCut )
{
Float_t tmax = (maxrev - first); // local tmax estimate
TGraph *graph = new TGraph( nsamples, fXaxis, &fReversed[first] );
fTf1->SetParameter(0, maxf*fkEulerSquared );
- fTf1->SetParameter(1, tmax - fTau);
+ fTf1->SetParameter(1, tmax - fTau);
// set rather loose parameter limits
fTf1->SetParLimits(0, 0.5*maxf*fkEulerSquared, 2*maxf*fkEulerSquared );
- fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4);
-
+ fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4);
+
if (fFixTau) {
- fTf1->FixParameter(2, fTau);
+ fTf1->FixParameter(2, fTau);
}
else {
- fTf1->ReleaseParameter(2); // allow par. to vary
- fTf1->SetParameter(2, fTau);
+ fTf1->ReleaseParameter(2); // allow par. to vary
+ fTf1->SetParameter(2, fTau);
}
-
+
Short_t tmpStatus = 0;
try {
- tmpStatus = graph->Fit(fTf1, "Q0RW");
+ tmpStatus = graph->Fit(fTf1, "Q0RW");
}
catch (const std::exception & e) {
- AliError( Form("TGraph Fit exception %s", e.what()) );
- return AliCaloFitResults( maxamp, ped, Ret::kNoFit, maxf, timebinOffset,
- timebinOffset, Ret::kDummy, Ret::kDummy, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
+ AliError( Form("TGraph Fit exception %s, fit status %d", e.what(),tmpStatus) );
+ return AliCaloFitResults( maxamp, ped, Ret::kNoFit, maxf, timebinOffset,
+ timebinOffset, Ret::kDummy, Ret::kDummy, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
}
-
+
if( fVerbose == true )
- {
- AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
- PrintFitResult( fTf1 ) ;
- }
+ {
+ AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
+ PrintFitResult( fTf1 ) ;
+ }
// global tmax
tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0
- + fTf1->GetParameter(2); // +tau, makes sum tmax
+ + fTf1->GetParameter(2); // +tau, makes sum tmax
- delete graph;
- return AliCaloFitResults( maxamp, ped , Ret::kFitPar,
- fTf1->GetParameter(0)/fkEulerSquared,
- tmax,
- timebinOffset,
- fTf1->GetChisquare(),
- fTf1->GetNDF(),
- Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
+ delete graph;
+ return AliCaloFitResults( maxamp, ped , Ret::kFitPar,
+ fTf1->GetParameter(0)/fkEulerSquared,
+ tmax,
+ timebinOffset,
+ fTf1->GetChisquare(),
+ fTf1->GetNDF(),
+ Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
- // delete graph;
-
+ // delete graph;
+
}
- else
+ else
{
Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
Int_t ndf = last - first - 1; // nsamples - 2
return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset,
- timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
+ timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
}
- } // ampcut
- }
+ } // ampcut
+ }
return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
}
if(l == 0) return;
if(ind>=0 && ind < l->GetSize()){
h = dynamic_cast<TProfile *>(l->At(ind));
- if(h>0) h->Fill(x,y,w);
+ if(h) h->Fill(x,y,w);
}
}
// Note: should always return 0=OK, because otherwise all tracking
// is aborted for this event
- if (!esd) {
+ if (!esd)
+ {
AliError("NULL ESD passed");
return 1;
}
// step 1: collect clusters
Int_t okLoadClusters, nClusters;
- if (!fClusters || (fClusters && fClusters->IsEmpty())) {
+
+ if (!fClusters || (fClusters && fClusters->IsEmpty()))
okLoadClusters = LoadClusters(esd);
- }
+
nClusters = fClusters->GetEntries();
// step 2: collect ESD tracks
okLoadTracks = LoadTracks(esd);
nTracks = fTracks->GetEntries();
+ AliDebug(5,Form("Propagate back %d tracks ok %d, for %d clusters ok %d",
+ nTracks,okLoadTracks,nClusters,okLoadClusters));
+
// step 3: for each track, find the closest cluster as matched within residual cuts
Int_t index=-1;
- for (Int_t it = 0; it < nTracks; it++) {
+ for (Int_t it = 0; it < nTracks; it++)
+ {
AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
index = FindMatchedCluster(track);
- if (index>-1) {
+ if (index>-1)
+ {
AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
track->SetEMCALcluster(cluster->Index());
track->SetStatus(AliESDtrack::kEMCALmatch);
// Otherwise use the TPCInner point
AliExternalTrackParam *trkParam = 0;
- if (!fITSTrackSA) {
+ if (!fITSTrackSA)
+ {
const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
+
if (friendTrack && friendTrack->GetTPCOut())
trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
else if (track->GetInnerParam())
AliExternalTrackParam trkParamTmp(*trkParam);
Float_t eta, phi, pt;
- if (!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, fEMCalSurfaceDistance, track->GetMass(kTRUE), fStep, eta, phi, pt)) {
+ if (!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, fEMCalSurfaceDistance, track->GetMass(kTRUE), fStep, eta, phi, pt))
+ {
if (fITSTrackSA) delete trkParam;
return index;
}
+
track->SetTrackPhiEtaPtOnEMCal(phi,eta,pt);
- if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) {
+
+ if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
+ {
if (fITSTrackSA) delete trkParam;
return index;
}
Double_t trkPos[3];
trkParamTmp.GetXYZ(trkPos);
Int_t nclusters = fClusters->GetEntries();
- for (Int_t ic=0; ic<nclusters; ic++) {
+ for (Int_t ic=0; ic<nclusters; ic++)
+ {
AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
- Float_t clsPos[3] = {static_cast<Float_t>(cluster->X()),static_cast<Float_t>(cluster->Y()),static_cast<Float_t>(cluster->Z())};
+
+ Float_t clsPos[3] = {static_cast<Float_t>(cluster->X()),
+ static_cast<Float_t>(cluster->Y()),
+ static_cast<Float_t>(cluster->Z())};
+
Double_t dR = TMath::Sqrt(TMath::Power(trkPos[0]-clsPos[0],2)+TMath::Power(trkPos[1]-clsPos[1],2)+TMath::Power(trkPos[2]-clsPos[2],2));
//printf("\n dR=%f,wind=%f\n",dR,fClusterWindow); //MARCEL
+
if (dR > fClusterWindow) continue;
AliExternalTrackParam trkParTmp(trkParamTmp);
Float_t tmpEta, tmpPhi;
if (!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
- if (TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta)) {
+
+ if (TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
+ {
maxPhi=tmpPhi;
maxEta=tmpEta;
index=ic;
}
if (fITSTrackSA) delete trkParam;
+
return index;
}
{
// Return sum of amplidutes from EMCal
// Used calibration coefficeint for transition to energy
- return fAmpJetMatrix >0 ?fAmpJetMatrix->Sum() :0.0;
+ return fAmpJetMatrix ?fAmpJetMatrix->Sum() :0.0;
}
//____________________________________________________________________________
gROOT->cd();
TH1::AddDirectory(1);
- TH1 *h = 0, *hgid=0;
+ TH1 * hgid=0;
Int_t nphi=180, nmax=1100;
Double_t phimin=0.0, phimax=360.;
Double_t pmax=110.;
if(p>0.1) pmax = 1.1*p;
- h = new TH1F("00_hNPrim"," number of primary particles ", 10, 0.5, 10.5);
+ new TH1F("00_hNPrim"," number of primary particles ", 10, 0.5, 10.5);
hgid = new TH1F("01_hGidprim","Geant Id of primary particles", 16, 0.5, 16.5);
new TH1F("02_hPmomPrim","p of primary particles", nmax, 0.0, pmax);
new TH1F("03_hEtaPrim","#eta primary particles", 80, 0.0, 8.0);
if(!gn.Contains("WSUC")) { // ALICE
AliMatrix(fIdRotm, 90.-angle,180., 90.0,90.0, angle, 0.);
phiOK = mod->GetCenterOfModule().Phi()*180./TMath::Pi();
- // printf(" %2i | angle | %6.3f - %6.3f = %6.3f(eta %5.3f)\n",
- //iz+1, angle, phiOK, angle-phiOK, mod->GetEtaOfCenterOfModule());
+ AliDebug(4,Form(" %2i | angle | %6.3f - %6.3f = %6.3f(eta %5.3f)\n",
+ iz+1, angle, phiOK, angle-phiOK, mod->GetEtaOfCenterOfModule()));
xpos = mod->GetPosXfromR() + g->GetSteelFrontThickness() - fSmodPar0;
zpos = mod->GetPosZ() - fSmodPar2;
if(iz == 0) AliMatrix(fIdRotm, 0.,0., 90.,0., 90.,90.); // (x')z; y'(x); z'(y)
else AliMatrix(fIdRotm, 90-angle,270., 90.0,0.0, angle,90.);
phiOK = mod->GetCenterOfModule().Phi()*180./TMath::Pi();
- //printf(" %2i | angle -phiOK | %6.3f - %6.3f = %6.3f(eta %5.3f)\n",
- //iz+1, angle, phiOK, angle-phiOK, mod->GetEtaOfCenterOfModule());
+ AliDebug(4,Form(" %2i | angle -phiOK | %6.3f - %6.3f = %6.3f(eta %5.3f)\n",
+ iz+1, angle, phiOK, angle-phiOK, mod->GetEtaOfCenterOfModule()));
zpos = mod->GetPosZ() - fSmodPar2;
ypos = mod->GetPosXfromR() - fSmodPar1;
//printf(" zpos %7.2f ypos %7.2f fIdRotm %i\n xpos ", zpos, xpos, fIdRotm);
{ // flat in phi
xpos = g->GetPhiModuleSize()*(2*ix+1 - g->GetNPhi())/2.;
TVirtualMC::GetMC()->Gspos(child, ++nr, mother, xpos, ypos, zpos, fIdRotm, "ONLY") ;
- printf(" %7.2f ", xpos);
+ //printf(" %7.2f ", xpos);
}
- printf("\n");
+ //printf("\n");
}
}
AliDebug(2,Form(" Number of modules in Super Module(%s) %i \n", mother, nr));
name.Data(), pbShape.Data(), pbtiChonly.Data(), nr);
AliEMCALGeometry * g = GetGeometry();
- double par[5], parPB[5], parSC[5];
+ double par[5], parPB[5];//, parSC[5];
double xpos = 0.0, ypos = 0.0;
double zpos = -fSampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
if(name == "SCMX") { // common trapezoid - 11 parameters
double tanx = (parTRAP[1] - parTRAP[0]) / (2.*parTRAP[4]); // tanx = tany now
double tany = (parTRAP[3] - parTRAP[2]) / (2.*parTRAP[4]), ztmp=0.;
parPB[4] = g->GetECPbRadThick()/2.;
- parSC[2] = g->GetECScintThick()/2.;
+ //parSC[2] = g->GetECScintThick()/2.;
for(int iz=0; iz<g->GetNECLayers(); iz++){
ztmp = fSampleWidth*double(iz);
parPB[0] = parTRAP[0] + tanx*ztmp;