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