#include <TVectorD.h>
#include <TSystem.h>
#include <TFile.h>
-
+#include <TList.h>
+#include <TTree.h>
// ALIROOT includes
#include <TTreeStream.h>
#include <AliAnalysisManager.h>
#include "AliMCEventHandler.h"
#include <AliESD.h>
+#include "AliTrackComparison.h"
#include "AliMCTrackingTestTask.h"
#include "AliGenInfoMaker.h"
#include "AliHelix.h"
+#include "AliTrackPointArray.h"
+#include "AliESDCaloCluster.h"
//
#include "AliMCInfo.h"
+#include "AliComparisonObject.h"
#include "AliESDRecInfo.h"
#include "AliTPCParamSR.h"
#include "AliTracker.h"
#include "AliTPCseed.h"
+#include "AliCDBManager.h"
+#include "AliGRPManager.h"
+#include "AliGeomManager.h"
+
// STL includes
#include <iostream>
ClassImp(AliMCTrackingTestTask)
+#define USE_STREAMER 0
+#define DEBUG 0
+
//________________________________________________________________________
AliMCTrackingTestTask::AliMCTrackingTestTask() :
AliAnalysisTask(),
fMCinfo(0), //! MC event handler
fESD(0),
+ fCurrentRun(-1),
fDebugStreamer(0),
fStreamLevel(0),
fDebugLevel(0),
- fDebugOutputPath()
+ fDebugOutputPath(),
+ fOutList(NULL),
+ fPitList(NULL),
+ fCompList(NULL)
{
//
// Default constructor (should not be used)
AliAnalysisTask(),
fMCinfo(0), //! MC event handler
fESD(0),
- //
+ fCurrentRun(-1),
fDebugStreamer(0),
fStreamLevel(0),
fDebugLevel(),
- fDebugOutputPath()
+ fDebugOutputPath(),
+ fOutList(NULL),
+ fPitList(NULL),
+ fCompList(NULL)
{
//
// Default constructor
AliAnalysisTask(name, "AliMCTrackingTestTask"),
fMCinfo(0), //! MC event handler
fESD(0),
+ fCurrentRun(-1),
fDebugStreamer(0),
fStreamLevel(0),
fDebugLevel(0),
- fDebugOutputPath()
+ fDebugOutputPath(),
+ fOutList(NULL),
+ fPitList(NULL),
+ fCompList(NULL)
{
//
// Normal constructor
//
// Input slot #0 works with a TChain
DefineInput(0, TChain::Class());
- // Output slot #0 writes into a TList
- DefineOutput(0, TObjArray::Class());
+ DefineOutput(0, TList::Class());
+
+ // create the list for comparison objects
+ fCompList = new TList;
//
//
}
if (fDebugLevel>0) printf("AliMCTrackingTestTask::~AliMCTrackingTestTask\n");
if (fDebugStreamer) delete fDebugStreamer;
fDebugStreamer=0;
+ if (fOutList) delete fOutList; fOutList = 0;
+ if (fCompList) delete fCompList; fCompList = 0;
}
+//_____________________________________________________________________________
+Bool_t AliMCTrackingTestTask::AddComparisonObject(AliTrackComparison *cObj)
+{
+ // add comparison object to the list
+ if(cObj == 0) {
+ Printf("ERROR: Could not add comparison object");
+ return kFALSE;
+ }
+
+ // add object to the list
+ fCompList->AddLast(cObj);
+
+ return kTRUE;
+}
//________________________________________________________________________
void AliMCTrackingTestTask::ConnectInputData(Option_t *)
{
mcinfo->SetReadTR(kTRUE);
fMCinfo = mcinfo->MCEvent();
-
-
}
-
-
-
-
-
//________________________________________________________________________
void AliMCTrackingTestTask::CreateOutputObjects()
{
if(fDebugLevel>3)
cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
+ if (!fOutList)
+ fOutList = new TList;
+ fOutList->SetOwner(kTRUE);
+ fPitList = fOutList->MakeIterator();
+
+ // create output list
+
+ // add comparison objects to the output
+ AliTrackComparison *cObj=0;
+ Int_t count=0;
+ TIterator *pitCompList = fCompList->MakeIterator();
+ pitCompList->Reset();
+ while(( cObj = (AliTrackComparison*)pitCompList->Next()) != NULL) {
+ fOutList->Add(cObj);
+ count++;
+ }
+ Printf("UserCreateOutputObjects(): Number of output comparison objects: %d \n", count);
+
+ PostData(0, fOutList);
}
// Execute analysis for current event
//
+#if DEBUG
+ printf("New event!\n");
+#endif
+
if(fDebugLevel>3)
cout << "AliMCTrackingTestTask::Exec()" << endl;
- // If MC has been connected
+ // If MC has been connected
if (!fMCinfo){
cout << "Not MC info\n" << endl;
}else{
ProcessMCInfo();
- //mcinfo->Print();
+ // fMCinfo->Print();
//DumpInfo();
}
- //
-}
-
+ PostData(0, fOutList);
+}
//________________________________________________________________________
if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
if (fDebugStreamer) delete fDebugStreamer;
fDebugStreamer = 0;
+
+// AliTrackComparison *cObj=0;
+// TIterator *pitCompList = fCompList->MakeIterator();
+// pitCompList->Reset();
+// while(( cObj = (AliTrackComparison*)pitCompList->Next()) != NULL) {
+// for(Int_t i=0; i<6; i++)
+// cObj->MakeDistortionMap(438,i);
+// }
+
return;
}
-
+//________________________________________________________________________
TTreeSRedirector *AliMCTrackingTestTask::GetDebugStreamer(){
//
// Get Debug streamer
TString dsName;
dsName=GetName();
dsName+="Debug.root";
+
+ printf(" get debug streamer \n");
+
dsName.ReplaceAll(" ","");
fDebugStreamer = new TTreeSRedirector(dsName.Data());
return fDebugStreamer;
}
-
-
-
-
-
+//________________________________________________________________________
AliExternalTrackParam * AliMCTrackingTestTask::MakeTrack(const AliTrackReference* ref, TParticle*part)
{
//
return param;
}
+//________________________________________________________________________
Bool_t AliMCTrackingTestTask::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
//
// Propagate track to point xyz using
- // AliTracker::PropagateTo functionality
+ // AliTracker::PropagateToBxByBz functionality
//
// param - track parameters
// xyz - position to propagate
// mass - particle mass
// step - step to be used
Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
- AliTracker::PropagateTrackTo(param, radius+step, mass, step, kTRUE,0.99);
- AliTracker::PropagateTrackTo(param, radius+0.5, mass, step*0.1, kTRUE,0.99);
- Double_t sxyz[3]={0,0,0};
- AliESDVertex vertex(xyz,sxyz);
- Bool_t isOK = param->PropagateToDCA(&vertex,AliTracker::GetBz(),10);
+ Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
+
+ AliTracker::PropagateTrackToBxByBz(param, radius+step, mass, step, kTRUE,0.99,-1);
+ param->Rotate(alpha);
+ Bool_t isOK = param->PropagateTo(radius,AliTracker::GetBz());
+
return isOK;
}
-
+//________________________________________________________________________
void AliMCTrackingTestTask::ProcessMCInfo(){
//
//
//
//
- TParticle * particle= new TParticle;
+#if DEBUG
+ printf("ProcessMCInfo\n");
+#endif
+
+ const AliESDVertex *pVertex = fESD->GetPrimaryVertex();
+ if(!pVertex) AliError("No primary vertex found!\n");
+ Double_t vPos[3];
+ pVertex->GetXYZ(vPos);
+
+ TParticle * particle= new TParticle();
TClonesArray * trefs = new TClonesArray("AliTrackReference");
const Double_t kPcut=0.1;
//
//
//
+
for (Int_t ipart=0;ipart<npart;ipart++){
Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
- if (status<0) continue;
- if (!particle) continue;
- if (!trefs) continue;
+ if (status<0 || !particle || !trefs) continue;
Int_t nref = trefs->GetEntries();
if (nref<5) continue;
FitTrackRefs(particle,trefs);
+
AliTrackReference * tpcIn=0;
AliTrackReference * tpcOut=0;
AliTrackReference * trdIn=0;
AliTrackReference * trdOut=0;
AliTrackReference * itsIn=0;
AliTrackReference * itsOut=0;
+
+ AliTrackReference * tofIn=0;
+ AliTrackReference * tofOut=0;
+ AliTrackReference * hmpidIn=0;
+ AliTrackReference * hmpidOut=0;
+ AliTrackReference * emcalIn=0;
+ AliTrackReference * emcalOut=0;
+
Double_t rmax=0;
Double_t rmin=1000;
for (Int_t iref=0;iref<nref;iref++){
if (ref->R()>trdIn->R()) trdOut = ref;
}
}
+ //TOF
+ if (ref->DetectorId()==AliTrackReference::kTOF){
+ if (!tofIn) {
+ tofIn = ref;
+ }else{
+ if (ref->R()>tofIn->R()) tofOut = ref;
+ }
+ }
+
+ //HMPID
+ if (ref->DetectorId()==AliTrackReference::kHMPID){
+ if (!hmpidIn) {
+ hmpidIn = ref;
+ }else{
+ if (ref->R()>hmpidIn->R()) hmpidOut = ref;
+ }
+ }
+
+
+ //EMCAL
+ if (ref->DetectorId()==AliTrackReference::kEMCAL){
+ if (!emcalIn) {
+ emcalIn = ref;
+ }else{
+ if (ref->R()>emcalIn->R()) emcalOut = ref;
+ }
+ }
+
if (ref->R()<rmin) rmin=ref->R();
if (ref->R()>rmax) rmax=ref->R();
+ } // end ref loop
+
+ // -----------------------------------------------
+ // check for gammas
+
+ // electron
+ if ( TMath::Abs(particle->GetPdgCode()) == 11 ) {
+
+ // findable
+ if (IsFindable(ipart,70)) {
+
+ // is from gamma conversion
+ Int_t motherId = particle->GetFirstMother();
+ if (motherId > 0) {
+ if (motherId < npart) {
+ TParticle* mother = (fMCinfo->Stack())->Particle(motherId);
+ if (mother && TMath::Abs(mother->GetPdgCode()) == 22) {
+ Int_t nDaughters = mother->GetNDaughters();
+
+ for (Int_t idx=0; idx<nDaughters; ++idx) {
+ Int_t daughterId = mother->GetDaughter(idx);
+ if ( daughterId == ipart || daughterId >= npart )
+ continue;
+
+ TParticle* daughter = (fMCinfo->Stack())->Particle(daughterId);
+ if (daughter && TMath::Abs(daughter->GetPdgCode()) == 11) {
+ //Bool_t findable = IsFindable(daughterId,70);
+#if USE_STREAMER
+ TTreeSRedirector *pcstream = GetDebugStreamer();
+ if (pcstream){
+ (*pcstream)<<"MCgamma"<<
+ "triggerP="<<particle<< // trigger electron
+ "motherP="<<mother<< // mother gamma
+ "daughterP="<<daughter<< // daughter electron
+ "isFindable="<<findable<< // 2nd is findable
+ "\n";
+ }
+#endif
+ }
+ }
+ }
+ }
+ }
+ }
}
+
+ // -----------------------------------------------
if (tpcIn && tpcOut) {
ProcessRefTracker(tpcIn,tpcOut,particle,1);
ProcessRefTracker(tpcIn,tpcOut,particle,3);
}
- if (itsIn && itsOut) ProcessRefTracker(itsIn,itsOut,particle,0);
- if (trdIn && trdOut) ProcessRefTracker(trdIn,trdOut,particle,2);
+ if (itsIn && itsOut) ProcessRefTracker(itsIn, itsOut, particle, 0);
+ if (trdIn && trdOut) ProcessRefTracker(trdIn, trdOut, particle, 2);
+
+ if (tpcOut && trdIn) ProcessRefTracker(tpcOut,trdIn, particle, 4);
+ if (tpcOut && tofIn) ProcessRefTracker(tpcOut,tofIn, particle, 5);
+ if (tpcOut && hmpidIn) ProcessRefTracker(tpcOut,hmpidIn,particle, 6);
+ if (tpcOut && emcalIn) ProcessRefTracker(tpcOut,emcalIn,particle, 7);
+
+ if (tpcOut && trdOut) ProcessRefTracker(tpcOut,trdOut, particle, 8);
+ if (trdIn && tofIn) ProcessRefTracker(trdIn, tofIn, particle, 9);
+ if (tofIn && tofOut) ProcessRefTracker(tofIn, tofOut, particle,10);
+
+
+ // -----------------------------------------------
+ //Test for new base tracking class
+
+ AliTrackComparison *cObj=0;
+ AliTrackPoint *point=new AliTrackPoint();
+ AliESDCaloCluster *cluster=0;
+ AliExternalTrackParam *tpcOuter=0;
+
+ Double_t eMax=0;
+ Int_t clsIndex=-1;
+
+ for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
+ {
+ cluster = fESD->GetCaloCluster(iCluster);
+ if(!cluster->IsEMCAL()) continue;
+ if(cluster->GetLabel()==ipart && cluster->E()>eMax)
+ {
+ clsIndex=iCluster;
+ eMax=cluster->E();
+ }
+ }
+
+ if(clsIndex>-1)
+ {
+ if(cluster) cluster=0;
+ cluster = fESD->GetCaloCluster(clsIndex);
+ Float_t clusterPos[3];
+ cluster->GetPosition(clusterPos);
+ point->SetXYZ(clusterPos[0],clusterPos[1],clusterPos[2],0);
+ }
+
+ if(tpcOut)
+ tpcOuter = MakeTrack(tpcOut,particle);
+
+
+ Double_t mass = particle->GetMass();
+ Int_t charge = TMath::Nint(particle->GetPDG()->Charge()/3.);
+ fPitList->Reset();
+ while(( cObj = (AliTrackComparison *)fPitList->Next()) != NULL) {
+ TString objName(cObj->GetName());
+ if(!objName.CompareTo("TPCOutToEMCalInElecCls"))
+ {
+ if(TMath::Abs(particle->GetPdgCode())==11 && tpcOuter && point && cluster && emcalIn)
+ {
+ printf("\n\nTPCOutToEMCalInElecCls: ");
+ cout<<cObj->AddTracks(tpcOuter,point,mass,cluster->E(),vPos)<<endl;
+ }
+ }
+
+ if(!objName.CompareTo("TPCOutToEMCalInElec"))
+ {
+ if(TMath::Abs(particle->GetPdgCode())==11 && tpcOut && emcalIn)
+ {
+ printf("TPCOutToEMCalInElec: ");
+ cout<<cObj->AddTracks(tpcOut,emcalIn,mass,charge)<<endl;
+ }
+ }
+
+ if(!objName.CompareTo("TPCOutToEMCalInPion"))
+ {
+ if(TMath::Abs(particle->GetPdgCode())==211 && tpcOut && emcalIn)
+ cObj->AddTracks(tpcOut,emcalIn,mass,charge);
+ }
+
+ if(!objName.CompareTo("TPCOutToTOFIn"))
+ {
+ if(tpcOut && tofIn)
+ cObj->AddTracks(tpcOut,tofIn,mass,charge);
+ }
+
+ if(!objName.CompareTo("TPCOutToHMPIDIn"))
+ {
+ if(tpcOut && hmpidIn)
+ cObj->AddTracks(tpcOut,hmpidIn,mass,charge);
+ }
+ }
+ //End of the test for new base tracking class
+ // -----------------------------------------------
+ delete point;
+
}
-}
+ trefs->Clear("C");
+ //delete particle;
+ //delete tpcIn;
+}
void AliMCTrackingTestTask::ProcessRefTracker(AliTrackReference* refIn, AliTrackReference* refOut, TParticle*part,Int_t type){
//
// Test propagation from In to out
//
+
+#if DEBUG
+ printf("ProcessRefTracker\n");
+#endif
+
AliExternalTrackParam *param = 0;
AliExternalTrackParam *paramMC = 0;
- Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
+ AliExternalTrackParam *paramDebug = 0;
+
+ // Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
+ Double_t xyzOut[3]={refOut->X(),refOut->Y(), refOut->Z()};
Double_t mass = part->GetMass();
Double_t step=1;
//
- param=MakeTrack(refOut,part);
+ param=MakeTrack(refIn,part);
paramMC=MakeTrack(refOut,part);
+ paramDebug=MakeTrack(refIn,part);
+
if (!param) return;
- if (type<3) PropagateToPoint(param,xyzIn, mass, step);
- if (type==3) {
+ if (type!=3) PropagateToPoint(param,xyzOut, mass, step);
+ //
+#if 0
+ /*
+ if (type==3) {
AliTPCseed seed;
seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
- if(seed.Rotate(alpha-seed.GetAlpha())==kFALSE) return;
+ seed.Rotate(alpha-seed.GetAlpha());
seed.SetMass(mass);
- for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
+ for (Float_t xlayer= seed.GetX(); xlayer<refOut->R(); xlayer+=step){
seed.PropagateTo(xlayer);
}
- seed.PropagateTo(refIn->R());
+ seed.PropagateTo(refOut->R());
param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
}
+*/
+#endif
+#if USE_STREAMER
TTreeSRedirector *pcstream = GetDebugStreamer();
TVectorD gpos(3);
TVectorD gmom(3);
+ Bool_t isOK=kTRUE;
+ isOK&=param->Rotate(paramMC->GetAlpha());
+ isOK&=param->PropagateTo(paramMC->GetX(),AliTracker::GetBz());
param->GetXYZ(gpos.GetMatrixArray());
param->GetPxPyPz(gmom.GetMatrixArray());
if (pcstream){
(*pcstream)<<"MC"<<
- "type="<<type<<
- "step="<<step<<
- "refIn.="<<refIn<<
- "refOut.="<<refOut<<
- "p.="<<part<<
- "par.="<<param<<
- "parMC.="<<paramMC<<
- "gpos.="<<&gpos<<
- "gmom.="<<&gmom<<
+ "isOK="<<isOK<<
+ "type="<<type<< // detector matching type
+ "step="<<step<< // propagation step length
+ "refIn.="<<refIn<< // starting track refernce
+ "refOut.="<<refOut<< // outer track reference
+ "p.="<<part<< // particle desription (TParticle)
+ "par.="<<param<< // AliExternalTrackParam create at starting point propagated to outer track ref radius
+ "parMC.="<<paramMC<< // AliExternalTrackParam created at the outer point
+ "gpos.="<<&gpos<< // global position
+ "gmom.="<<&gmom<< // global momenta
"\n";
}
+#endif
+ delete param;
+ delete paramMC;
+ delete paramDebug;
+}
+
+
+Bool_t AliMCTrackingTestTask::IsFindable(Int_t label, Float_t minTrackLength ) {
+ //
+ // Find findable tracks
+ //
+
+ AliMCParticle *mcParticle = (AliMCParticle*) fMCinfo->GetTrack(label);
+ if(!mcParticle) return kFALSE;
+
+ Int_t counter;
+ Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(),0.05,counter,3.0);
+ //printf("tpcTrackLength %f \n", tpcTrackLength);
+
+ return (tpcTrackLength>minTrackLength);
}
//
//
//
+
+#if DEBUG
+ printf("FitTrackRefs\n");
+#endif
+
const Int_t kMinRefs=6;
Int_t nrefs = trefs->GetEntries();
if (nrefs<kMinRefs) return; // we should have enough references
covar[14]=1;
AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
- AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
+ //AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
AliExternalTrackParam *paramUpdate = MakeTrack(refIn,part);
paramUpdate->AddCovariance(covar);
- Double_t mass = part->GetMass();
- Double_t charge = part->GetPDG()->Charge()/3.;
+ //Double_t mass = part->GetMass();
+ //Double_t charge = part->GetPDG()->Charge()/3.;
/*
Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
Float_t radiusIn= refIn->R();
Double_t clcov[3] = { 0.005,0,0.005};
isOKU&= paramUpdate->Update(clpos, clcov);
}
+#if USE_STREAMER
TTreeSRedirector *pcstream = GetDebugStreamer();
if (pcstream){
TVectorD gposU(3);
"gmomP.="<<&gmomP<<
"\n";
}
+#endif
+ delete paramPropagate;
+ delete paramUpdate;
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff