PWGCF/FLOW/macros/GenerateEventsOnTheFly.C

   1 #include "TStopwatch.h"
   2 #include "Riostream.h"
   3 #include "TFile.h"
   4 #include "TTimer.h"
   5 #include "TPythia6.h"
   6 #include "TMath.h"
   7
   8 /*
   9  * macro to generate flow events on the fly
  10  * if the macro is called with kFALSE as argument,
  11  * an example function of how to do analysis on the output is called
  12  *
  13  * at the end of the macro, an example function is given that can be used
  14  * to connect the on the fly events to the AliAnalysistwoParticleResonanceFlowTask
  15  *
  16  * Redmer Alexander Bertens (rbertens@cern.ch)
  17  * Carlos Eugenio Perez Lara
  18  * Andrea Dubla
  19  */
  20
  21 class AliFlowOnTheFlyGenerator;
  22 class AliAnalysisTwoParticleResonanceFlowTask;
  23 class AliFlowAnalysisWithQCumulants;
  24 // main function
  25 //
  26 // you can call this macro in different modes.
  27 // a) default mode will launch the event genetor and write the output to file
  28 // b) called with kFALSE the macro will read events from file and do an example analysis
  29 // this could be expanded to support more modes (e.g. not writing to file but immediately
  30 // doing the desired v2 analysis)
  31 void GenerateEventsOnTheFly(Bool_t generate = kTRUE) {
  32   if(!generate) {
  33       LoadLibraries();
  34       AnalyseEventsFromFile();
  35       return;
  36   }
  37   TStopwatch aa;
  38   aa.Start();
  39   // load libraries
  40   if(!LoadLibraries()) {
  41       printf(" > problem loading libraries, maybe you need to fix your includes or pythia path <");
  42       exit(0);
  43   }
  44    ///========= GENERATOR INTERFACE ======///
  45   // this is the generator interface. here you can specify what your event should look like in terms of
  46   // mother particles (define the available species and their multiplicities and define a background)
  47   // see the decayer interface for decayer specifics
  48   // in this example macro, pythia is used as a decayer.
  49   // however, any decayer that inherits from TVirualMCDecayer can be used
  50   //
  51   // by default, all tracks are made rp's, so 'offline' tagging is necessary.
  52   // decayed particles are not saved in the output events.
  53   // output is built up from files containing AliFlowEventSimple objects, which in turn consist of
  54   // AliFlowSimpleTracks, in which a few things should be noted:
  55   // 1) Charge() tells whether a particle is a primary (-1) or a secondary (1)
  56   // 2) Weight() stores the longitudinal momentum (Pz)
  57   // 3) GetID() gives the pdg value of the genreated track
  58   // ===================================///
  59   //
  60   // a) select which species you want to have in your initial event (pdg code)
  61   Int_t nSpecies        = 7; // number of poi species
  62   Int_t species[]       = {     22, 333, 221, 311, 3122, 3334, 313}; // gamma, phi eta, k0, lambda, omega, kstar
  63   // b) specify their multiplicities per event
  64   Int_t multiplicities[] = {    10, 10, 10, 10, 10, 10, 10}; // so we get 10 particles of each species here
  65   // c) set the multiplicity for charged background (mixture of pi, k, pi)
  66   Int_t setChargedBackground = 300;
  67   // d) poissonian fluctuations for above multiplicities
  68   Bool_t fluctuate      = kFALSE;
  69   // e) which tracks should get a flow boost ?
  70   Bool_t addMotherV2    = kTRUE;
  71   Bool_t addDaughterV2  = kFALSE;
  72   // f) for each species (poi or rp, or species that will be created by the decayer) you can specify 'custom'
  73   //    functions for pt, v2 or v3. this will override defaults
  74   //    exapmle: modify the v2 for species 310
  75   TString       myFunction   = "0.05*TMath::Sin(x)+.07";
  76   Int_t         modifyV2    = 310;
  77   // and call - once an instance of the class is created (see below)
  78   // AliFlowOnTheFlyEventGenerator::SetPtDependentV2(myFunction.Data(), modifyV2)
  79
  80   // == done with the flow stuff, proceed to the decayer, qa and writing to file == //
  81
  82   // g) decayer modes. select which decays needs to be forced (in the example by pythia)
  83   // note that by default gamma's are 'decayed' to e+ e- pairs
  84   Int_t decayModes[]    = { TPythia6Decayer::kHadronicD };     // some decay modes
  85   // h) number of events you want to generate.
  86   const int events      = 50000;
  87   // i) write output to file
  88   Bool_t writeOutput    = kTRUE;
  89   // j) do qa. qa histograms will be created at the end of generation adn written to a rootfile in your working directory
  90   Bool_t qa             = kTRUE;
  91   // k) write analysis to an output file
  92   TString trunkName             =       "OnTheFlyEvent"; // trunk of output file name
  93   // l) specify the max number of events per file (new files will be generated automatically)
  94   const int maxEventsPerFile    =    1000;          // specify the maximum number of events that is stored per file
  95                                                      // note that events are stored temporarily in RAM memory before writing to file
  96                                                      // so setting this number to a very large value could lead to an unresponsive system
  97   ///====== END OF GENERATOR INTERFACE===///
  98   // from here there's nothing that needs to be changed for the event generation
  99
 100   int nFileCounter(0);
 101   TClonesArray *event = new TClonesArray("TParticle",5000);
 102   TObjArray *dataContainer = new TObjArray(maxEventsPerFile);
 103   dataContainer->SetOwner(kTRUE);
 104   TFile* _tmpfile = 0x0;
 105   // setup the decayer
 106   TPythia6Decayer* decayer = new TPythia6Decayer();
 107   Int_t nDecayModes(sizeof(decayModes)/sizeof(Int_t));
 108   // get an estimate of no of tracks per event
 109   Int_t nMult(setChargedBackground);
 110   for(Int_t i(0); i < nSpecies; i++) nMult += species[i];
 111   for(Int_t i(0); i < nDecayModes; i++) decayer->SetForceDecay(decayModes[i]);
 112   // get an instance of the class that we'll use to generate events
 113   AliFlowOnTheFlyEventGenerator* eventGenerator = new AliFlowOnTheFlyEventGenerator(    qa,             // make some QA plots
 114                                                                                         3,              // introduce flow fluctuations
 115                                                                                         nMult,          // set initial value for evnet multiplcity
 116                                                                                         decayer,        // specify the decayer (NULL is no decay)
 117                                                                                         addMotherV2,    // add v2 for mothers
 118                                                                                         NULL,           // add v3 for mothers  (not implemented yet)
 119                                                                                         addDaughterV2,  // add v2 for daughters
 120                                                                                         NULL);          // add v3 for daughters (not implemented yet)
 121   // example of how to pass a custom function for v2 for a certain species to the generator
 122   //  eventGenerator->SetPtDependentV2(myFunction.Data(), modifyV2);
 123   for(int i(0); i < events; i++) {      // event generator loop
 124      // prepare I/O necessities
 125      if(!_tmpfile) _tmpfile = new TFile(Form("%s_%i.root", trunkName.Data(), nFileCounter), "RECREATE");       // new file, overwrite if already exists
 126      // for each event, one can embed a TClonesArray of type TParticle*, e.g. do
 127      //       eventGenerator->EmbedEvent(myEmbeddedEvent);
 128      // now that the generator is prepared, generate the actual event
 129      AliFlowEventSimple* _tmpevent = eventGenerator->GenerateOnTheFlyEvent(event, nSpecies, species, multiplicities, setChargedBackground, fluctuate);
 130      if(i==0) printf(" \n\n > Generating events on the fly, please be patient ... < \n");
 131      _tmpevent->Write();
 132      // check if we need to open a new file, or if we've reached the end of the generator cycle
 133      if((i-nFileCounter*maxEventsPerFile)==(maxEventsPerFile-1) || i == (events-1)) {
 134          printf("   - writing memory buffer %i to output file ... \n", nFileCounter);
 135          dataContainer->Write();
 136          dataContainer->Clear();
 137          _tmpfile->Close();
 138          _tmpfile->Delete();
 139          _tmpfile = 0x0;
 140          nFileCounter++;
 141      }
 142   }
 143   if(qa) eventGenerator->DoGeneratorQA(kTRUE, kFALSE);
 144   delete eventGenerator;
 145   delete decayer;
 146   delete event;
 147   aa.Stop();
 148   printf("\n > events have been generated and written to file ! < \n");
 149   aa.Print();
 150 }
 151 //_____________________________________________________________________________
 152 void AnalyseEventsFromFile() {
 153     // example function: read events from file and do some analysis.
 154     Int_t nFiles(0);                                    // file counter
 155     // setup analysis methods
 156     AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 157     PrepareQCumulants(qc);
 158
 159     while(kTRUE) {
 160         TFile file(Form("OnTheFlyEvent_%i.root",nFiles));         // open the root file
 161         if(file.IsZombie()) break;
 162         TIter iter(file.GetListOfKeys());           // get a list of keys
 163
 164         // loop over all the events in the file
 165         while(kTRUE) {                              // infinite loop ...
 166             TKey* key = iter();                     // get the next key from the file
 167             if(!key) break;                         // ... and exit the loop if the key is empty
 168             AliFlowEventSimple* event = (AliFlowEventSimple*)key->ReadObj();        // cast to a flow event
 169             printf("   - read event with %i tracks\n", event->NumberOfTracks());
 170             // track loop
 171             for(int i(0); i < event->NumberOfTracks(); i++) {
 172                 // do some offline fun stuff with the events
 173                 AliFlowTrackSimple* track = event->GetTrack(i);
 174                 if(!track) continue;                // break the loop when we don't find a track
 175                                                     // do some next stuff with the tracks
 176                 // example: retag the primary kaons as poi's and the rest as reference
 177                 if(track->Charge()!=-1) continue;           // reject secondary tracks UGLY !!!
 178                 if(TMath::Abs(track->GetID()==321)) {       // find the kaons
 179                     if(track->InRPSelection()) {            // if the kaon is an rp (it should be by default)
 180                         track->SetForPOISelection(kTRUE);   // set it as poi ...
 181                         track->SetForRPSelection(kFALSE);   // ... and 'detag' it as rp
 182                         // make sure to correct the reference multiplicity
 183                         event->SetReferenceMultiplicity(-1+event->GetReferenceMultiplicity());
 184                     }
 185                     else track->SetForPOISelection(kTRUE);  // tag as poi
 186                 }
 187                 else {                                      // if track is not a kaon
 188                     if(!track->InRPSelection()) {           // check if it's an rp
 189                         track->SetForRPSelection(kTRUE);    // if not, tag as rp (it should be there, but just to be safe ... )
 190                         // and again correct the reference multiplicity
 191                         event->SetReferenceMultiplicity(1+event->GetReferenceMultiplicity());
 192                     }
 193                 }
 194             }
 195             // insert the poor guy into the flow anaysis method
 196             qc->Make(event);
 197             delete event;
 198         }
 199         nFiles++;
 200     }
 201     // prepare output for the flow package analyses
 202     TString outputFileName = "AnalysisResults.root";
 203     TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
 204     const Int_t nMethods(1);
 205     TString method[] = {"QC"};
 206     TDirectoryFile *dirFileFinal[nMethods] = {NULL};
 207     TString fileName[nMethods];
 208     for(Int_t i(0); i < nMethods; i++) {
 209          fileName[i]+="output";
 210          fileName[i]+=method[i].Data();
 211          fileName[i]+="analysis";
 212          dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
 213     }
 214     qc->Finish();qc->WriteHistograms(dirFileFinal[0]);
 215     outputFile->Close();
 216     delete outputFile;
 217 }
 218
 219 /* ====================================================
 220  *
 221  * the following two functions serve as an exmaple
 222  * of how to run the phi and kstar analysis on on the
 223  * fly events
 224  *
 225  * ====================================================
 226  */
 227
 228 //_____________________________________________________________________________
 229 void TwoParticleResonanceFlowOnTheFly() {
 230
 231     // load additional libraries
 232     gSystem->Load("libGeom");
 233     gSystem->Load("libVMC");
 234     gSystem->Load("libXMLIO");
 235     gSystem->Load("libPhysics");
 236     gSystem->Load("libANALYSIS");
 237     gSystem->Load("libANALYSISalice");
 238     gSystem->Load("libPWGflowBase");
 239     gSystem->Load("libPWGflowTasks");
 240
 241     TString dirName = "reconstruction";
 242     Int_t nFiles(0);
 243     // define the poi cuts, see AddTwoParticelResonanceFlowTask.C for more info
 244     const Int_t mb(30); // no of massbands available for analysis
 245     Float_t minMass(.99), maxMass(1.092);       // upper and lower bound
 246     AliFlowTrackSimpleCuts* POIfilterQC[mb];    // pointers to poi filters
 247     Double_t flowBands[2][mb];                  // define the invm regins
 248     Double_t _inc = (maxMass-minMass)/(float)mb;
 249     for (Int_t _mb = 0; _mb < mb; _mb++) {      // create the poi filters
 250        flowBands[0][_mb] = minMass + _mb * _inc;
 251        flowBands[1][_mb] = minMass + (_mb + 1) * _inc;
 252        POIfilterQC[_mb] = new AliFlowTrackSimpleCuts(Form("FilterPOIQC_MB%d", _mb));
 253        POIfilterQC[_mb]->SetEtaMin(-0.8);       // eta range
 254        POIfilterQC[_mb]->SetEtaMax(+0.8);
 255        POIfilterQC[_mb]->SetMassMin(flowBands[0][_mb]); // invm range lower bound
 256        POIfilterQC[_mb]->SetMassMax(flowBands[1][_mb]); // invm rnage upper bound
 257     }
 258     // do the flow analysis
 259     AliFlowAnalysisWithQCumulants* qc[mb];
 260     for(int i(0); i < mb; i++) {        // init the q-cumulant tasks, one for each invm bin
 261         qc[i] = new AliFlowAnalysisWithQCumulants();
 262         PrepareQCumulants(qc[i]);
 263         printf("  > init qc task %i < \n", i);
 264     }
 265     AliAnalysisTwoParticleResonanceFlowTask* task = new AliAnalysisTwoParticleResonanceFlowTask("onthefly");
 266     SetUpTask(task, minMass,maxMass);// setup the task
 267     // open files
 268     while(kTRUE) {      // infinite loop which we will break when we've looked through all the files
 269         TFile file(Form("OnTheFlyEvent_%i.root", nFiles));   // open the root file
 270         if(nFiles==0&&file.IsZombie()) {                     // something went wrong ...
 271             printf(" > cannot get input files ... exiting < \n");
 272             exit(0);
 273         }
 274         if(file.IsZombie()) break;                  // break the loop on the first empty file
 275         TIter iter(file.GetListOfKeys());           // get a list of keys
 276         while(kTRUE) {                              // infinite loop over events ...
 277             TKey* key = iter();                     // get the next key from the file
 278             if(!key) break;                         // ... and exit the loop if the key is empty
 279             AliFlowEventSimple* event = (AliFlowEventSimple*)key->ReadObj();        // cast to a flow event
 280             printf("   - read event with %i tracks from file %i \n       > task ", event->NumberOfTracks(), nFiles);
 281             task->DoAnalysisOnTheFly(event);                       // do the on the fly analysis
 282             AliFlowEventSimple* flowEvent = task->GetFlowEvent();  // retrieve the flow event
 283             for(int j(0); j < mb; j++) {                           // loop over all invm bands
 284                 flowEvent->TagPOI(POIfilterQC[j]);                 // 'offline' tagging of poi's in certain mass range
 285                 flowEvent->TagSubeventsInEta(-.8, 0, 0, .8);       // setup subevents
 286                 qc[j]->Make(flowEvent);                            // do qc analysis
 287                 printf(" %i", j);
 288             }
 289         }
 290         nFiles++;
 291     }
 292     // prepare output for the flow package analyses
 293     TFile *outputFile = new TFile("AnalysisResults.root","RECREATE");   // common outputfile
 294     TDirectoryFile* dirFileFinal[mb];                                   // tdirectory files
 295     TString fileName[mb];                                               // dir names
 296     for(Int_t i(0); i < mb; i++) {                                      // loop over all bands ...
 297          fileName[i]+=Form("QC_minv_%i", i);
 298          dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
 299          qc[i]->Finish();                               // finalize the method
 300          qc[i]->WriteHistograms(dirFileFinal[i]);       // and write it to file
 301     }
 302     // write the output of the phi reconstruction to the same file
 303     TDirectoryFile* dir = new TDirectoryFile(dirName.Data(), dirName.Data());
 304     task->DoAnalysisOnTheFly(dir);
 305     // end of analysis
 306     outputFile->Close();
 307     delete outputFile;
 308 }
 309 //_____________________________________________________________________________
 310 void SetUpTask(AliAnalysisTwoParticleResonanceFlowTask* task, Float_t minMass, Float_t maxMass)
 311 {
 312     // some magic which is necessary to 'trick' the analysis task into thinking
 313     // that we're actually doing an analysis on aod events.
 314     // some of these configurations are used (e.g. setupSpecies, common constants
 315     // and the binning in pt)
 316     // others have no meaning (PID, DCA, RP and POI cuts)
 317     // note that UserCreateOutputObjects is necessary since it initializes
 318     // the output of the analysis
 319     gSystem->Load("libPWGflowTasks");
 320     Float_t PIDconfig[] = {0,0,0,0,0,0,0.3};                    // not used
 321     task->SetPIDConfiguration(PIDconfig);                       // not used
 322     task->SetupSpeciesA(321, 1, 4.93676999e-01, 0.15, 5.);      // pid and charge
 323     task->SetupSpeciesB(-321, -1, 4.93676999e-01, 0.15, 5.);    // pid and charge
 324     task->SetRPCuts(new AliFlowTrackCuts("unused_rp_cuts"));    // not used
 325     task->SetPOICuts(new AliFlowTrackCuts("unused_poi_cuts"));  // not used
 326     Float_t POIDCA[] = {0,0,0,0,0};                             // not used
 327     task->SetPOIDCAXYZ(POIDCA);                                 // not used
 328     task->SetCommonConstants(30, minMass, maxMass);             // necesssary
 329     Float_t _pt[] = {0., 0.6, 1.2, 1.8, 2.4, 3., 4., 5., 6., 7.};       // same
 330     task->SetPtBins(_pt, (Int_t)(sizeof(_pt)/sizeof(_pt[1]))-1);        // same
 331     task->UserCreateOutputObjects();                            // necessary
 332 }
 333 //_____________________________________________________________________________
 334 void PrepareQCumulants(AliFlowAnalysisWithQCumulants* qc)
 335 {
 336     // prepare the cumulants task
 337     qc->SetHarmonic(2);
 338     qc->SetCalculateDiffFlow(kTRUE);
 339     qc->SetCalculate2DDiffFlow(kFALSE); // vs (pt,eta)
 340     qc->SetApplyCorrectionForNUA(kFALSE);
 341     qc->SetFillMultipleControlHistograms(kFALSE);
 342     qc->SetMultiplicityWeight("combinations"); // default (other supported options are "unit" and "multiplicity")
 343     qc->SetCalculateCumulantsVsM(kFALSE);
 344     qc->SetCalculateAllCorrelationsVsM(kFALSE); // calculate all correlations in mixed harmonics "vs M"
 345     qc->SetnBinsMult(10000);
 346     qc->SetMinMult(0);
 347     qc->SetMaxMult(10000);
 348     qc->SetBookOnlyBasicCCH(kFALSE); // book only basic common control histograms
 349     qc->SetCalculateDiffFlowVsEta(kTRUE); // if you set kFALSE only differential flow vs pt is calculated
 350     qc->Init();
 351 }
 352 //_____________________________________________________________________________
 353 Bool_t LoadLibraries() {
 354     // load libraries
 355     gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS  -I$ALICE_ROOT/OADB -I$ALICE_ROOT/PWGHF -I$ALICE_ROOT/PWGHF/base -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWG/FLOW/Base -I$ALICE_ROOT/PWG/FLOW/Tasks -g");
 356     gSystem->Load("libGeom");
 357     gSystem->Load("libVMC");
 358     gSystem->Load("libXMLIO");
 359     gSystem->Load("libPhysics");
 360     gSystem->Load("libEG");
 361     gSystem->Load("libPWGflowBase");
 362     if(gSystem->Load("$ALICE_ROOT/lib/tgt_linuxx8664gcc/libpythia6")!=0) {
 363         printf(" \n\n\n *** fatal error, couldn't load pythia, exiting ! *** \n\n\n");
 364         return kFALSE;
 365     }
 366     return kTRUE;
 367 }
 368 //_____________________________________________________________________________