Misalign method is public from now on.

[u/mrichter/AliRoot.git] / STEER / AliMCParticle.cxx

diff --git a/STEER/AliMCParticle.cxx b/STEER/AliMCParticle.cxx
index 4bf7f7f0e5f62076b54542ad7b2ed86415c7c240..aff6710d20d077e5784fbfefb6a5a81c2c987145 100644 (file)
--- a/STEER/AliMCParticle.cxx
+++ b/STEER/AliMCParticle.cxx
@@ -21,37 +21,148 @@
 //     Author: Andreas Morsch, CERN
 //-------------------------------------------------------------------------
 
+#include <TRefArray.h>
+
 #include "AliMCParticle.h"
+#include "AliExternalTrackParam.h" 
 
 
 ClassImp(AliMCParticle)
 
 AliMCParticle::AliMCParticle():
     AliVParticle(),
-    fParticle(0)
+    fParticle(0),
+    fTrackReferences(0),
+    fNTrackRef(0),
+    fLabel(-1),
+    fMother(-1),
+    fFirstDaughter(-1),
+    fLastDaughter(-1)
 {
     // Constructor
 }
 
-AliMCParticle::AliMCParticle(TParticle* part):
+    
+AliMCParticle::AliMCParticle(TParticle* part, TRefArray* rarray, Int_t index):
     AliVParticle(),
-    fParticle(part)
+    fParticle(part),
+    fTrackReferences(rarray),
+    fNTrackRef(0),
+    fLabel(index),
+    fMother(-1),
+    fFirstDaughter(-1),
+    fLastDaughter(-1)
 {
     // Constructor
+    if (rarray != 0) {
+      fNTrackRef = fTrackReferences->GetEntriesFast();
+    }
 }
     
     
 AliMCParticle::AliMCParticle(const AliMCParticle& mcPart) :
     AliVParticle(mcPart),
-    fParticle(0)
+    fParticle(0),    
+    fTrackReferences(0),
+    fNTrackRef(0),
+    fLabel(-1),
+    fMother(-1),
+    fFirstDaughter(-1),
+    fLastDaughter(-1)
 {
 // Copy constructor
 }
 
 AliMCParticle& AliMCParticle::operator=(const AliMCParticle& mcPart)
-{ if (this!=&mcPart) { 
-    AliVParticle::operator=(mcPart); 
+{ 
+
+  if (this!=&mcPart) { 
+    AliVParticle::operator=(mcPart);
   }
   
   return *this; 
 }
+
+AliMCParticle::~AliMCParticle()
+{ 
+  // delete the track references passed externally
+  // fParticle should be handled by the user
+  // AliStack in case of AliMCEventHandler
+  if(fTrackReferences){
+      fTrackReferences->Delete();
+      delete fTrackReferences;
+      fTrackReferences = 0;
+  }
+}
+
+
+
+Float_t  AliMCParticle::GetTPCTrackLength(Float_t bz,  Float_t ptmin, Int_t &counter, Float_t deadWidth){
+  //
+  // return track length in geometrically active volume of TPC.
+  // z nad rphi acceptance is included
+  // doesn't take into account dead channel and ExB  
+  // Intput:
+  // trackRefs
+  // bz - magnetic field
+  // deadWidth - dead zone in r-phi
+  // Additional output:
+  // counter   - number of circles
+ 
+    const Float_t kRMin = 90;
+    const Float_t kRMax = 245;
+    const Float_t kZMax = 250;
+    const Float_t kMinPt= ptmin; 
+    Float_t length =0;
+    Int_t nrefs = fNTrackRef;
+    AliExternalTrackParam param;
+    Double_t cv[21];
+    for (Int_t i = 0; i < 21; i++) cv[i]=0;
+    counter=0; 
+    //
+    //
+    AliTrackReference *ref0 = (AliTrackReference*) (fTrackReferences->At(0));
+    Float_t direction = 0;
+  //
+    for (Int_t iref = 1; iref < nrefs; iref++){
+       AliTrackReference *ref = (AliTrackReference*) (fTrackReferences->At(iref));
+       if (!ref) continue;
+       if (!ref0 || ref0->DetectorId()!= AliTrackReference::kTPC){
+           ref0 = ref;
+           direction = ((ref0->X() * ref0->Px() + ref0->Y() * ref0->Py()) > 0)? 1. : -1.;
+           continue;
+       }
+       
+       Float_t newdirection = ((ref->X() * ref->Px() + ref->Y() * ref->Py()) > 0)? 1. : -1.;
+       if (newdirection*direction<0) {
+           counter++;  //circle counter 
+           direction = newdirection;
+           continue;
+       }
+       if (counter>0) continue;
+       if (ref0->Pt() < kMinPt) break;
+       Float_t radius0 = TMath::Max(TMath::Min(ref0->R(),kRMax),kRMin);;
+       Float_t radius1 = TMath::Max(TMath::Min(ref->R(),kRMax),kRMin);
+       Double_t xyz[3] = {ref0->X(), ref0->Y(), ref0->Z()};
+       Double_t pxyz[3]= {ref0->Px(), ref0->Py(), ref0->Pz()};
+       Double_t alpha;
+       param.Set(xyz,pxyz,cv,TMath::Nint(fParticle->GetPDG()->Charge()/3.));
+       
+       for (Float_t radius = radius0; radius < radius1; radius+=1){
+           param.GetXYZAt(radius, bz, xyz);
+           if (TMath::Abs(xyz[2]) > kZMax) continue;
+           Float_t gradius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
+           if (gradius > kRMax) continue;
+           alpha = TMath::ATan2(xyz[1],xyz[0]);
+           if (alpha<0) alpha += TMath::TwoPi();
+           //
+           Int_t sector   = Int_t(9 * alpha / TMath::Pi());
+           Float_t lalpha = alpha - ((sector + 0.5) * TMath::Pi() / 9.);
+           Float_t dedge  = (TMath::Tan(TMath::Pi() / 18.) - TMath::Abs(TMath::Tan(lalpha))) * gradius; 
+           if (dedge>deadWidth) length++;
+       }
+       if (ref->DetectorId()!= AliTrackReference::kTPC) break; 
+       ref0 = ref;
+    }
+    return length;
+}