/* History of cvs commits:
*
* $Log$
+ * Revision 1.107 2007/05/25 14:12:26 policheh
+ * Local to tracking CS transformation added for CPV rec. points
+ *
+ * Revision 1.106 2007/05/24 13:01:22 policheh
+ * Local to tracking CS transformation invoked for each EMC rec.point
+ *
+ * Revision 1.105 2007/05/02 13:41:22 kharlov
+ * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
+ *
+ * Revision 1.104 2007/04/27 16:55:53 kharlov
+ * Calibration stops if PHOS CDB objects do not exist
+ *
+ * Revision 1.103 2007/04/11 11:55:45 policheh
+ * SetDistancesToBadChannels() added.
+ *
+ * Revision 1.102 2007/03/28 19:18:15 kharlov
+ * RecPoints recalculation in TSM removed
+ *
+ * Revision 1.101 2007/03/06 06:51:27 kharlov
+ * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
+ *
+ * Revision 1.100 2007/01/10 11:07:26 kharlov
+ * Raw digits writing to file (B.Polichtchouk)
+ *
* Revision 1.99 2006/11/07 16:49:51 kharlov
* Corrections for next event switch in case of raw data (B.Polichtchouk)
*
// --- AliRoot header files ---
#include "AliLog.h"
#include "AliRunLoader.h"
+#include "AliGenerator.h"
#include "AliPHOSGetter.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSClusterizerv1.h"
if(fToUnfold)
MakeUnfolding() ;
-
+
WriteRecPoints();
if(strstr(option,"deb"))
AliPHOSGetter * gime = AliPHOSGetter::Instance();
TClonesArray * digits = gime->Digits();
-
gMinuit->mncler(); // Reset Minuit's list of paramters
gMinuit->SetPrintLevel(-1) ; // No Printout
{
// Set calibration parameters:
// if calibration database exists, they are read from database,
- // otherwise, they are taken from digitizer.
+ // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
//
// It is a user responsilibity to open CDB before reconstruction, for example:
// AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
- AliPHOSGetter * gime = AliPHOSGetter::Instance();
- // fCalibData = new AliPHOSCalibData(gAlice->GetRunNumber()); //original
-
fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
-
- if(!fCalibData)
- {
- if ( !gime->Digitizer() )
- gime->LoadDigitizer();
- AliPHOSDigitizer * dig = gime->Digitizer();
- fADCchanelEmc = dig->GetEMCchannel() ;
- fADCpedestalEmc = dig->GetEMCpedestal();
-
- fADCchanelCpv = dig->GetCPVchannel() ;
- fADCpedestalCpv = dig->GetCPVpedestal() ;
- }
+ if (fCalibData->GetCalibDataEmc() == 0)
+ AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
+ if (fCalibData->GetCalibDataCpv() == 0)
+ AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
+
}
//____________________________________________________________________________
TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
TClonesArray * digits = gime->Digits() ;
-
Int_t index ;
//Evaluate position, dispersion and other RecPoint properties..
Int_t nEmc = emcRecPoints->GetEntriesFast();
for(index = 0; index < nEmc; index++){
- AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
+ AliPHOSEmcRecPoint * rp =
+ dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
rp->Purify(fEmcMinE) ;
- if(rp->GetMultiplicity()>0.) //If this RP is not empty
- rp->EvalAll(fW0,digits) ;
- else{
+ if(rp->GetMultiplicity()==0){
emcRecPoints->RemoveAt(index) ;
delete rp ;
+ continue;
}
+
+ // No vertex is available now, calculate corrections in PID
+ rp->EvalAll(fW0,digits) ;
+ TVector3 fakeVtx(0.,0.,0.) ;
+ rp->EvalAll(fW0,fakeVtx,digits) ;
+ rp->EvalLocal2TrackingCSTransform();
}
emcRecPoints->Compress() ;
- emcRecPoints->Sort() ;
+// emcRecPoints->Sort() ; //Can not sort until position is calculated!
// emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
for(index = 0; index < emcRecPoints->GetEntries(); index++){
dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
}
+ //For each rec.point set the distance to the nearest bad crystal (BVP)
+ SetDistancesToBadChannels();
+
//Now the same for CPV
for(index = 0; index < cpvRecPoints->GetEntries(); index++){
AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
rp->EvalAll(fW0CPV,digits) ;
+ rp->EvalLocal2TrackingCSTransform();
}
- cpvRecPoints->Sort() ;
+// cpvRecPoints->Sort() ;
for(index = 0; index < cpvRecPoints->GetEntries(); index++)
dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
}
//____________________________________________________________________________
-Double_t AliPHOSClusterizerv1::ShowerShape(Double_t r)
+Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
{
// Shape of the shower (see PHOS TDR)
// If you change this function, change also the gradient evaluation in ChiSquare()
- Double_t r4 = r*r*r*r ;
- Double_t r295 = TMath::Power(r, 2.95) ;
+ //for the moment we neglect dependence on the incident angle.
+
+ Double_t r2 = x*x + z*z ;
+ Double_t r4 = r2*r2 ;
+ Double_t r295 = TMath::Power(r2, 2.95/2.) ;
Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
return shape ;
}
Int_t nDigits = iniEmc->GetMultiplicity() ;
Float_t * efit = new Float_t[nDigits] ;
- Float_t xDigit=0.,zDigit=0.,distance=0. ;
+ Float_t xDigit=0.,zDigit=0. ;
Float_t xpar=0.,zpar=0.,epar=0. ;
Int_t relid[4] ;
AliPHOSDigit * digit = 0 ;
Int_t * emcDigits = iniEmc->GetDigitsList() ;
+ TVector3 vIncid ;
+
Int_t iparam ;
Int_t iDigit ;
for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
zpar = fitparameters[iparam+1] ;
epar = fitparameters[iparam+2] ;
iparam += 3 ;
- distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
- distance = TMath::Sqrt(distance) ;
- efit[iDigit] += epar * ShowerShape(distance) ;
+// geom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
+// efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
+ efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
}
}
zpar = fitparameters[iparam+1] ;
epar = fitparameters[iparam+2] ;
iparam += 3 ;
+// geom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
AliPHOSEmcRecPoint * emcRP = 0 ;
digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ;
geom->AbsToRelNumbering(digit->GetId(), relid) ;
geom->RelPosInModule(relid, xDigit, zDigit) ;
- distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
- distance = TMath::Sqrt(distance) ;
- ratio = epar * ShowerShape(distance) / efit[iDigit] ;
+// ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
+ ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
eDigit = emcEnergies[iDigit] * ratio ;
emcRP->AddDigit( *digit, eDigit ) ;
}
AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
-
-
+// TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(2)) ; //Vertex position
// AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
Float_t * emcEnergies = emcRP->GetEnergiesList() ;
const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
+// TVector3 vInc ;
fret = 0. ;
Int_t iparam ;
Int_t iParam = 0 ;
efit = 0 ;
while(iParam < nPar ){
- Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
+ Double_t dx = (xDigit - x[iParam]) ;
iParam++ ;
- distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
- distance = TMath::Sqrt( distance ) ;
+ Double_t dz = (zDigit - x[iParam]) ;
iParam++ ;
- efit += x[iParam] * ShowerShape(distance) ;
+// geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
+// efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
+ efit += x[iParam] * ShowerShape(dx,dz) ;
iParam++ ;
}
Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
Double_t xpar = x[iParam] ;
Double_t zpar = x[iParam+1] ;
Double_t epar = x[iParam+2] ;
+// geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
- Double_t shape = sum * ShowerShape(dr) ;
+// Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
+ Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
+//DP: No incident angle dependence in gradient yet!!!!!!
Double_t r4 = dr*dr*dr*dr ;
Double_t r295 = TMath::Power(dr,2.95) ;
Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
Double_t zpar = x[iparam+1] ;
Double_t epar = x[iparam+2] ;
iparam += 3 ;
- Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
- distance = TMath::Sqrt(distance) ;
- efit += epar * ShowerShape(distance) ;
+// geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
+// efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
+ efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
}
fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
fCpvLocMaxCut,
fW0CPV )) ;
}
-
-
+//____________________________________________________________________________
+//void AliPHOSClusterizerv1::GetVertex(void)
+//{ //Extracts vertex posisition
+//
+ //ESD
+//DP - todo if(){
+//
+// }
+
+// //MC Generator
+// if(gAlice && gAlice->GetMCApp() && gAlice->Generator()){
+// Float_t x,y,z ;
+// gAlice->Generator()->GetOrigin(x,y,z) ;
+// fVtx.SetXYZ(x,y,z) ;
+// return ;
+// }
+//
+// //No any source
+// fVtx[0]=fVtx[1]=fVtx[2]=0. ;
+//
+//}
//____________________________________________________________________________
void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
{
}
}
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::SetDistancesToBadChannels()
+{
+ //For each EMC rec. point set the distance to the nearest bad crystal.
+ //Author: Boris Polichtchouk
+
+ if(!fCalibData->GetNumOfEmcBadChannels()) return;
+ AliInfo(Form("%d bad channel(s) found.\n",fCalibData->GetNumOfEmcBadChannels()));
+
+ AliPHOSGetter* gime = AliPHOSGetter::Instance();
+ AliPHOSGeometry* geom = gime->PHOSGeometry();
+
+ TObjArray * emcRecPoints = gime->EmcRecPoints() ;
+
+ Int_t badIds[8000];
+ fCalibData->EmcBadChannelIds(badIds);
+
+ AliPHOSEmcRecPoint* rp;
+
+ TMatrixF gmat;
+ TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
+ TVector3 gposBadChannel; // global position of bad crystal
+ TVector3 dR;
+
+ Float_t dist,minDist;
+
+ for(Int_t iRP=0; iRP<emcRecPoints->GetEntries(); iRP++){
+ rp = (AliPHOSEmcRecPoint*)emcRecPoints->At(iRP);
+ minDist = 1.e+07;
+
+ for(Int_t iBad=0; iBad<fCalibData->GetNumOfEmcBadChannels(); iBad++) {
+ rp->GetGlobalPosition(gposRecPoint,gmat);
+ geom->RelPosInAlice(badIds[iBad],gposBadChannel);
+ AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
+ gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
+ gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
+ dR = gposBadChannel-gposRecPoint;
+ dist = dR.Mag();
+ if(dist<minDist) minDist = dist;
+ }
+
+ rp->SetDistanceToBadCrystal(minDist);
+ }
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff