Adding new AnalysisTask for evaluation of meterial budget
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 10 Aug 2009 08:05:32 +0000 (08:05 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 10 Aug 2009 08:05:32 +0000 (08:05 +0000)
For the moment only MC part implemented (test of principle)
To be extended

PWG1/AliMaterialBudget.cxx [new file with mode: 0644]
PWG1/AliMaterialBudget.h [new file with mode: 0644]
PWG1/PWG1LinkDef.h
PWG1/libPWG1.pkg

diff --git a/PWG1/AliMaterialBudget.cxx b/PWG1/AliMaterialBudget.cxx
new file mode 100644 (file)
index 0000000..6152c96
--- /dev/null
@@ -0,0 +1,1052 @@
+
+
+
+
+//
+//
+
+// This class estiamtes the material budget of the inner detectors in ALICE based
+// on the "upper/lower track matching"-method of the ALICE TPC.
+
+//
+//  
+//
+
+//
+
+
+
+
+// ROOT includes
+#include <TChain.h>
+#include <TMath.h>
+#include <TVectorD.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include "TGeoGlobalMagField.h"
+
+// ALIROOT includes
+#include <TTreeStream.h>
+#include <AliAnalysisManager.h>
+#include <AliESDInputHandler.h>
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliMCEventHandler.h"
+
+#include <AliESD.h>
+#include "AliMaterialBudget.h"
+#include "AliGenInfoMaker.h"
+#include "AliHelix.h"
+
+//
+#include "AliMCInfo.h"
+#include "AliComparisonObject.h"
+#include "AliESDRecInfo.h"
+#include "AliTPCParamSR.h"
+#include "AliTracker.h"
+#include "AliTPCseed.h"
+
+// STL includes
+#include <iostream>
+
+using namespace std;
+
+ClassImp(AliMaterialBudget)
+
+//________________________________________________________________________
+AliMaterialBudget::AliMaterialBudget() : 
+  AliAnalysisTask(), 
+  fMCinfo(0),     //! MC event handler
+  fESD(0),
+  fDebugStreamer(0),
+  fStreamLevel(0),
+  fDebugLevel(0),
+  fDebugOutputPath(),
+  fListHist(0),
+  fHistMult(0),
+  fCutMaxD(5),        // maximal distance in rfi ditection
+  fCutMaxDz(40),      // maximal distance in z ditection
+  fCutTheta(0.03),    // maximal distan theta
+  fCutMinDir(-0.99)   // direction vector products
+{
+  //
+  // Default constructor (should not be used)
+  //
+}
+
+AliMaterialBudget::AliMaterialBudget(const AliMaterialBudget& /*info*/) : 
+  AliAnalysisTask(), 
+  fMCinfo(0),     //! MC event handler
+  fESD(0),
+  //
+  fDebugStreamer(0),
+  fStreamLevel(0),
+  fDebugLevel(),
+  fDebugOutputPath(), 
+  fListHist(0),
+  fHistMult(0),
+  fCutMaxD(5),        // maximal distance in rfi ditection
+  fCutMaxDz(40),      // maximal distance in z ditection
+  fCutTheta(0.03),    // maximal distan theta
+  fCutMinDir(-0.99)   // direction vector products
+{
+  //
+  // Default constructor 
+  //
+}
+
+
+
+//________________________________________________________________________
+AliMaterialBudget::AliMaterialBudget(const char *name) : 
+  AliAnalysisTask(name, "AliMaterialBudget"), 
+  fMCinfo(0),     //! MC event handler
+  fESD(0),
+  fDebugStreamer(0),
+  fStreamLevel(0),
+  fDebugLevel(0),
+  fDebugOutputPath(),
+  fListHist(0),
+  fHistMult(0),
+  fCutMaxD(5),        // maximal distance in rfi ditection
+  fCutMaxDz(40),      // maximal distance in z ditection
+  fCutTheta(0.03),    // maximal distan theta
+  fCutMinDir(-0.99)   // direction vector products
+{
+  //
+  // Normal constructor
+  //
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #0 writes into a TList
+  DefineOutput(0, TList::Class());
+  //
+  //
+}
+
+AliMaterialBudget::~AliMaterialBudget(){
+  //
+  //
+  //
+  if (fDebugLevel>0)  printf("AliMaterialBudget::~AliMaterialBudget\n");
+  if (fDebugStreamer) delete fDebugStreamer;
+  fDebugStreamer=0;
+}
+
+
+//________________________________________________________________________
+void AliMaterialBudget::ConnectInputData(Option_t *) 
+{
+  //
+  // Connect the input data
+  //
+  if(fDebugLevel>3)
+    cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
+
+  TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
+  if (!tree) {
+    //Printf("ERROR: Could not read chain from input slot 0");
+  }
+  else {
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    if (!esdH) {
+      //Printf("ERROR: Could not get ESDInputHandler");
+    }
+    else {
+      esdH->SetActiveBranches("ESDfriend");
+      fESD = esdH->GetEvent();
+      //Printf("*** CONNECTED NEW EVENT ****");
+    }  
+  }
+  AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());  
+  mcinfo->SetReadTR(kTRUE);
+  
+  fMCinfo = mcinfo->MCEvent();
+
+
+}
+
+
+
+
+
+
+//________________________________________________________________________
+void AliMaterialBudget::CreateOutputObjects() 
+{
+  //
+  // Connect the output objects
+  //
+  if(fDebugLevel>3)
+    cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
+  //
+  fListHist = new TList();
+  fListHist->SetOwner();
+  //
+  fHistMult = new TH1F("HistMult", "Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
+  fListHist->Add(fHistMult);
+
+
+}
+
+
+//________________________________________________________________________
+void AliMaterialBudget::Exec(Option_t *) {
+  //
+  // Execute analysis for current event 
+  //
+
+  if(fDebugLevel>3)
+    cout << "AliMaterialBudget::Exec()" << endl;
+
+  fHistMult->Fill(fESD->GetNumberOfTracks());
+  //FindPairs(fESD); // nearly everything takes place in find pairs...
+
+  // If MC has been connected   
+
+  if (!fMCinfo){
+    cout << "Not MC info\n" << endl;
+  }else{
+    ProcessMCInfo();
+    //mcinfo->Print();
+    //DumpInfo();
+  }
+  //
+  PostData(0, fListHist);
+}      
+
+
+
+
+//________________________________________________________________________
+void AliMaterialBudget::Terminate(Option_t *) {
+    //
+    // Terminate loop
+    //
+  if(fDebugLevel>3)
+    printf("AliMaterialBudget: Terminate() \n");  
+  //
+  if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
+  if (fDebugStreamer) delete fDebugStreamer;
+  fDebugStreamer = 0;
+  return;
+}
+
+
+
+TTreeSRedirector *AliMaterialBudget::GetDebugStreamer(){
+  //
+  // Get Debug streamer
+  // In case debug streamer not yet initialized and StreamLevel>0 create new one
+  //
+  if (fStreamLevel==0) return 0;
+  if (fDebugStreamer) return fDebugStreamer;
+  TString dsName;
+  dsName=GetName();
+  dsName+="Debug.root";
+  dsName.ReplaceAll(" ","");
+  fDebugStreamer = new TTreeSRedirector(dsName.Data());
+  return fDebugStreamer;
+}
+
+
+
+
+
+
+AliExternalTrackParam * AliMaterialBudget::MakeTrack(const AliTrackReference* ref, TParticle*part)
+{
+  //
+  // Make track out of the track ref
+  // part - TParticle used to determine chargr
+  // the covariance matrix - equal 0 - starting from ideal MC position
+  if (!part->GetPDG()) return 0;
+  Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
+  Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
+  Double_t charge = TMath::Nint(part->GetPDG()->Charge()/3.);
+  if (ref->X()*ref->Px()+ref->Y()*ref->Py() <0){
+    pxyz[0]*=-1;
+    pxyz[1]*=-1;
+    pxyz[2]*=-1;
+    charge*=-1.;
+  }
+  Double_t cv[21];
+  for (Int_t i=0; i<21;i++) cv[i]=0;
+  AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,charge);
+  return param;
+}
+
+Bool_t  AliMaterialBudget::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
+  // 
+  // Propagate track to point xyz using 
+  // AliTracker::PropagateTo 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::PropagateTrackToBxByBz(param, radius+step, mass, step, kTRUE,0.99);
+  AliTracker::PropagateTrackToBxByBz(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);
+  return isOK;
+}
+
+
+void  AliMaterialBudget::ProcessMCInfo(){
+  //
+  //
+  //
+  //
+  TParticle * particle= new TParticle;
+  TClonesArray * trefs = new TClonesArray("AliTrackReference");
+  const Double_t kPcut=0.05;
+  const Double_t kMinDrITS = 2.;   // minimal distance between references
+  const Double_t kMinDrTRD = 8.;   // minimal distance between references
+  const Double_t kMinDrTOF = 10.;   // minimal distance between references
+  //
+  //
+  // Process tracks
+  //
+  Int_t npart = fMCinfo->GetNumberOfTracks();
+  if (npart==0) return;
+  Double_t vertex[4]={0,0,0,0};
+  fMCinfo->GetParticleAndTR(0, particle, trefs);
+  if (particle){
+    vertex[0]=particle->Vx();
+    vertex[1]=particle->Vy();
+    vertex[2]=particle->Vz();
+    vertex[3]=particle->R();
+  }
+  //
+  //
+  AliTrackReference dummy,*pdummy= &dummy;
+  AliExternalTrackParam epdummy,*pepdummy= &epdummy;
+  Int_t nRefITS =0;
+  Int_t nRefTPC =0;
+  Int_t nRefTRD =0;
+  Int_t nRefTOF =0;
+  AliTrackReference * refITS0, *refITS1;
+  AliTrackReference * refTPC0, *refTPC1;
+  AliTrackReference * refTPCIn0, *refTPCIn1;
+  AliTrackReference * refTRD0, *refTRD1;
+  AliTrackReference * refTOF0, *refTOF1;
+  AliTrackReference *refMinR;
+  //
+  for (Int_t ipart=0;ipart<npart;ipart++){
+    Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
+    AliMCParticle * pp = fMCinfo->GetTrack(ipart);
+    if (!pp) continue;
+    if (particle->P()<kPcut) continue;
+    Double_t mass = particle->GetMass();
+
+    // RESET
+    nRefITS =0;
+    nRefTPC =0;
+    nRefTRD =0;
+    nRefTOF =0;
+    refITS0=pdummy;    refITS1=pdummy;
+    refTPC0=pdummy;    refTPC1=pdummy;
+    refTPCIn0=pdummy;    refTPCIn1=pdummy;
+    refTRD0=pdummy;    refTRD1=pdummy;
+    refTOF0=pdummy;    refTOF1=pdummy;
+    refMinR = pdummy;
+    //
+    Int_t nref = pp->GetNumberOfTrackReferences();
+    if (nref==0) continue;
+    for (Int_t iref = 0; iref < nref; iref++) { 
+      AliTrackReference *ref = pp->GetTrackReference(iref);
+      if (ref->DetectorId()==AliTrackReference::kDisappeared) continue;      
+      //      if (ref.Px()*particle.Px()+ref.Py()*particle.Py()<0) break; // returning track
+      //
+      if (ref->DetectorId()==AliTrackReference::kITS){
+       if (TMath::Abs(ref->R()-refITS1->R())>kMinDrITS) {
+         refITS1 = ref;
+         nRefITS++;
+       }
+       if (refITS0==pdummy) refITS0=ref;
+      }
+      if (ref->DetectorId()==AliTrackReference::kTPC){ 
+       nRefTPC++;
+       refTPC1 = ref;
+       if (refTPC0==pdummy) refTPC0=ref;
+      }
+      if (ref->DetectorId()==AliTrackReference::kTRD){
+       if (TMath::Abs(ref->R()-refTRD1->R())>kMinDrTRD) {
+         refTRD1 = ref;
+         nRefTRD++;
+       }
+       if (refTRD0==pdummy) refTRD0=ref;
+      }
+      if (ref->DetectorId()==AliTrackReference::kTOF){
+       if (TMath::Abs(ref->X()-refTOF1->X()) + TMath::Abs(ref->Y()-refTOF1->Y()) + TMath::Abs(ref->Z()-refTOF1->Z())>kMinDrTOF) {
+         refTOF1 = ref;
+         nRefTOF++;
+       }
+       if (refTOF0==pdummy) refTOF0=ref;
+      }
+      //
+      // "find inner track ref"
+      if (ref->DetectorId()==AliTrackReference::kTPC){ 
+       if (ref->Px()*ref->X()+ref->Py()*ref->Y()<0){
+         //  track in
+         if (refTPCIn0 == pdummy) refTPCIn0=ref;
+         if (refTPCIn0 != pdummy &&  refTPCIn0->R()>ref->R())
+           refTPCIn0=ref; 
+       }
+       if (ref->Px()*ref->X()+ref->Py()*ref->Y()>0){
+         //  track in
+         if (refTPCIn1 == pdummy) refTPCIn1=ref;
+         if (refTPCIn1 != pdummy &&  refTPCIn1->R()>ref->R())
+           refTPCIn1=ref; 
+       }
+      }
+
+
+      if (refMinR==pdummy && ref->P()>0  ){
+       refMinR=ref;
+      }
+      if (refMinR->R()>ref->R() && ref->P()>0 ) refMinR=ref;
+      
+    }
+    //
+    AliExternalTrackParam * trackMC = pepdummy;
+    //track0->GetDZ(0,0,0,bz,dvertex0)
+    Float_t dist[2]={0,0};
+    AliMagF* field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
+    Double_t esdfield= fESD->GetMagneticField();
+    Double_t xyz[3]={0,0,0};
+    Double_t bxyz[3]={0,0,0};
+    field->Field(xyz,bxyz);
+    if (refMinR->P()>0) {
+      trackMC = MakeTrack(refMinR,particle); 
+      trackMC->GetDZ(0,0,0,bxyz[2],dist);
+    }
+    Double_t alphaTOF0 = TMath::ATan2(refTOF0->Y(),refTOF0->X());
+    Double_t alphaTOF1 = TMath::ATan2(refTOF1->Y(),refTOF1->X());
+    Int_t dsecTOF   = TMath::Nint(180*(alphaTOF0-alphaTOF1)/(TMath::Pi()*20.)-0.5);
+    //
+    // make the two different TPC tracks and propagate them to their DCA
+    //
+    Double_t dP = 0;
+    Bool_t statusProp = kFALSE;
+    Double_t dY = 0;
+    Double_t dZ = 0;
+    AliExternalTrackParam * track0  = pepdummy;
+    AliExternalTrackParam * track1  = pepdummy;
+    AliExternalTrackParam * otrack0 = pepdummy;
+    AliExternalTrackParam * otrack1 = pepdummy;
+    if (refTPCIn0!=pdummy && refTPCIn1!=pdummy) {
+      track0 = MakeTrack(refTPCIn0,particle); 
+      track1 = MakeTrack(refTPCIn1,particle);
+      otrack0 = MakeTrack(refTPCIn0,particle); 
+      otrack1 = MakeTrack(refTPCIn1,particle);
+      dP = track0->P() - track1->P(); // momentum loss
+      statusProp = AliMaterialBudget::PropagateCosmicToDCA(track0,track1,mass);
+      if (statusProp) {
+       dY = track0->GetY() - track1->GetY();
+       dZ = track0->GetZ() - track1->GetZ();
+      }
+    }
+    //
+    TTreeSRedirector *pcstream = GetDebugStreamer();
+    if (pcstream){      
+      char name[100];
+      for (Int_t id=0; id<3;id++){
+       
+       //      if (id==0) sprintf(name,"mcAll"); // all tracks: inconvenient to cut if on is only interest in tracks which reach the TPC
+       if (id==0) continue; // require TPC
+       if (id==1) sprintf(name,"mcITS");
+       if (id==2) sprintf(name,"mcTPC");
+       if (id==1&& nRefITS==0) continue;
+       if (id==2&& nRefTPC==0) continue;
+
+       (*pcstream)<<name<<
+         "ipart="<<ipart<<
+         "p.="<<particle<<
+         "mass="<<mass<<
+         "tbfield="<<bxyz[2]<<
+         "esdbfield="<<esdfield<<
+         // counter
+         "nref="<<nref<<
+         "nRefITS="<<nRefITS<<     
+         "nRefTPC="<<nRefTPC<<
+         "nRefTRD="<<nRefTRD<<
+         "nRefTOF="<<nRefTOF<<
+         //references
+         "refMinR.="<<refMinR<<
+         "refITS0.="<<refITS0<<
+         "refITS1.="<<refITS1<<
+         "refTPC0.="<<refTPC0<<
+         "refTPC1.="<<refTPC1<<
+         "refTPCIn0.="<<refTPCIn0<<
+         "refTPCIn1.="<<refTPCIn1<<
+         "refTRD0.="<<refTRD0<<
+         "refTRD1.="<<refTRD1<<
+         "refTOF0.="<<refTOF0<<
+         "refTOF1.="<<refTOF1<<
+         //trigger variables
+         "dsecTOF="<<dsecTOF<<    // delta TOF sectors
+         "aTOF0="<<alphaTOF0<<
+         "aTOF1="<<alphaTOF1<<
+         //
+         // track
+         "dr="<<dist[0]<<
+         "dz="<<dist[1]<<
+         //
+         // "two" TPC tracks 
+         "status="<<statusProp<<
+         "dP="<<dP<<
+         "otrack0.="<<otrack0<<
+         "otrack1.="<<otrack1<<
+         "track0.="<<track0<<
+         "track1.="<<track1<<
+         "dY="<<dY<<
+         "dZ="<<dZ<<
+         "\n";
+      }
+    }    
+  }
+}
+
+
+
+void AliMaterialBudget::ProcessRefTracker(AliTrackReference* refIn,  AliTrackReference* refOut, TParticle*part,Int_t type){
+  //
+  // Test propagation from In to out
+  //
+  AliExternalTrackParam *param = 0;
+  AliExternalTrackParam *paramMC = 0;
+  Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
+  Double_t mass = part->GetMass();
+  Double_t step=1;
+  //
+  param=MakeTrack(refOut,part);
+  paramMC=MakeTrack(refOut,part);
+  if (!param) return;
+  if (type<3) PropagateToPoint(param,xyzIn, mass, step);
+  if (type==3) {
+    AliTPCseed seed;
+    seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
+    Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
+    seed.Rotate(alpha-seed.GetAlpha());
+    seed.SetMass(mass);
+    for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
+      seed.PropagateTo(xlayer);
+    }
+    seed.PropagateTo(refIn->R());
+    param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
+  }
+  TTreeSRedirector *pcstream = GetDebugStreamer();
+  TVectorD gpos(3);
+  TVectorD gmom(3);
+  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<<
+      "\n";
+  }
+}
+
+
+void  AliMaterialBudget::FitTrackRefs(TParticle * part, TClonesArray * trefs){
+  //
+  //
+  //
+  //
+  const Int_t kMinRefs=6;
+  Int_t nrefs = trefs->GetEntries();
+  if (nrefs<kMinRefs) return; // we should have enough references
+  Int_t iref0 =-1;
+  Int_t iref1 =-1;
+  
+  for (Int_t iref=0; iref<nrefs; iref++){
+    AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
+    if (!ref) continue;    
+    Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
+    if (dir<0) break;
+    if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
+    if (iref0<0) iref0 = iref;
+    iref1 = iref;    
+  }
+  if (iref1-iref0<kMinRefs) return;
+  Double_t covar[14];
+  for (Int_t icov=0; icov<14; icov++) covar[icov]=0;
+  covar[0]=1; 
+  covar[2]=1; 
+  covar[5]=1;
+  covar[9]=1;
+  covar[14]=1;
+
+  AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
+  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.;
+/*
+  Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
+  Float_t radiusIn= refIn->R();
+  Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
+  Float_t radiusOut= refOut->R();
+*/
+  Bool_t isOKP=kTRUE;
+  Bool_t isOKU=kTRUE;
+  AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
+  for (Int_t iref = iref0; iref<=iref1; iref++){
+    AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
+    Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
+    Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
+    Double_t mag[3];
+    field->Field(pos,mag);
+    isOKP&=paramPropagate->Rotate(alphaC);
+    isOKU&=paramUpdate->Rotate(alphaC);
+    for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
+      isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
+      isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
+    }
+    isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
+    isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
+    Double_t clpos[2] = {0, ref->Z()};
+    Double_t clcov[3] = { 0.005,0,0.005};
+    isOKU&= paramUpdate->Update(clpos, clcov);  
+  }
+  TTreeSRedirector *pcstream = GetDebugStreamer();
+  if (pcstream){
+    TVectorD gposU(3);
+    TVectorD gmomU(3);
+    TVectorD gposP(3);
+    TVectorD gmomP(3);
+    paramUpdate->GetXYZ(gposU.GetMatrixArray());
+    paramUpdate->GetPxPyPz(gmomU.GetMatrixArray());
+    paramPropagate->GetXYZ(gposP.GetMatrixArray());
+    paramPropagate->GetPxPyPz(gmomP.GetMatrixArray());
+
+     (*pcstream)<<"MCupdate"<<
+       "isOKU="<<isOKU<<
+       "isOKP="<<isOKP<<
+       "m="<<mass<<
+       "q="<<charge<<
+       "part.="<<part<<
+       "refIn.="<<refIn<<
+       "refOut.="<<refOut<<
+       "pP.="<<paramPropagate<<
+       "pU.="<<paramUpdate<<
+       "gposU.="<<&gposU<<
+       "gmomU.="<<&gmomU<<
+       "gposP.="<<&gposP<<
+       "gmomP.="<<&gmomP<<
+       "\n";
+   }
+}
+
+
+
+void AliMaterialBudget::FindPairs(AliESDEvent * event) {
+  //
+  // This function matches the cosmic tracks and calculates the energy loss in the material.
+  // If accessible the "true" energy loss is determined with the MC track references.
+  //
+  //
+  // Find cosmic pairs
+  // 
+  // Track0 is choosen in upper TPC part
+  // Track1 is choosen in lower TPC part
+  //
+  if(fDebugLevel>3)
+    cout << "AliMaterialBudget::FindPairs()" << endl;
+
+
+  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+  Int_t ntracks=event->GetNumberOfTracks(); 
+  TObjArray  tpcSeeds(ntracks);
+  if (ntracks==0) return;
+  Double_t vtxx[3]={0,0,0};
+  Double_t svtxx[3]={0.000001,0.000001,100.};
+  AliESDVertex vtx(vtxx,svtxx);
+  //
+  //track loop
+  //
+  for (Int_t i=0;i<ntracks;++i) {
+   AliESDtrack *track = event->GetTrack(i);
+     const AliExternalTrackParam * trackIn = track->GetInnerParam();
+   const AliExternalTrackParam * trackOut = track->GetOuterParam();
+   if (!trackIn) continue;
+   if (!trackOut) continue;
+   AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
+   TObject *calibObject;
+   AliTPCseed *seed = 0;
+   for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
+     if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
+   }
+   if (seed) tpcSeeds.AddAt(seed,i);
+  }
+
+  if (ntracks<2) return;
+  //
+  // Find pairs
+  //
+  for (Int_t i=0;i<ntracks;++i) {
+    AliESDtrack *track0 = event->GetTrack(i);     
+    // track0 - choosen upper part
+    if (!track0) continue;
+    if (!track0->GetOuterParam()) continue;
+    if (track0->GetOuterParam()->GetAlpha()<0) continue;
+    Double_t dir0[3];
+    track0->GetDirection(dir0);    
+    for (Int_t j=0;j<ntracks;++j) {
+      if (i==j) continue;
+      AliESDtrack *track1 = event->GetTrack(j);   
+      //track 1 lower part
+      if (!track1) continue;
+      if (!track1->GetOuterParam()) continue;
+      if (track1->GetOuterParam()->GetAlpha()>0) continue;
+      //
+      Double_t dir1[3];
+      track1->GetDirection(dir1);
+      
+      AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
+      AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
+      if (! seed0) continue;
+      if (! seed1) continue;
+      //
+      Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
+      Float_t d0  = track0->GetLinearD(0,0);
+      Float_t d1  = track1->GetLinearD(0,0);
+      //
+      // conservative cuts - convergence to be guarantied
+      // applying before track propagation
+      if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
+      if (dir>fCutMinDir) continue;               // direction vector product
+      Float_t bz = AliTracker::GetBz();
+      Float_t dvertex0[2];   //distance to 0,0
+      Float_t dvertex1[2];   //distance to 0,0 
+      track0->GetDZ(0,0,0,bz,dvertex0);
+      track1->GetDZ(0,0,0,bz,dvertex1);
+      if (TMath::Abs(dvertex0[1])>250) continue;
+      if (TMath::Abs(dvertex1[1])>250) continue;
+      //
+      //
+      //
+      Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
+      AliExternalTrackParam param0(*track0);
+      AliExternalTrackParam param1(*track1);
+      //
+      // Propagate using Magnetic field and correct fo material budget
+      //
+      AliTracker::PropagateTrackToBxByBz(&param0,dmax+1,0.0005,3,kTRUE);
+      AliTracker::PropagateTrackToBxByBz(&param1,dmax+1,0.0005,3,kTRUE);
+      //
+      // Propagate rest to the 0,0 DCA - z should be ignored
+      //
+      Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
+      Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
+      //      
+      param0.GetDZ(0,0,0,bz,dvertex0);
+      param1.GetDZ(0,0,0,bz,dvertex1);
+      if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
+      //
+      Double_t xyz0[3];//,pxyz0[3];
+      Double_t xyz1[3];//,pxyz1[3];
+      param0.GetXYZ(xyz0);
+      param1.GetXYZ(xyz1);
+      Bool_t isPair = IsPair(&param0,&param1);
+      //
+      // HERE WE WILL PUT THE ACCESS TO THE MC TRACKS AND MATCH THESE !!!!
+      //
+      Int_t label0 = TMath::Abs(track0->GetLabel());
+      AliMCParticle *mcParticle0 = fMCinfo->GetTrack(label0);
+      TParticle *particle0 = mcParticle0->Particle();
+      AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle0); // get the first TPC track reference
+      if (!ref0) continue;
+      AliExternalTrackParam *paramMC0 = 0;
+      paramMC0 = MakeTrack(ref0, particle0);
+      //
+      Int_t label1 = TMath::Abs(track1->GetLabel());
+      AliMCParticle *mcParticle1 = fMCinfo->GetTrack(label1);
+      TParticle *particle1 = mcParticle1->Particle();
+      AliTrackReference *ref1 = GetFirstTPCTrackRef(mcParticle1); // get the first TPC track reference
+      if (!ref1) continue;
+      AliExternalTrackParam *paramMC1 = 0;
+      paramMC1 = MakeTrack(ref1, particle1);
+      //
+      // ACCESS TOF INFORMATION
+      Int_t nTrackRefTOF0 = 0;
+      Int_t nTrackRefITS0 = 0;
+      AliTrackReference * refLastTOF0 = 0;
+      AliTrackReference * refFirstTOF0 = GetAllTOFinfo(mcParticle0, nTrackRefTOF0, nTrackRefITS0);
+      Float_t alphaTOF0 = 0;
+      if (refFirstTOF0) alphaTOF0 = refFirstTOF0->Alpha();
+      //
+      Int_t nTrackRefTOF1 = 0;
+      Int_t nTrackRefITS1 = 0;
+      AliTrackReference  * refLastTOF1 = 0;
+      AliTrackReference  * refFirstTOF1 =GetAllTOFinfo(mcParticle1, nTrackRefTOF1, nTrackRefITS1);
+      Float_t alphaTOF1 = 0;
+      if (refFirstTOF1) alphaTOF1 = refFirstTOF1->Alpha();
+      //cout << " STATUS "<<nTrackRefTOF0<<" "<<refFirstTOF0<<" "<<refLastTOF0<<" " <<refFirstTOF1<<" "<<refLastTOF1<<endl;
+      //
+      //
+      //
+      if (fStreamLevel>0){
+       TTreeSRedirector * cstream =  GetDebugStreamer();
+       AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
+       AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
+       AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
+       AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
+       //
+       //
+       //
+       if (cstream) {
+         (*cstream) << "Track0" <<
+           "dir="<<dir<<               //  direction
+           "OK="<<isPair<<             //  will be accepted
+           "b0="<<b0<<                 //  propagate status
+           "b1="<<b1<<                 //  propagate status
+           //
+           "Particle.="<<particle0<<             // TParticle with generated momentum
+           "NTrackRefTOF0="<<nTrackRefTOF0<<     // Number of TOF track references upper
+           "NTrackRefTOF1="<<nTrackRefTOF1<<     // Number of TOF track references lower
+           "NTrackRefITS0="<<nTrackRefITS0<<     // Number of ITS track references upper
+           "NTrackRefITS1="<<nTrackRefITS1<<     // Number of ITS track references lower
+           "Alpha0="<<alphaTOF0<<                   // alpha upper
+           "Alpha1="<<alphaTOF1<<                   // alpha lower
+           "RefFirstTOF0.="<<refFirstTOF0<<        // first tof reference upper
+           "RefLastTOF0.="<<refLastTOF0<<          // last tof reference upper
+           "RefFirstTOF1.="<<refFirstTOF1<<        // first tof reference lower
+           "RefLastTOF1.="<<refLastTOF1<<          // last tof reference lower
+           //
+           "Orig0.=" << track0 <<      //  original track  0
+           "Orig1.=" << track1 <<      //  original track  1
+           "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
+           "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0        
+           "Ip0.="<<ip0<<              //  inner param - upper
+           "Ip1.="<<ip1<<              //  inner param - lower
+           "Op0.="<<op0<<              //  outer param - upper
+           "Op1.="<<op1<<              //  outer param - lower
+           //
+           "paramTrackRef0.="<<paramMC0<< //  "ideal" MC true track parameters from track references - upper
+           "paramTrackRef1.="<<paramMC1<< //  "ideal" MC true track parameters from track references - lower
+           //
+           "v00="<<dvertex0[0]<<       //  distance using kalman
+           "v01="<<dvertex0[1]<<       // 
+           "v10="<<dvertex1[0]<<       //
+           "v11="<<dvertex1[1]<<       // 
+           "d0="<<d0<<                 //  linear distance to 0,0
+           "d1="<<d1<<                 //  linear distance to 0,0
+           //
+           //
+           //
+           "x00="<<xyz0[0]<<           // global position close to vertex
+           "x01="<<xyz0[1]<<
+           "x02="<<xyz0[2]<<
+           //
+           "x10="<<xyz1[0]<<           // global position close to vertex
+           "x11="<<xyz1[1]<<
+           "x12="<<xyz1[2]<<
+           //
+           "dir00="<<dir0[0]<<           // direction upper
+           "dir01="<<dir0[1]<<
+           "dir02="<<dir0[2]<<
+           //
+           "dir10="<<dir1[0]<<           // direction lower
+           "dir11="<<dir1[1]<<
+           "dir12="<<dir1[2]<<
+           //
+           //
+           "Seed0.=" << seed0 <<       //  original seed 0
+           "Seed1.=" << seed1 <<       //  original seed 1
+           //
+           "\n";
+       }
+      }
+      delete paramMC0;
+      delete paramMC1;
+    }
+  } 
+
+  return;
+
+
+}
+
+Bool_t  AliMaterialBudget::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
+  //
+  //
+  /*
+  // 0. Same direction - OPOSITE  - cutDir +cutT    
+  TCut cutDir("cutDir","dir<-0.99")
+  // 1. 
+  TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
+  //
+  // 2. The same rphi 
+  TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
+  //
+  //
+  //
+  TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");  
+  // 1/Pt diff cut
+  */
+  const Double_t *p0 = tr0->GetParameter();
+  const Double_t *p1 = tr1->GetParameter();
+  if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
+  if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
+  if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
+  
+  Double_t d0[3], d1[3];
+  tr0->GetDirection(d0);    
+  tr1->GetDirection(d1);       
+  if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
+  //
+  return kTRUE;  
+}
+
+AliTrackReference * AliMaterialBudget::GetFirstTPCTrackRef(AliMCParticle *mcParticle) 
+{
+  // return first TPC track reference 
+  if(!mcParticle) return 0;
+
+  // find first track reference 
+  // check direction to select proper reference point for looping tracks
+  Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
+  AliTrackReference *ref = 0;
+  AliTrackReference *refIn = 0;
+  for (Int_t iref = 0; iref < nTrackRef; iref++) { 
+    ref = mcParticle->GetTrackReference(iref);
+    if(ref && (ref->DetectorId()==AliTrackReference::kTPC))
+    {
+      //Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
+      //if(dir < 0.) break;
+
+      refIn = ref;
+      break;
+    }
+  }
+  
+  return refIn;
+}
+
+
+AliTrackReference *  AliMaterialBudget::GetAllTOFinfo(AliMCParticle *mcParticle, Int_t &nTrackRef, Int_t &nTrackRefITS, Int_t retValue) {
+  //
+  //
+  //
+
+  if(!mcParticle) return 0;
+  Int_t counter = 0;
+  nTrackRef = 0;
+  nTrackRefITS = 0;
+  AliTrackReference *ref = 0;
+  for (Int_t iref = 0; iref < mcParticle->GetNumberOfTrackReferences(); iref++) { 
+    ref = mcParticle->GetTrackReference(iref);
+    if(ref && (ref->DetectorId()==AliTrackReference::kTOF)) {
+      counter = iref;
+      nTrackRef++;
+    }
+    if(ref && (ref->DetectorId()==AliTrackReference::kITS)) nTrackRefITS++;    
+  }
+  if (nTrackRef ==0) return 0;
+  if (retValue == 0) return mcParticle->GetTrackReference(counter - nTrackRef +1);
+  return mcParticle->GetTrackReference(counter);
+}
+
+
+void AliMaterialBudget::FinishTaskOutput()
+{
+  //
+  // According description in AliAnalisysTask this method is call
+  // on the slaves before sending data
+  //
+  Terminate("slave");
+  gSystem->Exec("pwd");
+  RegisterDebugOutput();
+
+}
+
+
+void AliMaterialBudget::RegisterDebugOutput(){
+  //
+  //
+  //
+  //
+  // store  - copy debug output to the destination position
+  // currently ONLY for local copy
+  TString dsName;
+  dsName=GetName();
+  dsName+="Debug.root";
+  dsName.ReplaceAll(" ","");
+  TString dsName2=fDebugOutputPath.Data();
+  gSystem->MakeDirectory(dsName2.Data());
+  dsName2+=gSystem->HostName();
+  gSystem->MakeDirectory(dsName2.Data());
+  dsName2+="/";
+  dsName2+=gSystem->BaseName(gSystem->pwd());
+  dsName2+="/";
+  gSystem->MakeDirectory(dsName2.Data());
+  dsName2+=dsName;
+  AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
+  printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
+  TFile::Cp(dsName.Data(),dsName2.Data());
+}
+
+
+Bool_t AliMaterialBudget::PropagateCosmicToDCA(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass){
+  //
+  // param0 - upper part of cosmic track
+  // param1 - lower part of cosmic track
+  //
+  // 0. Propagate both tracks to DCA to (0,0,0)
+  // 1. After propagation to DCA rotate track param1 to corrdinate system of track1 <-> rotate param0 to coordinate system of param 1 ????
+  // 2. Propagate track 1 to refernce x from track0
+  //
+
+  // step 0.
+
+  Float_t d0  = param0->GetLinearD(0,0);
+  Float_t d1  = param1->GetLinearD(0,0);
+  Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
+  //
+  // propagate in the beginning taking all material into account
+  //
+  AliTracker::PropagateTrackToBxByBz(param0,dmax+1.,mass,0.5,kTRUE,0.99,-1.);
+  AliTracker::PropagateTrackToBxByBz(param1,dmax+1.,mass,0.5,kTRUE,0.99,1.);
+  //
+  Double_t vtxx[3]={0,0,0};
+  Double_t svtxx[3]={0.000001,0.000001,100.};
+  AliESDVertex vtx(vtxx,svtxx);
+  //
+  Bool_t b0 = param0->PropagateToDCA(&vtx,AliTracker::GetBz(),1000);
+  Bool_t b1 = param1->PropagateToDCA(&vtx,AliTracker::GetBz(),1000);
+
+  if (!(b0 && b1)) return 0;
+
+  // step 1.
+    
+  Float_t dAlpha = param0->GetAlpha();
+  param1->Rotate(dAlpha);
+
+  // step 2.
+
+  Float_t refX = param0->GetX();
+  param1->PropagateTo(refX,AliTracker::GetBz());
+
+  return kTRUE;
+  
+  
+}
+
+
diff --git a/PWG1/AliMaterialBudget.h b/PWG1/AliMaterialBudget.h
new file mode 100644 (file)
index 0000000..0fb1272
--- /dev/null
@@ -0,0 +1,93 @@
+#ifndef ALIMATERIALBUDGET_H
+#define ALIMATERIALBUDGET_H
+
+// ROOT includes
+#include <TList.h>
+#include <TH1.h>
+#include <TH2.h>
+
+// AliRoot includes
+#include <AliAnalysisTask.h>
+#include <AliESDEvent.h>
+#include <AliESDfriend.h>
+#include <AliESDtrack.h>
+#include <AliESDfriendTrack.h>
+#include <AliTPCseed.h>
+#include <TString.h>
+class AliGenInfoMaker;
+class TTreeSRedirector;
+class AliMCEventHadnler;
+class TParticle;
+class AliMCInfo;
+class AliESDRecInfo;
+class AliESDEvent;
+class AliMCEvent;
+class AliComparisonObject;
+
+class AliMaterialBudget : public AliAnalysisTask {
+ public:
+ AliMaterialBudget();
+ AliMaterialBudget(const char *name);
+  virtual ~AliMaterialBudget();
+  
+  virtual void   ConnectInputData(Option_t *);
+  virtual void   CreateOutputObjects();
+  virtual void   Exec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  virtual void FinishTaskOutput();
+  void         SetDebugOuputhPath(const char * name){fDebugOutputPath=name;}
+
+  //
+  void           FindPairs(AliESDEvent * event);
+  Bool_t         IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1);
+  //
+  void           ProcessMCInfo();
+  void           ProcessRefTracker(AliTrackReference* refIn, AliTrackReference* refOut, TParticle*part, Int_t type);
+  
+  void           FitTrackRefs(TParticle * part, TClonesArray * trefs);
+
+  //
+  // debug streamer part
+  //
+  TTreeSRedirector *GetDebugStreamer();
+  void       SetStreamLevel(Int_t streamLevel){fStreamLevel=streamLevel;}
+  void       SetDebugLevel(Int_t level) {fDebugLevel = level;}
+  Int_t      GetStreamLevel() const {return fStreamLevel;}
+  Int_t      GetDebugLevel() const {return fDebugLevel;}
+  //
+  static Bool_t PropagateCosmicToDCA(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass);
+  static AliExternalTrackParam * MakeTrack(const AliTrackReference* ref, TParticle*part);
+  static Bool_t  PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass,  Float_t step);
+  //
+  AliTrackReference * GetFirstTPCTrackRef(AliMCParticle *mcParticle);
+  AliTrackReference * GetAllTOFinfo(AliMCParticle *mcParticle, Int_t & nTrackRef, Int_t &nTrackRefITS, Int_t retValue =0);
+ protected:
+  void RegisterDebugOutput();
+  AliMaterialBudget(const AliMaterialBudget& /*info*/);
+  AliMaterialBudget& operator=(const AliMaterialBudget& /*info*/) { return *this;}
+  AliMCEvent  * fMCinfo;          //! MC event handler
+  AliESDEvent * fESD;             //! current esd event
+  //
+  //
+  //
+  TTreeSRedirector *fDebugStreamer;     //! debug streamer
+  Int_t  fStreamLevel;                  //  debug stream level 
+  Int_t  fDebugLevel;                   //  debug level
+  TString      fDebugOutputPath; // debug output path
+  //
+  // histogran
+  //
+  TList * fListHist;     // list for histograms
+  TH1F  * fHistMult;     // track multiplicity histograms
+  //
+  // cuts
+  //
+  Float_t fCutMaxD;     // maximal distance in rfi ditection
+  Float_t fCutMaxDz;    // maximal distance in z ditection
+  Float_t fCutTheta;    // maximal distance in theta ditection
+  Float_t fCutMinDir;   // direction vector products
+  //
+  ClassDef(AliMaterialBudget, 1); // Analysis task base class for tracks
+};
+
+#endif
index 8babf31..4ec06fb 100644 (file)
@@ -44,6 +44,8 @@
 
 #pragma link C++ class AliAlignmentDataFilterITS+;
 #pragma link C++ class AliAnalysisTaskV0QA+;
+#pragma link C++ class AliMaterialBudget+;
+
 
 
 #endif
index 167c7b2..f88b0ea 100644 (file)
@@ -31,7 +31,8 @@ SRCS:=        AliTreeDraw.cxx \
        AliPerformanceTPC.cxx \
        AliPerformanceMC.cxx \
        AliAlignmentDataFilterITS.cxx \
-       AliAnalysisTaskV0QA.cxx
+       AliAnalysisTaskV0QA.cxx \
+       AliMaterialBudget.cxx
 
 HDRS:= $(SRCS:.cxx=.h)