PWG2/UNICOR/AliDLoop.cxx

   1 //Author: Dariusz Miskowiec 2007
   2
   3 //=============================================================================
   4 // simple event loop manager
   5 // First, events are classified according to centrality, reaction plane angle,
   6 // and z-vertex. Then, single particles and true pairs are processed. Finally,
   7 // event mixing is done within event classes. For this, the tree and the event
   8 // passed as arguments are cloned.
   9 // The class is, at present, deliberately left in a macro-like design, with all
  10 // the essential things happening inside the function AliDLoop::Run. This is
  11 // subject to future polishing if the class proves useful.
  12 //=============================================================================
  13
  14 #include <TROOT.h>
  15 #include <TMath.h>
  16 #include <TTree.h>
  17 #include <TSystem.h>
  18 #include <TStopwatch.h>
  19 #include "AliDLoop.h"
  20 #include "AliDAnalGlobal.h"
  21 #include "AliDAnalSingle.h"
  22 #include "AliDAnalCorrel.h"
  23 #include "AliDAnalPtfluc.h"
  24 #include "AliDEvent.h"
  25
  26 ClassImp(AliDLoop)
  27
  28 //=============================================================================
  29 AliDLoop::AliDLoop(TTree *tr, AliDEvent *ev0, char *outfil):
  30   TObject(),
  31   fTree0(tr), fTree1(0x0),
  32   fEv0(ev0), fEv1(0x0),
  33   fOutputFilename(outfil)
  34 {
  35   // constructor
  36   // Clone the tree such that during the event mixing the two events get
  37   // populated from two different trees. This costs some memory (37 MB for
  38   // ceres) but allows to avoid switching one tree between two different
  39   // events, especially difficult for alice because alice trees do not like
  40   // to be attached more than once.
  41
  42   fTree1 = (TTree*) tr->Clone();
  43   fEv1 = (AliDEvent*) ev0->Clone();
  44   fEv0->AttachTree(fTree0);
  45   fEv1->AttachTree(fTree1);
  46 }
  47 //=============================================================================
  48 AliDLoop::AliDLoop(const AliDLoop &loop) :
  49   TObject(loop), fTree0(loop.fTree0), fTree1(loop.fTree1), fEv0(loop.fEv0),
  50   fEv1(loop.fEv1), fOutputFilename(loop.fOutputFilename)
  51 {
  52   //copy constructor, shallow copy, just to make the compiler happy
  53 }
  54 //=============================================================================
  55 AliDLoop &AliDLoop::operator=(const AliDLoop &source)
  56 {
  57   // substitution operator, shallow copy, just to make the compiler happy
  58
  59   TObject::operator=(source);
  60   fTree0 = source.fTree0;
  61   fTree1 = source.fTree1;
  62   fEv0 = source.fEv0;
  63   fEv1 = source.fEv1;
  64   fOutputFilename = source.fOutputFilename;
  65   return *this;
  66 }
  67 //=============================================================================
  68 Int_t AliDLoop::Mem() const
  69 {
  70   // return memory in MB occupied by this process
  71
  72   ProcInfo_t info;
  73   gSystem->GetProcInfo(&info);
  74   Int_t memmb = info.fMemVirtual/1024;
  75   return memmb;
  76 }
  77 //=============================================================================
  78 void AliDLoop::Run(int nwish) const
  79 {
  80   // process nwish events from fTree0; nwish=-1 (default) means process all
  81
  82   // control parameters
  83
  84
  85   //  const int kcen = 10; // number of centrality bins for mixing
  86   //  const int kphi = 8; // number of phi bins for mixing
  87   //  const int kzve = 13; // number of z-vertex bins for mixing
  88   const int kcen = 1; // number of centrality bins for mixing
  89   const int kphi = 1; // number of phi bins for mixing
  90   const int kzve = 1; // number of z-vertex bins for mixing
  91   int mixingFactor = 1; // number of wished event pairs/number of events
  92   int mixingStep[10] = {2,3,5,7,11,13,17,19,23,29};
  93   if (mixingFactor>=10) {printf("mixing factor too high\n"); exit(-1);}
  94   int minn = 20*mixingFactor; // minimum number of events per bin
  95   int pr = !gROOT->IsBatch(); // more  printing in interactive mode
  96
  97   // initialize analysis
  98
  99   Double_t etami = fEv0->Etamin();
 100   Double_t etama = fEv0->Etamax();
 101   AliDAnalGlobal *dag = new AliDAnalGlobal("dag");
 102   AliDAnalSingle *all = new AliDAnalSingle("all",etami,etama,0);
 103   AliDAnalSingle *pim = new AliDAnalSingle("pim",etami,etama,-211);
 104   AliDAnalSingle *pip = new AliDAnalSingle("pip",etami,etama, 211);
 105   AliDAnalCorrel *cnn = new AliDAnalCorrel("cnn",etami,etama,-211,-211);
 106   AliDAnalCorrel *cpp = new AliDAnalCorrel("cpp",etami,etama, 211, 211);
 107   AliDAnalPtfluc *ptf = new AliDAnalPtfluc("ptf",-211,-211);
 108
 109   // initialize number of events etc
 110
 111   printf("\ndetermining number of events in chain...");
 112   int nmax = fTree0->GetEntries();
 113   int nact = nmax;
 114   if (nwish>-1 && nwish<nmax) nact = nwish;
 115   printf("analyzing %d events of %d\n\n",nact,nmax);
 116   TStopwatch sto;
 117
 118   sto.Reset();
 119   sto.Start();
 120
 121   // find groups of of similar events
 122   // first allocate nact entries to each array; later reduce
 123
 124   printf("%4d MB; booking   event index arrays\n",Mem());
 125   TArrayI *ar[kcen][kphi][kzve];
 126   int nar[kcen][kphi][kzve];
 127   for (int i=0; i<kcen; i++)
 128     for (int j=0; j<kphi; j++)
 129       for (int k=0; k<kzve; k++) {
 130         ar[i][j][k] = new TArrayI(nact);
 131         nar[i][j][k] = 0;
 132       }
 133   printf("%4d MB; filling   event index arrays\n",Mem());
 134   for (int i=0; i<nact; i++) {
 135     fTree0->GetEvent(i);
 136     if (!fEv0->Good()) continue;
 137     int icen = (int) (kcen*fEv0->Centrality());
 138     int iphi = (int) (kphi*(fEv0->RPphi()/TMath::Pi()+1)/2);
 139     int izve = (int) (kzve*(fEv0->Zver()+1.0)/2.0);
 140     if (pr) printf("\revent %5d   %2d%2d%2d",i,icen,iphi,izve);
 141     if (icen<0 || icen>=kcen) continue;
 142     if (iphi<0 || iphi>=kphi) continue;
 143     if (izve<0 || izve>=kzve) continue;
 144     ar[icen][iphi][izve]->AddAt(i,nar[icen][iphi][izve]);
 145     nar[icen][iphi][izve]++;
 146   }
 147   if (pr) printf("\n");
 148
 149   printf("%4d MB; shrinking event index arrays\n",Mem());
 150   int nall = 0;
 151   int nsup = 0;
 152   for (int i=0; i<kcen; i++)
 153     for (int j=0; j<kphi; j++)
 154       for (int k=0; k<kzve; k++) {
 155         nall += nar[i][j][k];
 156         // suppress scarse bins
 157         if (nar[i][j][k]<minn) {
 158           nsup += nar[i][j][k];
 159           nar[i][j][k]=0;
 160         }
 161         ar[i][j][k]->Set(nar[i][j][k]);
 162       }
 163   printf("%4d MB; %d events sorted into bins\n",Mem(),nall);
 164   printf("%d events suppressed because bin occupancy was below %d\n",nsup,minn);
 165   printf("sorting events took %.1f s CPU %.1f s real %.1f ms real per event\n\n",
 166          sto.CpuTime(),sto.RealTime(),1000*sto.RealTime()/nall);
 167   if ((nall-=nsup) < 1) return;
 168
 169   // actual loop
 170
 171   printf("%4d MB; starting the main loop\n",Mem());
 172   sto.Reset();
 173   sto.Start();
 174   for (int i=0; i<kcen; i++)
 175     for (int j=0; j<kphi; j++)
 176       for (int k=0; k<kzve; k++) {
 177
 178         // true pairs
 179
 180         int nevents = ar[i][j][k]->GetSize(); // number of events in this bin
 181         for (int l=0; l<nevents; l++) {
 182           if (pr) printf("\revent %6d of %6d  %3d%3d%3d",l,nevents,i,j,k);
 183           int m = ar[i][j][k]->At(l);
 184           fTree0->GetEvent(m);
 185           dag->Process(fEv0);
 186           all->Process(fEv0);
 187           pim->Process(fEv0);
 188           pip->Process(fEv0);
 189           cnn->Process(0,fEv0,fEv0,0);
 190           cnn->Process(2,fEv0,fEv0,TMath::DegToRad()*180);
 191           cpp->Process(0,fEv0,fEv0,0);
 192           cpp->Process(2,fEv0,fEv0,TMath::DegToRad()*180);
 193           ptf->Process(0,fEv0,fEv0);
 194         }
 195         if (pr && nevents) printf("\n");
 196
 197         // mixed pairs
 198         // loop somewhat optimized to reduce jumping back and forth
 199
 200         for (int imix=0; imix<mixingFactor; imix++) {
 201           int step = mixingStep[imix];
 202           for (int phase=0; phase<step; phase++) for (int l=phase; l<nevents; l+=step) {
 203             int m0 = ar[i][j][k]->At(l);
 204             int next = (l+step)%nevents;
 205             int m1 = ar[i][j][k]->At(next);
 206             if (pr) printf("\rmixing %5d and %5d",m0,m1);
 207             fTree0->GetEvent(m0);
 208             fTree1->GetEvent(m1);
 209             //printf("   event0 %f   event1 %f\n",fEv0->ParticleP(0),fEv1->ParticleP(0));
 210             cnn->Process(1,fEv0,fEv1,0);
 211             cpp->Process(1,fEv0,fEv1,0);
 212             ptf->Process(1,fEv0,fEv1);
 213           }
 214         if (pr && nevents) printf("\r");
 215         }
 216       }
 217
 218   sto.Stop();
 219   printf("event loop took %.1f s CPU %.1f s real %.1f ms real per event\n\n",
 220          sto.CpuTime(),sto.RealTime(),1000*sto.RealTime()/nall);
 221
 222   // save analysis results
 223
 224   dag->Save(fOutputFilename.Data(),"recreate");
 225   all->Save(fOutputFilename.Data());
 226   pim->Save(fOutputFilename.Data());
 227   pip->Save(fOutputFilename.Data());
 228   cnn->Save(fOutputFilename.Data());
 229   cpp->Save(fOutputFilename.Data());
 230   ptf->Save(fOutputFilename.Data());
 231 }
 232 //=============================================================================