[u/mrichter/AliRoot.git] / HBTAN / hbtanalysis.C

  #include "AliHBTAnalysis.h"
//  #include "alles.h"
  #include "AliHBTReader.h"
  #include "AliHBTReaderKineTree.h"
  #include "AliHBTReaderITSv2.h"
  #include "AliHBTReaderITSv1.h"
  #include "AliHBTReaderTPC.h"
  #include "AliHBTReaderInternal.h"
  #include "AliHBTParticleCut.h"
  #include "AliHBTEvent.h"
  #include "AliHBTPairCut.h"
  #include "AliHBTQResolutionFctns.h"
  #include "AliHBTQDistributionFctns.h"
  #include "AliHBTTwoTrackEffFctn.h"
  #include "AliHBTCorrelFctn.h"
  #include "TSystem.h"
  #include "TObjString.h"
  #include "TString.h"
  #include "TPDGCode.h"
  #include "AliHBTMonDistributionFctns.h"
  #include "AliHBTMonResolutionFctns.h"
  
void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles", 
                Int_t first = -1,Int_t last = -1, 
                char *outfile = "hbtanalysis.root")
 {
   
  AliHBTTrackPoints::SetDebug(0);
  AliHBTParticle::SetDebug(0);
  AliLoader::SetDebug(0);
  
  AliHBTParticleCut c1;
  AliHBTParticleCut c2 = c1;

//HBT Anlysis Macro
//Anlyzes TPC recontructed tracks and simulated particles that corresponds to them

//datatype defines type of data to be read
//  Kine  - analyzes Kine Tree: simulated particles
//  TPC   - analyzes TPC   tracking + particles corresponding to these tracks
//  ITSv1 - analyzes ITSv1 ----------------------//--------------------------
//  ITSv2 - analyzes ITSv2 ----------------------//--------------------------

//processopt defines option passed to AliHBTAnlysis::Process method
// default: TracksAndParticles - process both recontructed tracks and sim. particles corresponding to them 
//          Tracks - process only recontructed tracks
//          Particles - process only simulated particles

//Reads data from diroctories from first to last(including)
// For examples if first=3 and last=5 it reads from
//  ./3/
//  ./4/
//  ./5/
//if first or last is negative (or both), it reads from current directory
//
//these names I use when analysis is done directly from CASTOR files via RFIO
//  const char* basedir="rfio:/castor/cern.ch/user/s/skowron/";
//  const char* serie="standard";
//  const char* field = "0.4";

  const char* basedir=".";
  const char* serie="";
  const char* field = "";
 
  Bool_t threeDcuts = kFALSE;

  // Dynamically link some shared libs                    
//  gROOT->LoadMacro("loadlibs.C");
//  loadlibs();
  
  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libHBTAN");
  

  /***********************************************************/
  //Create analysis and reader
  AliHBTAnalysis * analysis = new AliHBTAnalysis();
  analysis->SetCutsOnTracks();
  
  AliHBTReader* reader;
  Int_t kine = strcmp(datatype,"Kine");
  Int_t ESD = strcmp(datatype,"ESD");
  Int_t TPC = strcmp(datatype,"TPC");
  Int_t ITSv1 = strcmp(datatype,"ITSv1");
  Int_t ITSv2 = strcmp(datatype,"ITSv2");
  Int_t intern = strcmp(datatype,"Intern");
  
  if(!kine)
   {
    reader = new AliHBTReaderKineTree();
    processopt="Particles"; //this reader by definition reads only simulated particles
   }
  else if(!ESD)
   {
    AliHBTReaderESD* esdreader = new AliHBTReaderESD();
    esdreader->ReadParticles(kTRUE);
    esdreader->SetNumberOfTrackPoints(5,30.);
    esdreader->SetClusterMap();
    reader = esdreader;
   }
  else if(!TPC)
   {
    reader = new AliHBTReaderTPC();
    //((AliHBTReaderTPC*)reader)->SetNumberOfTrackPoints(5,30.);
    ((AliHBTReaderTPC*)reader)->SetClusterMap();
   }
  else if(!ITSv1)
   {
    reader = new AliHBTReaderITSv1();
   }
  else if(!ITSv2)
   {
    reader = new AliHBTReaderITSv2();
   }
  else if(!intern)
   {
    reader = new AliHBTReaderInternal("data.root");
   }
  else
   {
    cerr<<"Option "<<datatype<<"  not recognized. Exiting"<<endl;
    return;
   }

  TObjArray* dirs=0;
  if ( (first >= 0) && (last>=0) && ( (last-first)>=0 ) )
   {//read from many dirs dirs
     char buff[50];
     dirs = new TObjArray(last-first+1);
     for (Int_t i = first; i<=last; i++)
      {
        sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i);
        TObjString *odir= new TObjString(buff);
        dirs->Add(odir);
      }
    }
   reader->SetDirs(dirs);
   
  /***********************************************************/    
  /*****   R E A D E R ***************************************/    
  /***********************************************************/    
  
  //we want have only low pt pi+ so set a cut to reader
  AliHBTParticleCut* readerpartcut= new AliHBTParticleCut();
  readerpartcut->SetPtRange(0.0,10.0);
  readerpartcut->SetPID(kPiPlus);
  reader->AddParticleCut(readerpartcut);//read this particle type with this cut
  
//  readerpartcut->SetPtRange(0.0,10.0);
//  readerpartcut->SetPID(kPiMinus);
//  reader->AddParticleCut(readerpartcut);//read this particle type with this cut
  
  analysis->SetReader(reader);

  AliHBTPairCut *paircut = new AliHBTPairCut();
  Float_t qinvmin = 0.0;
  Float_t qinvmax = 0.05;//50MeV
  paircut->SetQInvRange(qinvmin,qinvmax);  
  
 // paircut->SetAvSeparationRange(10.); //AntiMerging
  paircut->SetClusterOverlapRange(-.5,1.0);//AntiSplitting
  
//  AliHBTParticleCut* partcut= new AliHBTParticleCut();
//  partcut->SetPID(kPiPlus);
//  paircut->SetFirstPartCut(partcut);
//  partcut->SetPID(kPiMinus);
//  paircut->SetSecondPartCut(partcut);
  
  analysis->SetGlobalPairCut(paircut);
  
  /***********************************************************/    
  /*****   W E I G H T S        ******************************/    
  /***********************************************************/    

  AliHBTLLWeights::Instance()->SetParticlesTypes(kPiPlus,kPiPlus);
  AliHBTLLWeights::Instance()->SetColoumb(kFALSE);
  AliHBTLLWeights::Instance()->SetStrongInterSwitch(kFALSE);

  AliHBTLLWeights::Instance()->SetRandomPosition(kFALSE);
//AliHBTLLWeights::Instance()->SetR1dw(8.0);

  AliHBTLLWeights::Instance()->Init();
  AliHBTLLWeights::Instance()->Set();

  //example function
  AliHBTWeightTheorQInvFctn  *wqinvcfP = new AliHBTWeightTheorQInvFctn(100,qinvmax);
  wqinvcfP->SetNumberOfBinsToScale(30);
  wqinvcfP->Rename("wqinvcfP","Lednicky Q_{inv} Theoretical Correlation Function");
  analysis->AddParticleFunction(wqinvcfP);

  /************************************************************/
  /****   Q INV Correlation Function   ************************/
  /************************************************************/
   

  AliHBTQInvCorrelFctn * qinvcfT = new AliHBTQInvCorrelFctn(100,qinvmax);
  AliHBTQInvCorrelFctn * qinvcfP = new AliHBTQInvCorrelFctn(100,qinvmax);
  
  analysis->AddTrackFunction(qinvcfT);
  analysis->AddParticleFunction(qinvcfP);
  
  
  qinvcfP->Rename("qinvcfP","Particle (simulated) Qinv CF \\pi^{+} \\pi^{+}");
  qinvcfT->Rename("qinvcfT","Track (recontructed) Qinv CF \\pi^{+} \\pi^{+}");
  
  /************************************************************/
  /****   Q OUT Correlation Function   ************************/
  /************************************************************/
  
  AliHBTQOutLCMSCorrelFctn* qoutP = new AliHBTQOutLCMSCorrelFctn();
  qoutP->Rename("qoutP","Particle (simulated) Q_{out} CF \\pi^{+} \\pi^{+}");
  AliHBTQOutLCMSCorrelFctn* qoutT = new AliHBTQOutLCMSCorrelFctn(); 
  qoutT->Rename("qoutT","Track (recontructed) Q_{out} CF \\pi^{+} \\pi^{+}");

  AliHBTPairCut *outPairCut = new AliHBTPairCut();
  outPairCut->SetQOutCMSLRange(0.0,0.15);
  outPairCut->SetQSideCMSLRange(0.0,0.02);
  outPairCut->SetQLongCMSLRange(0.0,0.02);
  qoutP->SetPairCut(outPairCut);
  qoutT->SetPairCut(outPairCut);

  //analysis->AddTrackFunction(qoutP);
  //analysis->AddParticleFunction(qoutT);

  /************************************************************/
  /****   Q SIDE Correlation Function   ***********************/
  /************************************************************/
  
  AliHBTQSideLCMSCorrelFctn* qsideP = new AliHBTQSideLCMSCorrelFctn(100,qinvmax); 
  qsideP->Rename("qsideP","Particle (simulated) Q_{side} CF \\pi^{+} \\pi^{+}");
  AliHBTQSideLCMSCorrelFctn* qsideT = new AliHBTQSideLCMSCorrelFctn(100,qinvmax); 
  qsideT->Rename("qsideT","Track (recontructed) Q_{side} CF \\pi^{+} \\pi^{+}");

  AliHBTPairCut *sidePairCut = new AliHBTPairCut();
  sidePairCut->SetQOutCMSLRange(0.0,0.02);
  sidePairCut->SetQSideCMSLRange(0.0,0.15);
  sidePairCut->SetQLongCMSLRange(0.0,0.02);
  qsideP->SetPairCut(sidePairCut);
  qsideT->SetPairCut(sidePairCut);

  //analysis->AddTrackFunction(qsideP);
  //analysis->AddParticleFunction(qsideT);

  /************************************************************/
  /****   Q LONG Correlation Function   ***********************/
  /************************************************************/
    
  AliHBTQLongLCMSCorrelFctn* qlongP = new AliHBTQLongLCMSCorrelFctn(100,qinvmax);
  qlongP->Rename("qlongP","Particle (simulated) Q_{long} CF \\pi^{+} \\pi^{+}");
  AliHBTQLongLCMSCorrelFctn* qlongT = new AliHBTQLongLCMSCorrelFctn(100,qinvmax); 
  qlongT->Rename("qlongT","Track (recontructed) Q_{long} CF \\pi^{+} \\pi^{+}");

  AliHBTPairCut *longPairCut = new AliHBTPairCut();
  longPairCut->SetQOutCMSLRange(0.0,0.02);
  longPairCut->SetQSideCMSLRange(0.0,0.02);
  longPairCut->SetQLongCMSLRange(0.0,0.15);
  qlongP->SetPairCut(longPairCut);
  qlongT->SetPairCut(longPairCut);

  //analysis->AddTrackFunction(qlongT);
  //analysis->AddParticleFunction(qlongT);


  /************************************************************/
  /************************************************************/
  /************************************************************/
  /****   R E S O L U T I O N S         ***********************/
  /************************************************************/
  /************************************************************/
  /************************************************************/


  AliHBTQInvResolVsKtFctn qinvVsktF(200,1.0,0.0,300,0.015,-0.015);    //QInvLCMS  Res   Vs   Kt
  if (threeDcuts)  qinvVsktF.SetPairCut(paircut);
  //analysis->AddResolutionFunction(&qinvVsktF);

  AliHBTQOutResolVsKtFctn qoutVsktF(200,1.0,0.0,300,0.05,-0.05);    //QOutLCMS  Res   Vs   Kt
  if (threeDcuts)  qoutVsktF.SetPairCut(outPairCut);
  //analysis->AddResolutionFunction(&qoutVsktF);

  AliHBTQSideResolVsKtFctn qsideVsktF(200,1.0,0.0,300,0.015,-0.015);   //QSideLCMS Res   Vs   Kt
  if (threeDcuts)  qsideVsktF.SetPairCut(sidePairCut);
  //analysis->AddResolutionFunction(&qsideVsktF);

  AliHBTQLongResolVsKtFctn qlongVsktF(200,1.0,0.0,300,0.015,-0.015);   //QLongLCMS Res   Vs   Kt
  if (threeDcuts)  qlongVsktF.SetPairCut(longPairCut);
  //analysis->AddResolutionFunction(&qlongVsktF);

  AliHBTQInvResolVsQInvFctn qInvVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QInvLCMS Res   Vs   QInvLCMS
  //analysis->AddResolutionFunction(&qInvVsqinvF);

  AliHBTQOutResolVsQInvFctn qoutVsqinvF(200,qinvmax,qinvmin,300,0.05,-0.05);  //QOutLCMS  Res   Vs   QInvLCMS
  if (threeDcuts)  qoutVsqinvF.SetPairCut(outPairCut);
  //analysis->AddResolutionFunction(&qoutVsqinvF);

  AliHBTQSideResolVsQInvFctn qsideVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QSideLCMS Res   Vs   QInvLCMS
  if (threeDcuts)  qsideVsqinvF.SetPairCut(sidePairCut);
  //analysis->AddResolutionFunction(&qsideVsqinvF);

  AliHBTQLongResolVsQInvFctn qlongVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QLongLCMS Res   Vs   QInvLCMS
  if (threeDcuts)  qlongVsqinvF.SetPairCut(longPairCut);
  //analysis->AddResolutionFunction(&qlongVsqinvF);


  AliHBTPairThetaResolVsQInvFctn pairThetaVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
  if (threeDcuts) pairThetaVsqinv.SetPairCut(longPairCut);
  //analysis->AddResolutionFunction(&pairThetaVsqinv);

  AliHBTPairThetaResolVsKtFctn pairThetaVsKt(200,1.0,0.0,300,0.015,-0.015);
  if (threeDcuts) pairThetaVsKt.SetPairCut(longPairCut);
  //analysis->AddResolutionFunction(&pairThetaVsKt);

  AliHBTPairCut phipc;
  phipc.SetQLongCMSLRange(0.0,0.02);

  AliHBTPairPhiResolVsQInvFctn pairPhiVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
  if (threeDcuts) pairPhiVsqinv.SetPairCut(&phipc);
  //analysis->AddResolutionFunction(&pairPhiVsqinv);

  AliHBTPairPhiResolVsKtFctn pairPhiVsKt(200,1.0,0.0,300,0.015,-0.015);
  if (threeDcuts) pairPhiVsKt.SetPairCut(&phipc);
  //analysis->AddResolutionFunction(&pairPhiVsKt);


  AliHBTQOutResolVsQOutFctn qoutVsqoutF;  //QOutLCMS  Res   Vs   QOut
  if (threeDcuts) qoutVsqoutF.SetPairCut(outPairCut);
  //analysis->AddResolutionFunction(&qoutVsqoutF);

  AliHBTQSideResolVsQSideFctn qsideVsqsideF;//QSideLCMS Res   Vs   QSide
  if (threeDcuts) qsideVsqsideF.SetPairCut(sidePairCut);
  //analysis->AddResolutionFunction(&qsideVsqsideF);

  
  AliHBTQLongResolVsQLongFctn qlongVsqlongF;//QLongLCMS Res   Vs   QLong
  if (threeDcuts) qlongVsqlongF.SetPairCut(longPairCut);
  //analysis->AddResolutionFunction(&qlongVsqlongF);


  /************************************************************/
  /************************************************************/
  /************************************************************/
  /****   D I S T R B U T I O N S         ***********************/
  /************************************************************/
  /************************************************************/
  /************************************************************/

  /************************************************************/
  /****   OUT VS KT                 ***********************/
  /************************************************************/
  AliHBTQOutDistributionVsKtFctn qoutVsKtDistP;
  qoutVsKtDistP.Rename("qoutVsKtDistP",qoutVsKtDistP.GetTitle());
  AliHBTQOutDistributionVsKtFctn qoutVsKtDistT;
  qoutVsKtDistT.Rename("qoutVsKtDistT",qoutVsKtDistT.GetTitle());
  if (threeDcuts)  qoutVsKtDistP.SetPairCut(outPairCut);
  if (threeDcuts)  qoutVsKtDistT.SetPairCut(outPairCut);
  //analysis->AddParticleFunction(&qoutVsKtDistP);
  //analysis->AddTrackFunction(&qoutVsKtDistT);
  /************************************************************/
  /****   Side VS KT                 ***********************/
  /************************************************************/
  AliHBTQSideDistributionVsKtFctn qSideVsKtDistP;
  qSideVsKtDistP.Rename("qSideVsKtDistP",qSideVsKtDistP.GetTitle());
  AliHBTQSideDistributionVsKtFctn qSideVsKtDistT;
  qSideVsKtDistT.Rename("qSideVsKtDistT",qSideVsKtDistT.GetTitle());
  if (threeDcuts)  qSideVsKtDistP.SetPairCut(sidePairCut);
  if (threeDcuts)  qSideVsKtDistT.SetPairCut(sidePairCut);
  //analysis->AddParticleFunction(&qSideVsKtDistP);
  //analysis->AddTrackFunction(&qSideVsKtDistT);
  /************************************************************/
  /****   Long VS KT                 ***********************/
  /************************************************************/
  AliHBTQLongDistributionVsKtFctn qLongVsKtDistP;
  qLongVsKtDistP.Rename("qLongVsKtDistP",qLongVsKtDistP.GetTitle());
  AliHBTQLongDistributionVsKtFctn qLongVsKtDistT;
  qLongVsKtDistT.Rename("qLongVsKtDistT",qLongVsKtDistT.GetTitle());
  if (threeDcuts)  qLongVsKtDistP.SetPairCut(longPairCut);
  if (threeDcuts)  qLongVsKtDistT.SetPairCut(longPairCut);
  //analysis->AddParticleFunction(&qLongVsKtDistP);
  //analysis->AddTrackFunction(&qLongVsKtDistT);

  /************************************************************/
  /************************************************************/
  /************************************************************/
  /****    M O N I T O R  R E S O L U T I O N S   *************/
  /************************************************************/
  /************************************************************/
  /************************************************************/
  
  AliHBTMonPxResolutionVsPtFctn PxVsPtResF(200,1.0,0.0,300,0.05,-0.05);
  //analysis->AddParticleAndTrackMonitorFunction(&PxVsPtResF);
  
  AliHBTMonPyResolutionVsPtFctn PyVsPtResF(200,1.0,0.0,300,0.05,-0.05);
  //analysis->AddParticleAndTrackMonitorFunction(&PyVsPtResF);
  
  AliHBTMonPzResolutionVsPtFctn PzVsPtResF(200,1.0,0.0,300,0.05,-0.05);
  //analysis->AddParticleAndTrackMonitorFunction(&PzVsPtResF);

  AliHBTMonPResolutionVsPtFctn PVsPtResF(200,1.0,0.0,300,0.05,-0.05);
  //analysis->AddParticleAndTrackMonitorFunction(&PVsPtResF);

  AliHBTMonPtResolutionVsPtFctn PtVsPtResF(400,4.0,0.0,300,0.07,-0.07);
  analysis->AddParticleAndTrackMonitorFunction(&PtVsPtResF);

  AliHBTMonPhiResolutionVsPtFctn PhiVsPtResF(200,1.0,0.0,300,0.02,-0.02);
  analysis->AddParticleAndTrackMonitorFunction(&PhiVsPtResF);

  AliHBTMonThetaResolutionVsPtFctn ThetaVsPtResF(200,1.0,0.0,300,0.02,-0.02);
  analysis->AddParticleAndTrackMonitorFunction(&ThetaVsPtResF);
  
  
  /************************************************************/
  /****   P R O C E S S                 ***********************/
  /************************************************************/
  analysis->SetBufferSize(2);
  analysis->Process(processopt);
  
  TFile histoOutput(outfile,"recreate"); 
  analysis->WriteFunctions();
  histoOutput.Close();
  
  delete qinvcfP;
  delete qinvcfT;
  delete paircut;
  delete readerpartcut;
  if (dirs) 
   {
    dirs->SetOwner();
    delete dirs;
   }
  delete reader;
  delete analysis;
  
 }
Commit	Line	Data
f1e7ff24	1	#include "AliHBTAnalysis.h"
6e0d7f90	2	// #include "alles.h"
f1e7ff24	3	#include "AliHBTReader.h"
	4	#include "AliHBTReaderKineTree.h"
	5	#include "AliHBTReaderITSv2.h"
	6	#include "AliHBTReaderITSv1.h"
	7	#include "AliHBTReaderTPC.h"
	8	#include "AliHBTReaderInternal.h"
	9	#include "AliHBTParticleCut.h"
	10	#include "AliHBTEvent.h"
	11	#include "AliHBTPairCut.h"
	12	#include "AliHBTQResolutionFctns.h"
	13	#include "AliHBTQDistributionFctns.h"
	14	#include "AliHBTTwoTrackEffFctn.h"
	15	#include "AliHBTCorrelFctn.h"
	16	#include "TSystem.h"
	17	#include "TObjString.h"
	18	#include "TString.h"
6e0d7f90	19	#include "TPDGCode.h"
f1e7ff24	20	#include "AliHBTMonDistributionFctns.h"
	21	#include "AliHBTMonResolutionFctns.h"
	22
e477958c	23	void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
	24	Int_t first = -1,Int_t last = -1,
	25	char *outfile = "hbtanalysis.root")
0a6dde4c	26	{
9616170a	27
88378f71	28	AliHBTTrackPoints::SetDebug(0);
66d1d1a4	29	AliHBTParticle::SetDebug(0);
88378f71	30	AliLoader::SetDebug(0);
66d1d1a4	31
88378f71	32	AliHBTParticleCut c1;
	33	AliHBTParticleCut c2 = c1;
	34
0a6dde4c	35	//HBT Anlysis Macro
0a6dde4c	36	//Anlyzes TPC recontructed tracks and simulated particles that corresponds to them
e477958c	37
	38	//datatype defines type of data to be read
	39	// Kine - analyzes Kine Tree: simulated particles
	40	// TPC - analyzes TPC tracking + particles corresponding to these tracks
	41	// ITSv1 - analyzes ITSv1 ----------------------//--------------------------
	42	// ITSv2 - analyzes ITSv2 ----------------------//--------------------------
	43
	44	//processopt defines option passed to AliHBTAnlysis::Process method
	45	// default: TracksAndParticles - process both recontructed tracks and sim. particles corresponding to them
	46	// Tracks - process only recontructed tracks
	47	// Particles - process only simulated particles
	48
0a6dde4c	49	//Reads data from diroctories from first to last(including)
	50	// For examples if first=3 and last=5 it reads from
	51	// ./3/
	52	// ./4/
	53	// ./5/
	54	//if first or last is negative (or both), it reads from current directory
	55	//
	56	//these names I use when analysis is done directly from CASTOR files via RFIO
	57	// const char* basedir="rfio:/castor/cern.ch/user/s/skowron/";
	58	// const char* serie="standard";
	59	// const char* field = "0.4";
9616170a	60
0a6dde4c	61	const char* basedir=".";
	62	const char* serie="";
	63	const char* field = "";
	64
f1e7ff24	65	Bool_t threeDcuts = kFALSE;
f1e7ff24	66
0a6dde4c	67	// Dynamically link some shared libs
f1e7ff24	68	// gROOT->LoadMacro("loadlibs.C");
f1e7ff24	69	// loadlibs();
0a6dde4c	70
71094efa	71	gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libHBTAN");
0a6dde4c	72
	73
	74	/***********************************************************/
	75	//Create analysis and reader
	76	AliHBTAnalysis * analysis = new AliHBTAnalysis();
66d1d1a4	77	analysis->SetCutsOnTracks();
0a6dde4c	78
0a6dde4c	79	AliHBTReader* reader;
e477958c	80	Int_t kine = strcmp(datatype,"Kine");
6e0d7f90	81	Int_t ESD = strcmp(datatype,"ESD");
e477958c	82	Int_t TPC = strcmp(datatype,"TPC");
	83	Int_t ITSv1 = strcmp(datatype,"ITSv1");
	84	Int_t ITSv2 = strcmp(datatype,"ITSv2");
91923a38	85	Int_t intern = strcmp(datatype,"Intern");
91923a38	86
e477958c	87	if(!kine)
	88	{
	89	reader = new AliHBTReaderKineTree();
	90	processopt="Particles"; //this reader by definition reads only simulated particles
	91	}
6e0d7f90	92	else if(!ESD)
	93	{
	94	AliHBTReaderESD* esdreader = new AliHBTReaderESD();
	95	esdreader->ReadParticles(kTRUE);
88378f71	96	esdreader->SetNumberOfTrackPoints(5,30.);
88378f71	97	esdreader->SetClusterMap();
6e0d7f90	98	reader = esdreader;
6e0d7f90	99	}
e477958c	100	else if(!TPC)
	101	{
	102	reader = new AliHBTReaderTPC();
66d1d1a4	103	//((AliHBTReaderTPC*)reader)->SetNumberOfTrackPoints(5,30.);
66d1d1a4	104	((AliHBTReaderTPC*)reader)->SetClusterMap();
e477958c	105	}
	106	else if(!ITSv1)
	107	{
	108	reader = new AliHBTReaderITSv1();
	109	}
	110	else if(!ITSv2)
	111	{
	112	reader = new AliHBTReaderITSv2();
	113	}
91923a38	114	else if(!intern)
91923a38	115	{
6e0d7f90	116	reader = new AliHBTReaderInternal("data.root");
91923a38	117	}
e477958c	118	else
	119	{
	120	cerr<<"Option "<<datatype<<" not recognized. Exiting"<<endl;
	121	return;
	122	}
	123
0a6dde4c	124	TObjArray* dirs=0;
e477958c	125	if ( (first >= 0) && (last>=0) && ( (last-first)>=0 ) )
	126	{//read from many dirs dirs
	127	char buff[50];
	128	dirs = new TObjArray(last-first+1);
	129	for (Int_t i = first; i<=last; i++)
	130	{
	131	sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i);
	132	TObjString *odir= new TObjString(buff);
	133	dirs->Add(odir);
	134	}
0a6dde4c	135	}
e477958c	136	reader->SetDirs(dirs);
0a6dde4c	137
0a6dde4c	138	/***********************************************************/
f1e7ff24	139	/*** R E A D E R *************************************/
f1e7ff24	140	/***********************************************************/
0a6dde4c	141
	142	//we want have only low pt pi+ so set a cut to reader
	143	AliHBTParticleCut* readerpartcut= new AliHBTParticleCut();
6e0d7f90	144	readerpartcut->SetPtRange(0.0,10.0);
0a6dde4c	145	readerpartcut->SetPID(kPiPlus);
	146	reader->AddParticleCut(readerpartcut);//read this particle type with this cut
	147
6e0d7f90	148	// readerpartcut->SetPtRange(0.0,10.0);
	149	// readerpartcut->SetPID(kPiMinus);
	150	// reader->AddParticleCut(readerpartcut);//read this particle type with this cut
	151
0a6dde4c	152	analysis->SetReader(reader);
f1e7ff24	153
9616170a	154	AliHBTPairCut *paircut = new AliHBTPairCut();
	155	Float_t qinvmin = 0.0;
	156	Float_t qinvmax = 0.05;//50MeV
	157	paircut->SetQInvRange(qinvmin,qinvmax);
	158
66d1d1a4	159	// paircut->SetAvSeparationRange(10.); //AntiMerging
66d1d1a4	160	paircut->SetClusterOverlapRange(-.5,1.0);//AntiSplitting
9616170a	161
	162	// AliHBTParticleCut* partcut= new AliHBTParticleCut();
	163	// partcut->SetPID(kPiPlus);
	164	// paircut->SetFirstPartCut(partcut);
	165	// partcut->SetPID(kPiMinus);
	166	// paircut->SetSecondPartCut(partcut);
	167
	168	analysis->SetGlobalPairCut(paircut);
9616170a	169
16f9289f	170	/***********************************************************/
	171	/*** W E I G H T S ****************************/
	172	/***********************************************************/
	173
	174	AliHBTLLWeights::Instance()->SetParticlesTypes(kPiPlus,kPiPlus);
	175	AliHBTLLWeights::Instance()->SetColoumb(kFALSE);
	176	AliHBTLLWeights::Instance()->SetStrongInterSwitch(kFALSE);
	177
	178	AliHBTLLWeights::Instance()->SetRandomPosition(kFALSE);
	179	//AliHBTLLWeights::Instance()->SetR1dw(8.0);
	180
	181	AliHBTLLWeights::Instance()->Init();
	182	AliHBTLLWeights::Instance()->Set();
	183
	184	//example function
9616170a	185	AliHBTWeightTheorQInvFctn *wqinvcfP = new AliHBTWeightTheorQInvFctn(100,qinvmax);
	186	wqinvcfP->SetNumberOfBinsToScale(30);
	187	wqinvcfP->Rename("wqinvcfP","Lednicky Q_{inv} Theoretical Correlation Function");
16f9289f	188	analysis->AddParticleFunction(wqinvcfP);
16f9289f	189
f1e7ff24	190	/************************************************************/
f1e7ff24	191	/** Q INV Correlation Function **********************/
0a6dde4c	192	/************************************************************/
6e0d7f90	193
0a6dde4c	194
f1e7ff24	195	AliHBTQInvCorrelFctn * qinvcfT = new AliHBTQInvCorrelFctn(100,qinvmax);
f1e7ff24	196	AliHBTQInvCorrelFctn * qinvcfP = new AliHBTQInvCorrelFctn(100,qinvmax);
0a6dde4c	197
	198	analysis->AddTrackFunction(qinvcfT);
	199	analysis->AddParticleFunction(qinvcfP);
	200
0a6dde4c	201
	202	qinvcfP->Rename("qinvcfP","Particle (simulated) Qinv CF \\pi^{+} \\pi^{+}");
	203	qinvcfT->Rename("qinvcfT","Track (recontructed) Qinv CF \\pi^{+} \\pi^{+}");
91923a38	204
f1e7ff24	205	/************************************************************/
	206	/** Q OUT Correlation Function **********************/
	207	/************************************************************/
	208
78d7c6d3	209	AliHBTQOutLCMSCorrelFctn* qoutP = new AliHBTQOutLCMSCorrelFctn();
91923a38	210	qoutP->Rename("qoutP","Particle (simulated) Q_{out} CF \\pi^{+} \\pi^{+}");
78d7c6d3	211	AliHBTQOutLCMSCorrelFctn* qoutT = new AliHBTQOutLCMSCorrelFctn();
91923a38	212	qoutT->Rename("qoutT","Track (recontructed) Q_{out} CF \\pi^{+} \\pi^{+}");
	213
	214	AliHBTPairCut *outPairCut = new AliHBTPairCut();
	215	outPairCut->SetQOutCMSLRange(0.0,0.15);
	216	outPairCut->SetQSideCMSLRange(0.0,0.02);
	217	outPairCut->SetQLongCMSLRange(0.0,0.02);
	218	qoutP->SetPairCut(outPairCut);
	219	qoutT->SetPairCut(outPairCut);
f1e7ff24	220
6e0d7f90	221	//analysis->AddTrackFunction(qoutP);
6e0d7f90	222	//analysis->AddParticleFunction(qoutT);
f1e7ff24	223
	224	/************************************************************/
	225	/** Q SIDE Correlation Function *********************/
	226	/************************************************************/
91923a38	227
78d7c6d3	228	AliHBTQSideLCMSCorrelFctn* qsideP = new AliHBTQSideLCMSCorrelFctn(100,qinvmax);
91923a38	229	qsideP->Rename("qsideP","Particle (simulated) Q_{side} CF \\pi^{+} \\pi^{+}");
78d7c6d3	230	AliHBTQSideLCMSCorrelFctn* qsideT = new AliHBTQSideLCMSCorrelFctn(100,qinvmax);
91923a38	231	qsideT->Rename("qsideT","Track (recontructed) Q_{side} CF \\pi^{+} \\pi^{+}");
	232
	233	AliHBTPairCut *sidePairCut = new AliHBTPairCut();
	234	sidePairCut->SetQOutCMSLRange(0.0,0.02);
	235	sidePairCut->SetQSideCMSLRange(0.0,0.15);
	236	sidePairCut->SetQLongCMSLRange(0.0,0.02);
	237	qsideP->SetPairCut(sidePairCut);
	238	qsideT->SetPairCut(sidePairCut);
0a6dde4c	239
6e0d7f90	240	//analysis->AddTrackFunction(qsideP);
6e0d7f90	241	//analysis->AddParticleFunction(qsideT);
f1e7ff24	242
	243	/************************************************************/
	244	/** Q LONG Correlation Function *********************/
	245	/************************************************************/
91923a38	246
78d7c6d3	247	AliHBTQLongLCMSCorrelFctn* qlongP = new AliHBTQLongLCMSCorrelFctn(100,qinvmax);
91923a38	248	qlongP->Rename("qlongP","Particle (simulated) Q_{long} CF \\pi^{+} \\pi^{+}");
78d7c6d3	249	AliHBTQLongLCMSCorrelFctn* qlongT = new AliHBTQLongLCMSCorrelFctn(100,qinvmax);
91923a38	250	qlongT->Rename("qlongT","Track (recontructed) Q_{long} CF \\pi^{+} \\pi^{+}");
	251
	252	AliHBTPairCut *longPairCut = new AliHBTPairCut();
	253	longPairCut->SetQOutCMSLRange(0.0,0.02);
	254	longPairCut->SetQSideCMSLRange(0.0,0.02);
	255	longPairCut->SetQLongCMSLRange(0.0,0.15);
	256	qlongP->SetPairCut(longPairCut);
	257	qlongT->SetPairCut(longPairCut);
	258
6e0d7f90	259	//analysis->AddTrackFunction(qlongT);
6e0d7f90	260	//analysis->AddParticleFunction(qlongT);
f1e7ff24	261
	262
	263	/************************************************************/
	264	/************************************************************/
	265	/************************************************************/
	266	/** R E S O L U T I O N S *********************/
	267	/************************************************************/
	268	/************************************************************/
	269	/************************************************************/
	270
	271
	272
78d7c6d3	273	AliHBTQInvResolVsKtFctn qinvVsktF(200,1.0,0.0,300,0.015,-0.015); //QInvLCMS Res Vs Kt
f1e7ff24	274	if (threeDcuts) qinvVsktF.SetPairCut(paircut);
6e0d7f90	275	//analysis->AddResolutionFunction(&qinvVsktF);
f1e7ff24	276
78d7c6d3	277	AliHBTQOutResolVsKtFctn qoutVsktF(200,1.0,0.0,300,0.05,-0.05); //QOutLCMS Res Vs Kt
f1e7ff24	278	if (threeDcuts) qoutVsktF.SetPairCut(outPairCut);
6e0d7f90	279	//analysis->AddResolutionFunction(&qoutVsktF);
f1e7ff24	280
78d7c6d3	281	AliHBTQSideResolVsKtFctn qsideVsktF(200,1.0,0.0,300,0.015,-0.015); //QSideLCMS Res Vs Kt
f1e7ff24	282	if (threeDcuts) qsideVsktF.SetPairCut(sidePairCut);
6e0d7f90	283	//analysis->AddResolutionFunction(&qsideVsktF);
f1e7ff24	284
78d7c6d3	285	AliHBTQLongResolVsKtFctn qlongVsktF(200,1.0,0.0,300,0.015,-0.015); //QLongLCMS Res Vs Kt
f1e7ff24	286	if (threeDcuts) qlongVsktF.SetPairCut(longPairCut);
6e0d7f90	287	//analysis->AddResolutionFunction(&qlongVsktF);
f1e7ff24	288
78d7c6d3	289	AliHBTQInvResolVsQInvFctn qInvVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QInvLCMS Res Vs QInvLCMS
6e0d7f90	290	//analysis->AddResolutionFunction(&qInvVsqinvF);
f1e7ff24	291
78d7c6d3	292	AliHBTQOutResolVsQInvFctn qoutVsqinvF(200,qinvmax,qinvmin,300,0.05,-0.05); //QOutLCMS Res Vs QInvLCMS
f1e7ff24	293	if (threeDcuts) qoutVsqinvF.SetPairCut(outPairCut);
6e0d7f90	294	//analysis->AddResolutionFunction(&qoutVsqinvF);
f1e7ff24	295
78d7c6d3	296	AliHBTQSideResolVsQInvFctn qsideVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QSideLCMS Res Vs QInvLCMS
f1e7ff24	297	if (threeDcuts) qsideVsqinvF.SetPairCut(sidePairCut);
6e0d7f90	298	//analysis->AddResolutionFunction(&qsideVsqinvF);
f1e7ff24	299
78d7c6d3	300	AliHBTQLongResolVsQInvFctn qlongVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QLongLCMS Res Vs QInvLCMS
f1e7ff24	301	if (threeDcuts) qlongVsqinvF.SetPairCut(longPairCut);
6e0d7f90	302	//analysis->AddResolutionFunction(&qlongVsqinvF);
f1e7ff24	303
	304
	305	AliHBTPairThetaResolVsQInvFctn pairThetaVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
	306	if (threeDcuts) pairThetaVsqinv.SetPairCut(longPairCut);
6e0d7f90	307	//analysis->AddResolutionFunction(&pairThetaVsqinv);
f1e7ff24	308
	309	AliHBTPairThetaResolVsKtFctn pairThetaVsKt(200,1.0,0.0,300,0.015,-0.015);
	310	if (threeDcuts) pairThetaVsKt.SetPairCut(longPairCut);
6e0d7f90	311	//analysis->AddResolutionFunction(&pairThetaVsKt);
f1e7ff24	312
	313	AliHBTPairCut phipc;
	314	phipc.SetQLongCMSLRange(0.0,0.02);
	315
	316	AliHBTPairPhiResolVsQInvFctn pairPhiVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
	317	if (threeDcuts) pairPhiVsqinv.SetPairCut(&phipc);
6e0d7f90	318	//analysis->AddResolutionFunction(&pairPhiVsqinv);
f1e7ff24	319
	320	AliHBTPairPhiResolVsKtFctn pairPhiVsKt(200,1.0,0.0,300,0.015,-0.015);
	321	if (threeDcuts) pairPhiVsKt.SetPairCut(&phipc);
6e0d7f90	322	//analysis->AddResolutionFunction(&pairPhiVsKt);
f1e7ff24	323
f1e7ff24	324
78d7c6d3	325	AliHBTQOutResolVsQOutFctn qoutVsqoutF; //QOutLCMS Res Vs QOut
f1e7ff24	326	if (threeDcuts) qoutVsqoutF.SetPairCut(outPairCut);
6e0d7f90	327	//analysis->AddResolutionFunction(&qoutVsqoutF);
f1e7ff24	328
78d7c6d3	329	AliHBTQSideResolVsQSideFctn qsideVsqsideF;//QSideLCMS Res Vs QSide
f1e7ff24	330	if (threeDcuts) qsideVsqsideF.SetPairCut(sidePairCut);
6e0d7f90	331	//analysis->AddResolutionFunction(&qsideVsqsideF);
f1e7ff24	332
91923a38	333
78d7c6d3	334	AliHBTQLongResolVsQLongFctn qlongVsqlongF;//QLongLCMS Res Vs QLong
f1e7ff24	335	if (threeDcuts) qlongVsqlongF.SetPairCut(longPairCut);
6e0d7f90	336	//analysis->AddResolutionFunction(&qlongVsqlongF);
f1e7ff24	337
	338
	339	/************************************************************/
	340	/************************************************************/
	341	/************************************************************/
	342	/** D I S T R B U T I O N S *********************/
	343	/************************************************************/
	344	/************************************************************/
	345	/************************************************************/
91923a38	346
f1e7ff24	347	/************************************************************/
	348	/** OUT VS KT *********************/
	349	/************************************************************/
	350	AliHBTQOutDistributionVsKtFctn qoutVsKtDistP;
	351	qoutVsKtDistP.Rename("qoutVsKtDistP",qoutVsKtDistP.GetTitle());
	352	AliHBTQOutDistributionVsKtFctn qoutVsKtDistT;
	353	qoutVsKtDistT.Rename("qoutVsKtDistT",qoutVsKtDistT.GetTitle());
	354	if (threeDcuts) qoutVsKtDistP.SetPairCut(outPairCut);
	355	if (threeDcuts) qoutVsKtDistT.SetPairCut(outPairCut);
6e0d7f90	356	//analysis->AddParticleFunction(&qoutVsKtDistP);
6e0d7f90	357	//analysis->AddTrackFunction(&qoutVsKtDistT);
f1e7ff24	358	/************************************************************/
	359	/** Side VS KT *********************/
	360	/************************************************************/
	361	AliHBTQSideDistributionVsKtFctn qSideVsKtDistP;
	362	qSideVsKtDistP.Rename("qSideVsKtDistP",qSideVsKtDistP.GetTitle());
	363	AliHBTQSideDistributionVsKtFctn qSideVsKtDistT;
	364	qSideVsKtDistT.Rename("qSideVsKtDistT",qSideVsKtDistT.GetTitle());
	365	if (threeDcuts) qSideVsKtDistP.SetPairCut(sidePairCut);
	366	if (threeDcuts) qSideVsKtDistT.SetPairCut(sidePairCut);
6e0d7f90	367	//analysis->AddParticleFunction(&qSideVsKtDistP);
6e0d7f90	368	//analysis->AddTrackFunction(&qSideVsKtDistT);
f1e7ff24	369	/************************************************************/
	370	/** Long VS KT *********************/
	371	/************************************************************/
	372	AliHBTQLongDistributionVsKtFctn qLongVsKtDistP;
	373	qLongVsKtDistP.Rename("qLongVsKtDistP",qLongVsKtDistP.GetTitle());
	374	AliHBTQLongDistributionVsKtFctn qLongVsKtDistT;
	375	qLongVsKtDistT.Rename("qLongVsKtDistT",qLongVsKtDistT.GetTitle());
	376	if (threeDcuts) qLongVsKtDistP.SetPairCut(longPairCut);
	377	if (threeDcuts) qLongVsKtDistT.SetPairCut(longPairCut);
6e0d7f90	378	//analysis->AddParticleFunction(&qLongVsKtDistP);
6e0d7f90	379	//analysis->AddTrackFunction(&qLongVsKtDistT);
f1e7ff24	380
	381	/************************************************************/
	382	/************************************************************/
	383	/************************************************************/
	384	/** M O N I T O R R E S O L U T I O N S ***********/
	385	/************************************************************/
	386	/************************************************************/
	387	/************************************************************/
	388
	389	AliHBTMonPxResolutionVsPtFctn PxVsPtResF(200,1.0,0.0,300,0.05,-0.05);
6e0d7f90	390	//analysis->AddParticleAndTrackMonitorFunction(&PxVsPtResF);
f1e7ff24	391
f1e7ff24	392	AliHBTMonPyResolutionVsPtFctn PyVsPtResF(200,1.0,0.0,300,0.05,-0.05);
6e0d7f90	393	//analysis->AddParticleAndTrackMonitorFunction(&PyVsPtResF);
f1e7ff24	394
f1e7ff24	395	AliHBTMonPzResolutionVsPtFctn PzVsPtResF(200,1.0,0.0,300,0.05,-0.05);
6e0d7f90	396	//analysis->AddParticleAndTrackMonitorFunction(&PzVsPtResF);
f1e7ff24	397
f1e7ff24	398	AliHBTMonPResolutionVsPtFctn PVsPtResF(200,1.0,0.0,300,0.05,-0.05);
6e0d7f90	399	//analysis->AddParticleAndTrackMonitorFunction(&PVsPtResF);
f1e7ff24	400
6e0d7f90	401	AliHBTMonPtResolutionVsPtFctn PtVsPtResF(400,4.0,0.0,300,0.07,-0.07);
f1e7ff24	402	analysis->AddParticleAndTrackMonitorFunction(&PtVsPtResF);
	403
	404	AliHBTMonPhiResolutionVsPtFctn PhiVsPtResF(200,1.0,0.0,300,0.02,-0.02);
	405	analysis->AddParticleAndTrackMonitorFunction(&PhiVsPtResF);
	406
	407	AliHBTMonThetaResolutionVsPtFctn ThetaVsPtResF(200,1.0,0.0,300,0.02,-0.02);
	408	analysis->AddParticleAndTrackMonitorFunction(&ThetaVsPtResF);
	409
	410
	411
	412
	413	/************************************************************/
	414	/** P R O C E S S *********************/
	415	/************************************************************/
6e0d7f90	416	analysis->SetBufferSize(2);
91923a38	417	analysis->Process(processopt);
91923a38	418
0a6dde4c	419	TFile histoOutput(outfile,"recreate");
	420	analysis->WriteFunctions();
	421	histoOutput.Close();
	422
	423	delete qinvcfP;
	424	delete qinvcfT;
d4ebbc5f	425	delete paircut;
0a6dde4c	426	delete readerpartcut;
	427	if (dirs)
	428	{
	429	dirs->SetOwner();
	430	delete dirs;
	431	}
	432	delete reader;
	433	delete analysis;
	434
	435	}
	436