[u/mrichter/AliRoot.git] / PWGPP / TPC / AliGenInfoTask.cxx

//
// This class is the task for connecting together 
// MC information and the RC information 
//
// The task is a wrapper over two components
// AliGenInfoMaker
// AliESDRecInfoMaker.h

// ROOT includes
#include <TChain.h>
#include <TMath.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 "AliGenInfoTask.h"
#include "AliGenInfoMaker.h"
#include "AliHelix.h"

//
#include "AliMCInfo.h"
#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>

using namespace std;

ClassImp(AliGenInfoTask)

//________________________________________________________________________
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
  fGenV0Array(0),      //clones array with filtered V0s
  fRecTracksArray(0),  //clones array with filtered particles
  fDebugStreamer(0),
  fStreamLevel(0),
  fDebugLevel(0),
  fDebugOutputPath()
{
  //
  // Default constructor (should not be used)
  //
}

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
  fGenV0Array(0),      //clones array with filtered V0s
  fRecTracksArray(0),  //clones array with filtered particles
  //
  fDebugStreamer(0),
  fStreamLevel(0),
  fDebugLevel(),
  fDebugOutputPath()
{
  //
  // Default constructor 
  //
}


//________________________________________________________________________
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
  fGenV0Array(0),      //clones array with filtered V0s
  fRecTracksArray(0),  //clones array with filtered particles
  fDebugStreamer(0),
  fStreamLevel(0),
  fDebugLevel(0),
  fDebugOutputPath()
{
  //
  // Normal constructor
  //
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TList
  DefineOutput(0, TObjArray::Class());
  //
  //
  fCompList = new TObjArray;
}

AliGenInfoTask::~AliGenInfoTask(){
  //
  //
  //
  if (fDebugLevel>0)  printf("AliGenInfoTask::~AliGenInfoTask\n");
  if (fDebugStreamer) delete fDebugStreamer;
  fDebugStreamer=0;
  if(fCompList)   delete fCompList;  
  fCompList =0; 
}


//________________________________________________________________________
void AliGenInfoTask::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 {
      fESD = esdH->GetEvent();
      fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
      fESD->SetESDfriend(fESDfriend);
      //Printf("*** CONNECTED NEW EVENT ****");
    }  
  }
  AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());  
  mcinfo->SetReadTR(kTRUE);
  
  fMCinfo = mcinfo->MCEvent();


}


//_____________________________________________________________________________
Bool_t AliGenInfoTask::AddComparisonObject(AliComparisonObject *pObj) 
{
  // add comparison object to the list
  if(pObj == 0) {
      Printf("ERROR: Could not add comparison object");
	  return kFALSE;
  }
  // add object to the list
  fCompList->AddLast(pObj);       
  return kTRUE;
}


//________________________________________________________________________
void AliGenInfoTask::CreateOutputObjects() 
{
  //
  // Connect the output objects
  //
  if(fDebugLevel>3)
    cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
  
}


//________________________________________________________________________
void AliGenInfoTask::Exec(Option_t *) {
  //
  // Execute analysis for current event 
  //

  if(fDebugLevel>3)
    cout << "AliGenInfoTask::Exec()" << endl;
    

  // If MC has been connected   
  if (fGenTracksArray) fGenTracksArray->Delete();
  if (fRecTracksArray) fRecTracksArray->Delete();

  if (!fMCinfo){
    cout << "Not MC info\n" << endl;
  }else{
    //mcinfo->Print();
    ProcessMCInfo();
    ProcessESDInfo();
    DumpInfo();
    ProcessComparison();
  }
  //
}      


//________________________________________________________________________
void AliGenInfoTask::Terminate(Option_t *) {
    //
    // Terminate loop
    //
  if(fDebugLevel>3)
    printf("AliGenInfoTask: Terminate() \n");  
  //
  if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
  if (fDebugStreamer) delete fDebugStreamer;
  fDebugStreamer = 0;
  return;
}


TTreeSRedirector *AliGenInfoTask::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;
}


AliMCInfo*  AliGenInfoTask::GetTrack(Int_t index, Bool_t force){
  //
  // Get the MC info for given track
  //
  if (!fGenTracksArray) fGenTracksArray = new TClonesArray("AliMCInfo",1000);
  if (index>fGenTracksArray->GetEntriesFast()) fGenTracksArray->Expand(index*2+10);
  AliMCInfo * info = (AliMCInfo*)fGenTracksArray->At(index);
  if (!force) return info;
  if (!info){
    info = new ((*fGenTracksArray)[index]) AliMCInfo;
  }
  return info;
}

AliESDRecInfo*  AliGenInfoTask::GetRecTrack(Int_t index, Bool_t force){
  //
  // Get the MC info for given track
  //
  if (!fRecTracksArray) fRecTracksArray = new TClonesArray("AliESDRecInfo",1000);
  if (index>fRecTracksArray->GetEntriesFast()) fRecTracksArray->Expand(index*2+10);
  AliESDRecInfo * info = (AliESDRecInfo*)fRecTracksArray->At(index);
  if (!force) return info;
  if (!info){
    info = new ((*fRecTracksArray)[index]) AliESDRecInfo;
  }
  return info;
}


void  AliGenInfoTask::ProcessMCInfo(){
  //
  // Dump information from MC to the array
  //
  //
  TParticle * particle= new TParticle;
  TClonesArray * trefs = new TClonesArray("AliTrackReference");
  //
  //
  // 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();
  }

  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 (!AcceptParticle(particle)) continue;
    //if (trefs->GetEntries()<1) continue;
    AliMCInfo * mcinfo = GetTrack(ipart,kTRUE);
    mcinfo->Update(particle,trefs,vertex,ipart);
    //
    TTreeSRedirector *pcstream = GetDebugStreamer();
    if (pcstream){
      (*pcstream)<<"MC"<<
	"p.="<<particle<<
	"MC.="<<mcinfo<<
	"\n";
    }
  }
}

void AliGenInfoTask::ProcessESDInfo(){
  //
  //
  //
  if (!fESD) return;
  static AliTPCParamSR param;
  fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
  if (!fESDfriend) {
    //Printf("ERROR: fESDfriend not available");
    return;
  }
  fESD->SetESDfriend(fESDfriend);
  //
  //
  Int_t ntracks = fESD->GetNumberOfTracks();
  for (Int_t itrack=0; itrack<ntracks; itrack++){
    AliESDtrack *track = fESD->GetTrack(itrack);
    Int_t label = TMath::Abs(track->GetLabel());
    AliMCInfo * mcinfo = GetTrack(label,kFALSE);
    if (!mcinfo) continue;
    AliESDRecInfo *recInfo= GetRecTrack(label,kTRUE);
    recInfo->AddESDtrack(track,mcinfo);
    recInfo->Update(mcinfo,&param,kTRUE);
  }
  //


} 


void AliGenInfoTask::ProcessComparison(){
  //
  //
  //
  static AliESDRecInfo dummy;
  Int_t npart = fMCinfo->GetNumberOfTracks();
  for (Int_t ipart=0;ipart<npart;ipart++){
    AliMCInfo * mcinfo = GetTrack(ipart,kFALSE);
    if (!mcinfo) continue;
    AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE);
    if (!recInfo) recInfo=&dummy; 
    //
    for (Int_t icomp = 0; icomp<fCompList->GetEntries(); icomp++){
      AliComparisonObject *pObj= (AliComparisonObject *)fCompList->At(icomp);
      if (pObj){
	pObj->Exec(mcinfo,recInfo);
      }
    }    
  }
  PostData(0, fCompList);
}

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++){
    AliMCInfo * mcinfo = GetTrack(ipart,kFALSE);
    if (!mcinfo) continue;
    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"<<	
	"MC.="<<mcinfo<<
	"RC.="<<recInfo<<
	"length="<<length<<
	"counter="<<counter<<
	//	"tr.="<<tpctrack<<
	"\n";
    }    
  }
}


Bool_t AliGenInfoTask::AcceptParticle(TParticle *part){
  //
  /*
    MC cuts
    TCut cutPt("p.Pt()>0.1");
    TCut cutZ("abs(p.Vz())<250");
    TCut cutR("abs(p.R())<250");
    //
  */
  //
  //
  if (part->Pt()<0.1) return kFALSE;
  if (TMath::Abs(part->Vz())>250) return kFALSE;
  if (part->R()>360)  return kFALSE;
  if (part->GetPDG()){
    if (part->GetPDG()->Charge()==0) return kFALSE;
  }
  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;
}
Commit	Line	Data
7cc34f08	1	//
	2	// This class is the task for connecting together
	3	// MC information and the RC information
	4	//
	5	// The task is a wrapper over two components
	6	// AliGenInfoMaker
	7	// AliESDRecInfoMaker.h
	8
	9	// ROOT includes
	10	#include <TChain.h>
	11	#include <TMath.h>
	12
	13	// ALIROOT includes
	14	#include <TTreeStream.h>
	15	#include <AliAnalysisManager.h>
	16	#include <AliESDInputHandler.h>
	17	#include "AliStack.h"
	18	#include "AliMCEvent.h"
	19	#include "AliMCEventHandler.h"
	20
	21	#include <AliESD.h>
	22	#include "AliGenInfoTask.h"
	23	#include "AliGenInfoMaker.h"
	24	#include "AliHelix.h"
	25
	26	//
	27	#include "AliMCInfo.h"
	28	#include "AliComparisonObject.h"
	29	#include "AliESDRecInfo.h"
	30	#include "AliTPCParamSR.h"
	31	#include "TSystem.h"
	32	#include "TTimeStamp.h"
	33	#include "TFile.h"
	34	#include "AliTPCseed.h"
	35
	36	// STL includes
	37	#include <iostream>
	38
	39	using namespace std;
	40
	41	ClassImp(AliGenInfoTask)
	42
	43	//________________________________________________________________________
	44	AliGenInfoTask::AliGenInfoTask() :
	45	AliAnalysisTask(),
	46	fMCinfo(0), //! MC event handler
	47	fESD(0),
	48	fESDfriend(0),
	49	fCompList(0), //array of comparison objects
	50	fGenTracksArray(0), //clones array with filtered particles
	51	fGenKinkArray(0), //clones array with filtered Kinks
	52	fGenV0Array(0), //clones array with filtered V0s
	53	fRecTracksArray(0), //clones array with filtered particles
	54	fDebugStreamer(0),
	55	fStreamLevel(0),
	56	fDebugLevel(0),
	57	fDebugOutputPath()
	58	{
	59	//
	60	// Default constructor (should not be used)
	61	//
	62	}
	63
	64	AliGenInfoTask::AliGenInfoTask(const AliGenInfoTask& /info/) :
65	AliAnalysisTask(),
66	fMCinfo(0), //! MC event handler
67	fESD(0),
68	fESDfriend(0),
69	fCompList(0),
70	fGenTracksArray(0), //clones array with filtered particles
71	fGenKinkArray(0), //clones array with filtered Kinks
72	fGenV0Array(0), //clones array with filtered V0s
73	fRecTracksArray(0), //clones array with filtered particles
74	//
75	fDebugStreamer(0),
76	fStreamLevel(0),
77	fDebugLevel(),
78	fDebugOutputPath()
79	{
80	//
81	// Default constructor
82	//
83	}
84
85
86
87	//________________________________________________________________________
88	AliGenInfoTask::AliGenInfoTask(const char *name) :
89	AliAnalysisTask(name, "AliGenInfoTask"),
90	fMCinfo(0), //! MC event handler
91	fESD(0),
92	fESDfriend(0),
93	fCompList(0),
94	fGenTracksArray(0), //clones array with filtered particles
95	fGenKinkArray(0), //clones array with filtered Kinks
96	fGenV0Array(0), //clones array with filtered V0s
97	fRecTracksArray(0), //clones array with filtered particles
98	fDebugStreamer(0),
99	fStreamLevel(0),
100	fDebugLevel(0),
101	fDebugOutputPath()
102	{
103	//
104	// Normal constructor
105	//
106	// Input slot #0 works with a TChain
107	DefineInput(0, TChain::Class());
108	// Output slot #0 writes into a TList
109	DefineOutput(0, TObjArray::Class());
110	//
111	//
112	fCompList = new TObjArray;
113	}
114
115	AliGenInfoTask::~AliGenInfoTask(){
116	//
117	//
118	//
119	if (fDebugLevel>0) printf("AliGenInfoTask::~AliGenInfoTask\n");
120	if (fDebugStreamer) delete fDebugStreamer;
121	fDebugStreamer=0;
122	if(fCompList) delete fCompList;
123	fCompList =0;
124	}
125
126
127	//________________________________________________________________________
128	void AliGenInfoTask::ConnectInputData(Option_t *)
129	{
130	//
131	// Connect the input data
132	//
133	if(fDebugLevel>3)
134	cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
135
136	TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
137	if (!tree) {
138	//Printf("ERROR: Could not read chain from input slot 0");
139	}
140	else {
141	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
142	if (!esdH) {
143	//Printf("ERROR: Could not get ESDInputHandler");
144	}
145	else {
146	fESD = esdH->GetEvent();
147	fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
148	fESD->SetESDfriend(fESDfriend);
149	//Printf("* CONNECTED NEW EVENT **");
150	}
151	}
152	AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
153	mcinfo->SetReadTR(kTRUE);
154
155	fMCinfo = mcinfo->MCEvent();
156
157
158	}
159
160
161	//_____________________________________________________________________________
162	Bool_t AliGenInfoTask::AddComparisonObject(AliComparisonObject *pObj)
163	{
164	// add comparison object to the list
165	if(pObj == 0) {
166	Printf("ERROR: Could not add comparison object");
167	return kFALSE;
168	}
169	// add object to the list
170	fCompList->AddLast(pObj);
171	return kTRUE;
172	}
173
174
175
176
177	//________________________________________________________________________
178	void AliGenInfoTask::CreateOutputObjects()
179	{
180	//
181	// Connect the output objects
182	//
183	if(fDebugLevel>3)
184	cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
185
186	}
187
188
189	//________________________________________________________________________
190	void AliGenInfoTask::Exec(Option_t *) {
191	//
192	// Execute analysis for current event
193	//
194
195	if(fDebugLevel>3)
196	cout << "AliGenInfoTask::Exec()" << endl;
197
198
199	// If MC has been connected
200	if (fGenTracksArray) fGenTracksArray->Delete();
201	if (fRecTracksArray) fRecTracksArray->Delete();
202
203	if (!fMCinfo){
204	cout << "Not MC info\n" << endl;
205	}else{
206	//mcinfo->Print();
207	ProcessMCInfo();
208	ProcessESDInfo();
209	DumpInfo();
210	ProcessComparison();
211	}
212	//
213	}
214
215
216
217
218	//________________________________________________________________________
219	void AliGenInfoTask::Terminate(Option_t *) {
220	//
221	// Terminate loop
222	//
223	if(fDebugLevel>3)
224	printf("AliGenInfoTask: Terminate() \n");
225	//
226	if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
227	if (fDebugStreamer) delete fDebugStreamer;
228	fDebugStreamer = 0;
229	return;
230	}
231
232
233
234
235
236	TTreeSRedirector *AliGenInfoTask::GetDebugStreamer(){
237	//
238	// Get Debug streamer
239	// In case debug streamer not yet initialized and StreamLevel>0 create new one
240	//
241	if (fStreamLevel==0) return 0;
242	if (fDebugStreamer) return fDebugStreamer;
243	TString dsName;
244	dsName=GetName();
245	dsName+="Debug.root";
246	dsName.ReplaceAll(" ","");
247	fDebugStreamer = new TTreeSRedirector(dsName.Data());
248	return fDebugStreamer;
249	}
250
251
252
253	AliMCInfo* AliGenInfoTask::GetTrack(Int_t index, Bool_t force){
254	//
255	// Get the MC info for given track
256	//
257	if (!fGenTracksArray) fGenTracksArray = new TClonesArray("AliMCInfo",1000);
258	if (index>fGenTracksArray->GetEntriesFast()) fGenTracksArray->Expand(index*2+10);
259	AliMCInfo * info = (AliMCInfo*)fGenTracksArray->At(index);
260	if (!force) return info;
261	if (!info){
262	info = new ((*fGenTracksArray)[index]) AliMCInfo;
263	}
264	return info;
265	}
266
267	AliESDRecInfo* AliGenInfoTask::GetRecTrack(Int_t index, Bool_t force){
268	//
269	// Get the MC info for given track
270	//
271	if (!fRecTracksArray) fRecTracksArray = new TClonesArray("AliESDRecInfo",1000);
272	if (index>fRecTracksArray->GetEntriesFast()) fRecTracksArray->Expand(index*2+10);
273	AliESDRecInfo * info = (AliESDRecInfo*)fRecTracksArray->At(index);
274	if (!force) return info;
275	if (!info){
276	info = new ((*fRecTracksArray)[index]) AliESDRecInfo;
277	}
278	return info;
279	}
280
281
282
283
284	void AliGenInfoTask::ProcessMCInfo(){
285	//
286	// Dump information from MC to the array
287	//
288	//
289	TParticle * particle= new TParticle;
290	TClonesArray * trefs = new TClonesArray("AliTrackReference");
291	//
292	//
293	// Process tracks
294	//
295	Int_t npart = fMCinfo->GetNumberOfTracks();
296	if (npart==0) return;
297	Double_t vertex[4]={0,0,0,0};
298	fMCinfo->GetParticleAndTR(0, particle, trefs);
299	if (particle){
300	vertex[0]=particle->Vx();
301	vertex[1]=particle->Vy();
302	vertex[2]=particle->Vz();
303	vertex[3]=particle->R();
304	}
305
306	for (Int_t ipart=0;ipart<npart;ipart++){
307	Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
308	if (status<0) continue;
309	if (!particle) continue;
310	if (!trefs) continue;
311	if (!AcceptParticle(particle)) continue;
312	//if (trefs->GetEntries()<1) continue;
313	AliMCInfo * mcinfo = GetTrack(ipart,kTRUE);
314	mcinfo->Update(particle,trefs,vertex,ipart);
315	//
316	TTreeSRedirector *pcstream = GetDebugStreamer();
317	if (pcstream){
318	(*pcstream)<<"MC"<<
319	"p.="<<particle<<
320	"MC.="<<mcinfo<<
321	"\n";
322	}
323	}
324	}
325
326	void AliGenInfoTask::ProcessESDInfo(){
327	//
328	//
329	//
4a9220d4	330	if (!fESD) return;
7cc34f08	331	static AliTPCParamSR param;
	332	fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
	333	if (!fESDfriend) {
	334	//Printf("ERROR: fESDfriend not available");
	335	return;
	336	}
	337	fESD->SetESDfriend(fESDfriend);
	338	//
	339	//
7cc34f08	340	Int_t ntracks = fESD->GetNumberOfTracks();
	341	for (Int_t itrack=0; itrack<ntracks; itrack++){
	342	AliESDtrack *track = fESD->GetTrack(itrack);
	343	Int_t label = TMath::Abs(track->GetLabel());
	344	AliMCInfo * mcinfo = GetTrack(label,kFALSE);
	345	if (!mcinfo) continue;
	346	AliESDRecInfo *recInfo= GetRecTrack(label,kTRUE);
	347	recInfo->AddESDtrack(track,mcinfo);
	348	recInfo->Update(mcinfo,&param,kTRUE);
	349	}
	350	//
	351
	352
	353	}
	354
	355
	356	void AliGenInfoTask::ProcessComparison(){
	357	//
	358	//
	359	//
	360	static AliESDRecInfo dummy;
	361	Int_t npart = fMCinfo->GetNumberOfTracks();
	362	for (Int_t ipart=0;ipart<npart;ipart++){
	363	AliMCInfo * mcinfo = GetTrack(ipart,kFALSE);
	364	if (!mcinfo) continue;
	365	AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE);
	366	if (!recInfo) recInfo=&dummy;
	367	//
	368	for (Int_t icomp = 0; icomp<fCompList->GetEntries(); icomp++){
	369	AliComparisonObject pObj= (AliComparisonObject )fCompList->At(icomp);
	370	if (pObj){
	371	pObj->Exec(mcinfo,recInfo);
	372	}
	373	}
	374	}
	375	PostData(0, fCompList);
	376	}
	377
	378	void AliGenInfoTask::DumpInfo(){
	379	//
	380	//
	381	//
	382	TParticle * particle= new TParticle;
	383	TClonesArray * trefs = new TClonesArray("AliTrackReference");
	384	//
	385
	386	static AliESDRecInfo dummy;
	387	Int_t npart = fMCinfo->GetNumberOfTracks();
	388	for (Int_t ipart=0;ipart<npart;ipart++){
	389	AliMCInfo * mcinfo = GetTrack(ipart,kFALSE);
	390	if (!mcinfo) continue;
	391	AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE);
	392	if (!recInfo) recInfo=&dummy;
	393	TTreeSRedirector *pcstream = GetDebugStreamer();
	394
	395	fMCinfo->GetParticleAndTR(ipart, particle, trefs);
	396	Int_t counter=0;
	397	Float_t length=0;
	398	if (trefs!=0 && particle!=0){
	399	length = GetTPCTrackLength(*trefs, particle , fESD->GetMagneticField(), counter, 3.0);
	400	}
	401
	402	if (pcstream){
	403	(*pcstream)<<"RC"<<
404	"MC.="<<mcinfo<<
405	"RC.="<<recInfo<<
406	"length="<<length<<
407	"counter="<<counter<<
408	// "tr.="<<tpctrack<<
409	"\n";
410	}
411	}
412	}
413
414
415
416	Bool_t AliGenInfoTask::AcceptParticle(TParticle *part){
417	//
418	/*
419	MC cuts
420	TCut cutPt("p.Pt()>0.1");
421	TCut cutZ("abs(p.Vz())<250");
422	TCut cutR("abs(p.R())<250");
423	//
424	*/
425	//
426	//
427	if (part->Pt()<0.1) return kFALSE;
428	if (TMath::Abs(part->Vz())>250) return kFALSE;
429	if (part->R()>360) return kFALSE;
430	if (part->GetPDG()){
431	if (part->GetPDG()->Charge()==0) return kFALSE;
432	}
433	return kTRUE;
434	}
435
436
437
438	void AliGenInfoTask::FinishTaskOutput()
439	{
440	//
441	// According description in AliAnalisysTask this method is call
442	// on the slaves before sending data
443	//
444	Terminate("slave");
445	RegisterDebugOutput(fDebugOutputPath.Data());
446
447	}
448
449
450
451	void AliGenInfoTask::RegisterDebugOutput(const char /path*/){
452	//
453	// store - copy debug output to the destination position
454	// currently ONLY for local copy
455	TString dsName;
456	dsName=GetName();
457	dsName+="Debug.root";
458	dsName.ReplaceAll(" ","");
459	TString dsName2=fDebugOutputPath.Data();
460	gSystem->MakeDirectory(dsName2.Data());
461	dsName2+="/";;
462	dsName2+=gSystem->HostName();
463	gSystem->MakeDirectory(dsName2.Data());
464	dsName2+="/";
465	dsName2+=gSystem->BaseName(gSystem->pwd());
466	dsName2+="/";
467	gSystem->MakeDirectory(dsName2.Data());
468	dsName2+=dsName;
469	AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
470	printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
471	TFile::Cp(dsName.Data(),dsName2.Data());
472	}
473
474
475	Float_t AliGenInfoTask::GetTPCTrackLength(const TClonesArray& trackRefs, TParticle*part, Float_t bz, Int_t &counter, Float_t deadWidth){
476	//
477	// return track length in geometrically active volume of TPC.
478	// z nad rphi acceptance is included
479	// doesn't take into account dead channel and ExB
480	// Intput:
481	// trackRefs
482	// bz - magnetic field
483	// deadWidth - dead zone in r-phi
484	// Additional output:
485	// counter - number of circles
486	const Float_t kRMin = 90;
487	const Float_t kRMax = 245;
488	const Float_t kZMax = 250;
489	const Float_t kMinPt= 0.01;
490	Float_t length =0;
491	Int_t nrefs = trackRefs.GetEntriesFast();
492	AliExternalTrackParam param;
493	Double_t cv[21];
494	for (Int_t i=0; i<21;i++) cv[i]=0;
495	counter=0;
496	//
497	//
498	AliTrackReference ref0 = (AliTrackReference)trackRefs.At(0);
499	Float_t direction = 0;
500	//
501	for (Int_t iref=1; iref<nrefs;iref++){
502	AliTrackReference ref = (AliTrackReference)trackRefs.At(iref);
503	if (!ref) continue;
504	if (!ref0 \|\| ref0->DetectorId()!= AliTrackReference::kTPC){
505	ref0 = ref;
506	direction = ((ref0->X()ref0->Px()+ref0->Y()ref0->Py())>0)? 1.:-1.;
507	continue;
508	}
509	Float_t newdirection = ((ref->X()ref->Px()+ref->Y()ref->Py())>0)? 1.:-1.;
510	if (newdirection*direction<0) {
511	counter++; //circle counter
512	direction = newdirection;
513	continue;
514	}
515	if (counter>0) continue;
516	if (ref0->Pt()<kMinPt) break;
517	Float_t radius0 = TMath::Max(TMath::Min(ref0->R(),kRMax),kRMin);;
518	Float_t radius1 = TMath::Max(TMath::Min(ref->R(),kRMax),kRMin);
519	Double_t xyz[3]={ref0->X(),ref0->Y(),ref0->Z()};
520	Double_t pxyz[3]={ref0->Px(),ref0->Py(),ref0->Pz()};
521	Double_t alpha;
522	param.Set(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
523
524	for (Float_t radius = radius0; radius<radius1; radius+=1){
525	param.GetXYZAt(radius, bz, xyz);
526	if (TMath::Abs(xyz[2])>kZMax) continue;
527	Float_t gradius = TMath::Sqrt(xyz[1]xyz[1]+xyz[0]xyz[0]);
528	if (gradius>kRMax) continue;
529	alpha = TMath::ATan2(xyz[1],xyz[0]);
530	if (alpha<0) alpha+=TMath::TwoPi();
531	//
532	Int_t sector = Int_t(9*alpha/TMath::Pi());
533	Float_t lalpha = alpha-((sector+0.5)*TMath::Pi()/9.);
534	Float_t dedge = (TMath::Tan(TMath::Pi()/18.)-TMath::Abs(TMath::Tan(lalpha)))*gradius;
535	if (dedge>deadWidth) length++;
536	}
537	if (ref->DetectorId()!= AliTrackReference::kTPC) break;
538	ref0 = ref;
539	}
540	return length;
541	}