PWG2/FLOW/other/runFlowAnalysis.C

   1 #include "TStopwatch.h"
   2 #include "TObjArray"
   3 #include "Riostream.h"
   4
   5 //RUN SETTINGS
   6 //flow analysis method can be: (set to kTRUE or kFALSE)
   7 Bool_t SP    = kTRUE;
   8 Bool_t LYZ1  = kTRUE;
   9 Bool_t LYZ2  = kFALSE;
  10 Bool_t LYZEP = kFALSE;
  11 Bool_t GFC   = kTRUE;
  12 Bool_t QC    = kTRUE;
  13 Bool_t FQD   = kTRUE;
  14 Bool_t MCEP  = kFALSE; //does not work yet 24/12/08
  15
  16 Int_t offset = 0;
  17
  18 //CUT SETTINGS
  19 //integrated selection
  20 Double_t ptMaxInt  = 10.;
  21 Double_t ptMinInt  = 0.;
  22 Double_t etaMaxInt = 1.;
  23 Double_t etaMinInt = -1.;
  24 Double_t phiMaxInt = 7.5;
  25 Double_t phiMinInt = 0.;
  26 Int_t PIDInt       = 211;
  27
  28 //differential selection
  29 Double_t ptMaxDiff  = 10.;
  30 Double_t ptMinDiff  = 0.;
  31 Double_t etaMaxDiff = 1.;
  32 Double_t etaMinDiff = -1.;
  33 Double_t phiMaxDiff = 7.5;
  34 Double_t phiMinDiff = 0.;
  35 Int_t PIDDiff       = 211;
  36
  37 int runFlowAnalysis(Int_t aRuns = 100, const char*
  38                           //                      dir="/data/alice1/kolk/KineOnly3/")
  39                           dir="/Users/snelling/alice_data/KineOnly3/")
  40 {
  41   TStopwatch timer;
  42   timer.Start();
  43
  44   if (LYZ1 && LYZ2) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(); }
  45
  46   if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  47
  48   if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  49
  50
  51   cout<<endl;
  52   cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
  53   cout<<endl;
  54
  55   //  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  56   gSystem->AddIncludePath("-I$ROOTSYS/include");
  57
  58   // constants
  59   gROOT->LoadMacro("code/AliFlowCommonConstants.cxx+");
  60   gROOT->LoadMacro("code/AliFlowLYZConstants.cxx+");
  61   gROOT->LoadMacro("code/AliFlowCumuConstants.cxx+");
  62
  63   // flow event
  64   gROOT->LoadMacro("code/AliFlowVector.cxx+");
  65   gROOT->LoadMacro("code/AliFlowTrackSimple.cxx+");
  66   gROOT->LoadMacro("code/AliFlowEventSimple.cxx+");
  67
  68   // cuts
  69   gROOT->LoadMacro("code/AliFlowTrackSimpleCuts.cxx+");
  70
  71   // output histosgrams
  72   gROOT->LoadMacro("code/AliFlowCommonHist.cxx+");
  73   gROOT->LoadMacro("code/AliFlowCommonHistResults.cxx+");
  74   gROOT->LoadMacro("code/AliFlowLYZHist1.cxx+");
  75   gROOT->LoadMacro("code/AliFlowLYZHist2.cxx+");
  76
  77   // functions needed for various methods
  78   gROOT->LoadMacro("code/AliCumulantsFunctions.cxx+");
  79   gROOT->LoadMacro("code/AliQCumulantsFunctions.cxx+");
  80   gROOT->LoadMacro("code/AliFittingFunctionsForQDistribution.cxx+");
  81   gROOT->LoadMacro("code/AliFlowLYZEventPlane.cxx+");
  82
  83   // Flow Analysis code for various methods
  84   gROOT->LoadMacro("code/AliFlowAnalysisWithMCEventPlane.cxx+");
  85   gROOT->LoadMacro("code/AliFlowAnalysisWithScalarProduct.cxx+");
  86   gROOT->LoadMacro("code/AliFlowAnalysisWithLYZEventPlane.cxx+");
  87   gROOT->LoadMacro("code/AliFlowAnalysisWithLeeYangZeros.cxx+");
  88   gROOT->LoadMacro("code/AliFlowAnalysisWithCumulants.cxx+");
  89   gROOT->LoadMacro("code/AliFlowAnalysisWithQCumulants.cxx+");
  90   gROOT->LoadMacro("code/AliFittingQDistribution.cxx+");
  91
  92   // Method to fill the FlowEvent
  93   gROOT->LoadMacro("code/FlowEventSimpleMaker.cxx+");
  94
  95   //load needed libraries
  96   gSystem->Load("libTree.so");
  97
  98   cout << "finished loading macros!" << endl;
  99   //  gSystem->Load("libANALYSIS.so");
 100   //  gSystem->Load("libPWG2flow.so");
 101
 102   //------------------------------------------------------------------------
 103   //cuts
 104   AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
 105   cutsInt->SetPtMax(ptMaxInt);
 106   cutsInt->SetPtMin(ptMinInt);
 107   cutsInt->SetEtaMax(etaMaxInt);
 108   cutsInt->SetEtaMin(etaMinInt);
 109   cutsInt->SetPhiMax(phiMaxInt);
 110   cutsInt->SetPhiMin(phiMinInt);
 111   cutsInt->SetPID(PIDInt);
 112
 113   AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
 114   cutsDiff->SetPtMax(ptMaxDiff);
 115   cutsDiff->SetPtMin(ptMinDiff);
 116   cutsDiff->SetEtaMax(etaMaxDiff);
 117   cutsDiff->SetEtaMin(etaMinDiff);
 118   cutsDiff->SetPhiMax(phiMaxDiff);
 119   cutsDiff->SetPhiMin(phiMinDiff);
 120   cutsDiff->SetPID(PIDDiff);
 121
 122   //------------------------------------------------------------------------
 123   //flow event
 124   FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker();
 125
 126   AliFlowAnalysisWithQCumulants    *qc    = NULL;
 127   AliFlowAnalysisWithCumulants     *gfc   = NULL;
 128   AliFittingQDistribution          *fqd   = NULL;
 129   AliFlowAnalysisWithLeeYangZeros  *lyz1  = NULL;
 130   AliFlowAnalysisWithLeeYangZeros  *lyz2  = NULL;
 131   AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
 132   AliFlowAnalysisWithScalarProduct *sp    = NULL;
 133   AliFlowAnalysisWithMCEventPlane  *mcep  = NULL;
 134
 135   //flow methods:
 136   //MCEP = monte carlo event plane
 137   AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 138   if (MCEP)
 139     {
 140       mcep->Init();
 141     }
 142
 143   //QC = Q-cumulants
 144   if(QC)
 145     {
 146       AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 147       qc->CreateOutputObjects();
 148     }
 149
 150   //GFC = Generating Function Cumulants
 151   if(GFC)
 152     {
 153       AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 154       gfc->CreateOutputObjects();
 155     }
 156
 157   //FQD = Fitting q-distribution
 158   AliFittingQDistribution* fqd = new AliFittingQDistribution();
 159   if(FQD)
 160     {
 161       fqd->CreateOutputObjects();
 162     }
 163
 164   //LYZ1 = Lee-Yang Zeroes first run
 165   AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
 166   if(LYZ1)
 167     {
 168       lyz1->SetFirstRun(kTRUE);
 169       lyz1->SetUseSum(kTRUE);
 170       lyz1->Init();
 171     }
 172
 173   //LYZ2 = Lee-Yang Zeroes second run
 174   AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
 175   if(LYZ2)
 176     {
 177       // read the input file from the first run
 178       TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
 179       TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
 180       if(!inputFileLYZ2 || inputFileLYZ2->IsZombie()) {
 181         cerr << " ERROR: NO First Run file... " << endl ; }
 182       else {
 183         TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");
 184         if (!inputListLYZ2) {cout<<"list is NULL pointer!"<<endl;
 185         }
 186         else {
 187           cout<<"LYZ2 input file/list read..."<<endl;
 188           lyz2->SetFirstRunList(inputListLYZ2);
 189           lyz2->SetFirstRun(kFALSE);
 190           lyz2->SetUseSum(kTRUE);
 191           lyz2->Init();
 192         }
 193       }
 194     }
 195
 196  //LYZEP = Lee-Yang Zeroes event plane
 197   AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 198   AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 199   if(LYZEP)
 200     {
 201       // read the input file from the second lyz run
 202       TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
 203       TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 204       if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 205         cerr << " ERROR: NO Second Run file... " << endl ; }
 206       else {
 207         TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");
 208         if (!inputListLYZEP) {cout<<"list is NULL pointer!"<<endl;
 209         }
 210         else {
 211           cout<<"LYZEP input file/list read..."<<endl;
 212           ep   ->SetSecondRunList(inputListLYZEP);
 213           lyzep->SetSecondRunList(inputListLYZEP);
 214           ep   ->Init();
 215           lyzep->Init();
 216         }
 217       }
 218     }
 219
 220   //SP = Scalar Product
 221   AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 222   if(SP)
 223     {
 224       sp->Init();
 225     }
 226
 227
 228
 229   //------------------------------------------------------------------------
 230
 231   //standard code
 232   Int_t fCount = 0;
 233   TString execDir(gSystem->pwd());
 234   TString targetDir(dir);
 235   TSystemDirectory* baseDir = new TSystemDirectory(".", dir);
 236   TList* dirList           = baseDir->GetListOfFiles();
 237   if (!dirList) {
 238     cout << endl << "No input files in: " << targetDir.Data() << endl;
 239     break;
 240   }
 241   Int_t nDirs              = dirList->GetEntries();
 242   cout<<endl;
 243   cout<<"Int_t nDirs = "<<nDirs<<endl;
 244   gSystem->cd(execDir);
 245
 246   for(Int_t iDir=0;iDir<nDirs;++iDir)
 247     {
 248       TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
 249       if(!presentDir || !presentDir->IsDirectory() ||
 250          strcmp(presentDir->GetName(), ".") == 0 ||
 251          strcmp(presentDir->GetName(), "..") == 0)
 252         {
 253           cout << endl;
 254           cout << "Directory (" << iDir << "): " << presentDir->GetName() <<
 255             " - Skipping ... " << endl;
 256           continue ;
 257         }
 258
 259       if(offset > 0) { --offset ; continue ; }
 260       if((aRuns > 0) && (fCount >= aRuns)) { break ; }
 261
 262       TString presentDirName(dir); // aDataDir
 263       presentDirName += presentDir->GetName();
 264       presentDirName += "/";
 265       //cerr<<" presentDirName = "<<presentDirName<<endl;
 266
 267       TString fileName = presentDirName;
 268       fileName += "galice.root";
 269       Long_t *id, *size, *flags, *modtime;
 270       if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime))
 271         {
 272           cout << " File : " << fileName << " does NOT exist ! - Skipping ... "
 273                << endl;
 274           continue;
 275         }
 276       cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName
 277            << " ... " << endl;
 278
 279       //loop (simulations in the present dir)
 280       TSystemDirectory* evtsDir = new TSystemDirectory(".",
 281                                                        presentDirName.Data());
 282       TList* fileList       = evtsDir->GetListOfFiles();
 283       Int_t nFiles                  = fileList->GetEntries();
 284       cout<<" Int_t nFiles = "<<nFiles<<endl;
 285       gSystem->cd(execDir);
 286       for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
 287         {
 288           TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
 289           TString presentFileName(presentDirName);
 290           presentFileName += presentFile->GetName();
 291
 292           if(!(presentFileName.Contains("Kinematics") &&
 293                presentFileName.Contains("root"))) { continue ; }
 294
 295           cout << " found: " << presentFileName.Data() << endl;
 296
 297           TFile* kineFile = new TFile(presentFileName.Data(), "READ");
 298           // kineFile->ls();
 299           Int_t nEvts = kineFile->GetNkeys() ;
 300           cout << "  . found: " << nEvts << " KineTree(s) in " <<
 301             presentFileName.Data() << endl;
 302           TList* kineEventsList = (TList*)kineFile->GetListOfKeys();
 303           TTree* kTree;
 304           TIter next(kineEventsList);
 305           TKey* key;
 306
 307           //loop over the events
 308           while( key=(TKey *)next() )
 309             {
 310               TDirectory* tDir = (TDirectory*)key->ReadObj();
 311               if(!tDir) break;
 312
 313               TString evtDir(tDir->GetName());
 314               cout << "  . . found: " << tDir->GetName() << endl;
 315
 316               kTree = (TTree *)tDir->Get("TreeK");
 317               if(!kTree) break;
 318
 319               Int_t nPart = kTree->GetEntries();
 320               cout << "  . . . kTree " << fCount << " has " << nPart <<
 321                 " particles " << endl;
 322
 323
 324
 325               //-----------------------------------------------------------
 326               //fill and save the flow event
 327
 328               AliFlowEventSimple* fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff);
 329
 330               //pass the flow event to flow methods for analysis
 331               //MCEP
 332               if (MCEP)
 333                 {
 334                   //mcep->Make(fEvent,fEP);  //fix fEP
 335                   cout<<"  --> MCEP analysis..."<<endl;
 336                 }
 337               //QC
 338               if(QC)
 339                 {
 340                   qc->Make(fEvent);
 341                   cout<<"  --> QC analysis..."<<endl;
 342                 }
 343               //GFC
 344               if(GFC)
 345                 {
 346                   gfc->Make(fEvent);
 347                   cout<<"  --> GFC analysis..."<<endl;
 348                 }
 349               //FQD
 350               if(FQD)
 351                 {
 352                   fqd->Make(fEvent);
 353                   cout<<"  --> FQD analysis..."<<endl;
 354                 }
 355               //LYZ1
 356               if(LYZ1)
 357                 {
 358                   lyz1->Make(fEvent);
 359                   cout<<"  --> LYZ1 analysis..."<<endl;
 360                 }
 361               //LYZ2
 362               if(LYZ2)
 363                 {
 364                   lyz2->Make(fEvent);
 365                   cout<<"  --> LYZ2 analysis..."<<endl;
 366                 }
 367               //LYZEP
 368               if(LYZEP)
 369                 {
 370                   lyzep->Make(fEvent,ep);
 371                   cout<<"  --> LYZEP analysis..."<<endl;
 372                 }
 373               //SP
 374               if(SP)
 375                 {
 376                   sp->Make(fEvent);
 377                   cout<<"  --> SP analysis..."<<endl;
 378                 }
 379
 380
 381
 382               //-----------------------------------------------------------
 383
 384               fCount++;
 385               delete kTree;
 386               delete fEvent;
 387             }
 388           delete kineFile ;
 389         }
 390       delete evtsDir ;
 391     }
 392
 393   //--------------------------------------------------------------
 394   //calculating and storing the final results of flow analysis
 395   //MCEP
 396   if (MCEP)
 397     {
 398       mcep->Finish();
 399       TString *outputFileNameMCEP = new TString("outputMCEPanalysis.root");
 400       //mcep->WriteHistograms(outputFileNameMCEP); //add this method to MCEP classes
 401       delete outputFileNameMCEP;
 402     }
 403   //QC
 404   if(QC)
 405     {
 406       qc->Finish();
 407       TString *outputFileNameQC = new TString("outputQCanalysis.root");
 408       qc->WriteHistograms(outputFileNameQC);
 409       delete outputFileNameQC;
 410     }
 411   //GFC
 412   if(GFC)
 413     {
 414       gfc->Finish();
 415       TString *outputFileNameGFC = new TString("outputGFCanalysis.root");
 416       gfc->WriteHistograms(outputFileNameGFC);
 417       delete outputFileNameGFC;
 418     }
 419   //FQD
 420   if(FQD)
 421     {
 422       fqd->Finish();
 423       TString *outputFileNameFQD = new TString("outputFQDanalysis.root");
 424       fqd->WriteHistograms(outputFileNameFQD);
 425       delete outputFileNameFQD;
 426     }
 427   //LYZ1
 428   if(LYZ1)
 429     {
 430       lyz1->Finish();
 431       TString *outputFileNameLYZ1 = new TString("outputLYZ1analysis.root");
 432       lyz1->WriteHistograms(outputFileNameLYZ1);
 433       delete outputFileNameLYZ1;
 434     }
 435   //LYZ2
 436   if(LYZ2)
 437     {
 438       lyz2->Finish();
 439       TString *outputFileNameLYZ2 = new TString("outputLYZ2analysis.root");
 440       lyz2->WriteHistograms(outputFileNameLYZ2);
 441       delete outputFileNameLYZ2;
 442     }
 443   //LYZEP
 444   if(LYZEP)
 445     {
 446       lyzep->Finish();
 447       TString *outputFileNameLYZEP = new TString("outputLYZEPanalysis.root");
 448       lyzep->WriteHistograms(outputFileNameLYZEP);
 449       delete outputFileNameLYZEP;
 450     }
 451   //SP
 452   if(SP)
 453     {
 454       sp->Finish();
 455       TString *outputFileNameSP = new TString("outputSPanalysis.root");
 456       sp->WriteHistograms(outputFileNameSP);
 457       delete outputFileNameSP;
 458     }
 459
 460
 461
 462   //--------------------------------------------------------------
 463
 464   cout << endl;
 465   cout << " Finished ... " << endl;
 466   cout << endl;
 467
 468   timer.Stop();
 469   cout << endl;
 470   timer.Print();
 471 }