HBTAN/hbtanalysis.C

   1 void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
   2                 Int_t first = -1,Int_t last = -1,
   3                 char *outfile = "hbtanalysis.root")
   4  {
   5 //HBT Anlysis Macro
   6 //Anlyzes TPC recontructed tracks and simulated particles that corresponds to them
   7
   8 //datatype defines type of data to be read
   9 //  Kine  - analyzes Kine Tree: simulated particles
  10 //  TPC   - analyzes TPC   tracking + particles corresponding to these tracks
  11 //  ITSv1 - analyzes ITSv1 ----------------------//--------------------------
  12 //  ITSv2 - analyzes ITSv2 ----------------------//--------------------------
  13
  14 //processopt defines option passed to AliHBTAnlysis::Process method
  15 // default: TracksAndParticles - process both recontructed tracks and sim. particles corresponding to them
  16 //          Tracks - process only recontructed tracks
  17 //          Particles - process only simulated particles
  18
  19 //Reads data from diroctories from first to last(including)
  20 // For examples if first=3 and last=5 it reads from
  21 //  ./3/
  22 //  ./4/
  23 //  ./5/
  24 //if first or last is negative (or both), it reads from current directory
  25 //
  26 //these names I use when analysis is done directly from CASTOR files via RFIO
  27 //  const char* basedir="rfio:/castor/cern.ch/user/s/skowron/";
  28 //  const char* serie="standard";
  29 //  const char* field = "0.4";
  30   const char* basedir=".";
  31   const char* serie="";
  32   const char* field = "";
  33
  34   // Dynamically link some shared libs
  35   gROOT->LoadMacro("loadlibs.C");
  36   loadlibs();
  37
  38   gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libHBTAnalysis");
  39
  40
  41   /***********************************************************/
  42   //Create analysis and reader
  43   AliHBTAnalysis * analysis = new AliHBTAnalysis();
  44
  45   AliHBTReader* reader;
  46   Int_t kine = strcmp(datatype,"Kine");
  47   Int_t TPC = strcmp(datatype,"TPC");
  48   Int_t ITSv1 = strcmp(datatype,"ITSv1");
  49   Int_t ITSv2 = strcmp(datatype,"ITSv2");
  50
  51   if(!kine)
  52    {
  53     reader = new AliHBTReaderKineTree();
  54     processopt="Particles"; //this reader by definition reads only simulated particles
  55    }
  56   else if(!TPC)
  57    {
  58     reader = new AliHBTReaderTPC();
  59    }
  60   else if(!ITSv1)
  61    {
  62     reader = new AliHBTReaderITSv1();
  63    }
  64   else if(!ITSv2)
  65    {
  66     reader = new AliHBTReaderITSv2();
  67    }
  68   else
  69    {
  70     cerr<<"Option "<<datatype<<"  not recognized. Exiting"<<endl;
  71     return;
  72    }
  73
  74   TObjArray* dirs=0;
  75   if ( (first >= 0) && (last>=0) && ( (last-first)>=0 ) )
  76    {//read from many dirs dirs
  77      char buff[50];
  78      dirs = new TObjArray(last-first+1);
  79      for (Int_t i = first; i<=last; i++)
  80       {
  81         sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i);
  82         TObjString *odir= new TObjString(buff);
  83         dirs->Add(odir);
  84       }
  85     }
  86    reader->SetDirs(dirs);
  87
  88   /***********************************************************/
  89
  90   //we want have only low pt pi+ so set a cut to reader
  91   AliHBTParticleCut* readerpartcut= new AliHBTParticleCut();
  92   readerpartcut->SetPtRange(0.0,1.5);
  93   readerpartcut->SetPID(kPiPlus);
  94   reader->AddParticleCut(readerpartcut);//read this particle type with this cut
  95
  96   analysis->SetReader(reader);
  97   /************************************************************/
  98
  99   AliHBTPairCut *paircut = new AliHBTPairCut();
 100   paircut->SetQInvRange(0.0,0.20);
 101   analysis->SetGlobalPairCut(paircut);
 102
 103   AliHBTQInvCorrelFctn * qinvcfT= new AliHBTQInvCorrelFctn();
 104   AliHBTQInvCorrelFctn * qinvcfP= new AliHBTQInvCorrelFctn();
 105
 106   analysis->AddTrackFunction(qinvcfT);
 107   analysis->AddParticleFunction(qinvcfP);
 108
 109   analysis->Process(processopt);
 110
 111   qinvcfP->Rename("qinvcfP","Particle (simulated) Qinv CF \\pi^{+} \\pi^{+}");
 112   qinvcfT->Rename("qinvcfT","Track (recontructed) Qinv CF \\pi^{+} \\pi^{+}");
 113
 114   TFile histoOutput(outfile,"recreate");
 115   analysis->WriteFunctions();
 116   histoOutput.Close();
 117
 118   delete qinvcfP;
 119   delete qinvcfT;
 120   delete paircut;
 121   delete readerpartcut;
 122   if (dirs)
 123    {
 124     dirs->SetOwner();
 125     delete dirs;
 126    }
 127   delete reader;
 128   delete analysis;
 129
 130  }
 131