]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG1/AliGenInfoTask.cxx
adding the HLT offline QA skeleton
[u/mrichter/AliRoot.git] / PWG1 / AliGenInfoTask.cxx
index 4592a4f4f47be8cfdf8ca6da9a71e671aac82c66..72465dda9d332684be83d12628e58532c6470a64 100644 (file)
 #include "AliComparisonObject.h"
 #include "AliESDRecInfo.h"
 #include "AliTPCParamSR.h"
+#include "TSystem.h"
+#include "TTimeStamp.h"
+#include "TFile.h"
+#include "AliTPCseed.h"
 
 // STL includes
 #include <iostream>
@@ -41,6 +45,7 @@ AliGenInfoTask::AliGenInfoTask() :
   AliAnalysisTask(), 
   fMCinfo(0),     //! MC event handler
   fESD(0),
+  fESDfriend(0),
   fCompList(0),         //array of comparison objects
   fGenTracksArray(0),  //clones array with filtered particles
   fGenKinkArray(0),    //clones array with filtered Kinks
@@ -48,7 +53,8 @@ AliGenInfoTask::AliGenInfoTask() :
   fRecTracksArray(0),  //clones array with filtered particles
   fDebugStreamer(0),
   fStreamLevel(0),
-  fDebugLevel(0)
+  fDebugLevel(0),
+  fDebugOutputPath()
 {
   //
   // Default constructor (should not be used)
@@ -59,6 +65,7 @@ AliGenInfoTask::AliGenInfoTask(const AliGenInfoTask& /*info*/) :
   AliAnalysisTask(), 
   fMCinfo(0),     //! MC event handler
   fESD(0),
+  fESDfriend(0),
   fCompList(0),
   fGenTracksArray(0),  //clones array with filtered particles
   fGenKinkArray(0),    //clones array with filtered Kinks
@@ -67,7 +74,8 @@ AliGenInfoTask::AliGenInfoTask(const AliGenInfoTask& /*info*/) :
   //
   fDebugStreamer(0),
   fStreamLevel(0),
-  fDebugLevel()
+  fDebugLevel(),
+  fDebugOutputPath()
 {
   //
   // Default constructor 
@@ -81,6 +89,7 @@ AliGenInfoTask::AliGenInfoTask(const char *name) :
   AliAnalysisTask(name, "AliGenInfoTask"), 
   fMCinfo(0),     //! MC event handler
   fESD(0),
+  fESDfriend(0),
   fCompList(0),
   fGenTracksArray(0),  //clones array with filtered particles
   fGenKinkArray(0),    //clones array with filtered Kinks
@@ -88,7 +97,8 @@ AliGenInfoTask::AliGenInfoTask(const char *name) :
   fRecTracksArray(0),  //clones array with filtered particles
   fDebugStreamer(0),
   fStreamLevel(0),
-  fDebugLevel(0)
+  fDebugLevel(0),
+  fDebugOutputPath()
 {
   //
   // Normal constructor
@@ -134,6 +144,8 @@ void AliGenInfoTask::ConnectInputData(Option_t *)
     }
     else {
       fESD = esdH->GetEvent();
+      fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
+      fESD->SetESDfriend(fESDfriend);
       //Printf("*** CONNECTED NEW EVENT ****");
     }  
   }
@@ -219,6 +231,8 @@ void AliGenInfoTask::Terminate(Option_t *) {
 
 
 
+
+
 TTreeSRedirector *AliGenInfoTask::GetDebugStreamer(){
   //
   // Get Debug streamer
@@ -314,6 +328,12 @@ void AliGenInfoTask::ProcessESDInfo(){
   //
   //
   static AliTPCParamSR param;
+  fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
+  if (!fESDfriend) {
+    //Printf("ERROR: fESDfriend not available");
+    return;
+  }
+  fESD->SetESDfriend(fESDfriend);
   //
   //
   if (!fESD) return;
@@ -328,31 +348,6 @@ void AliGenInfoTask::ProcessESDInfo(){
     recInfo->Update(mcinfo,&param,kTRUE);
   }
   //
-  //
-  //
-  Int_t ntracksMC = fMCinfo->GetNumberOfTracks();
-  for (Int_t imc=0; imc<ntracksMC; imc++){
-    AliMCInfo * mcinfo = GetTrack(imc,kFALSE);
-    if (!mcinfo) continue;
-    AliESDRecInfo *recInfo= GetRecTrack(imc,kFALSE);
-    if (recInfo) continue;
-    if (mcinfo->GetNTPCRef()<2) continue;
-    //
-    //
-    for (Int_t itrack=0; itrack<ntracks; itrack++){
-      AliESDtrack *track = fESD->GetTrack(itrack);
-      Int_t label = TMath::Abs(track->GetLabel());
-      if (label!=mcinfo->GetLabel()) continue;
-      
-      AliMCInfo * mcinfo2 = GetTrack(label,kFALSE);
-      if (!mcinfo2) continue;
-      AliESDRecInfo *recInfo= GetRecTrack(label,kTRUE);
-      recInfo->AddESDtrack(track,mcinfo2);
-      recInfo->Update(mcinfo2,&param,kTRUE);
-    }
-  }
-
-  
 
 
 } 
@@ -384,6 +379,10 @@ void AliGenInfoTask::DumpInfo(){
   //
   //
   //
+  TParticle * particle= new TParticle;
+  TClonesArray * trefs = new TClonesArray("AliTrackReference");
+  //
+
   static AliESDRecInfo dummy;
   Int_t npart = fMCinfo->GetNumberOfTracks();
   for (Int_t ipart=0;ipart<npart;ipart++){
@@ -392,10 +391,21 @@ void AliGenInfoTask::DumpInfo(){
     AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE);
     if (!recInfo) recInfo=&dummy; 
     TTreeSRedirector *pcstream = GetDebugStreamer();
+    
+    fMCinfo->GetParticleAndTR(ipart, particle, trefs);
+    Int_t counter=0;
+    Float_t length=0;
+    if (trefs!=0 && particle!=0){
+      length = GetTPCTrackLength(*trefs, particle , fESD->GetMagneticField(), counter, 3.0);
+    }
+
     if (pcstream){
-      (*pcstream)<<"RC"<<
+      (*pcstream)<<"RC"<<      
        "MC.="<<mcinfo<<
        "RC.="<<recInfo<<
+       "length="<<length<<
+       "counter="<<counter<<
+       //      "tr.="<<tpctrack<<
        "\n";
     }    
   }
@@ -423,3 +433,109 @@ Bool_t AliGenInfoTask::AcceptParticle(TParticle *part){
   return kTRUE;
 }
 
+
+
+void AliGenInfoTask::FinishTaskOutput()
+{
+  //
+  // According description in AliAnalisysTask this method is call
+  // on the slaves before sending data
+  //
+  Terminate("slave");
+  RegisterDebugOutput(fDebugOutputPath.Data());
+
+}
+
+
+
+void AliGenInfoTask::RegisterDebugOutput(const char */*path*/){
+  //
+  // 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+="/";;
+  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());
+}
+
+
+Float_t  AliGenInfoTask::GetTPCTrackLength(const TClonesArray& trackRefs, TParticle*part, Float_t bz,  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= 0.01; 
+  Float_t length =0;
+  Int_t nrefs = trackRefs.GetEntriesFast();
+  AliExternalTrackParam param;
+  Double_t cv[21];
+  for (Int_t i=0; i<21;i++) cv[i]=0;
+  counter=0;
+  //
+  //
+  AliTrackReference *ref0 = (AliTrackReference*)trackRefs.At(0);
+  Float_t direction = 0;
+  //
+  for (Int_t iref=1; iref<nrefs;iref++){
+    AliTrackReference *ref = (AliTrackReference*)trackRefs.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(part->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;
+}