simulate partial acceptance
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 13 Jun 2009 11:21:46 +0000 (11:21 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 13 Jun 2009 11:21:46 +0000 (11:21 +0000)
PWG2/FLOW/AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx
PWG2/FLOW/AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.h
PWG2/FLOW/macros/runFlowAnalysisOnTheFly.C

index 84d3edf..08f1f1c 100644 (file)
@@ -56,7 +56,13 @@ AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t iseed):
   fV4RP(0.), 
   fV4SpreadRP(0.), 
   fV2RPMax(0.), 
-  fPtCutOff(0.), 
+  fPtCutOff(0.),
+  fPhiMin1(0.),              
+  fPhiMax1(0.),             
+  fProbability1(0.),       
+  fPhiMin2(0.),   
+  fPhiMax2(0.),            
+  fProbability2(0.),          
   fPtSpectra(NULL),
   fPhiDistribution(NULL),
   fMyTRandom3(NULL),
@@ -182,6 +188,12 @@ AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly()
   Double_t dTmpEta = 0.;
   Double_t dTmpV2 = 0.;
   Double_t dTmpPhi = 0.;
+  Bool_t bUniformAcceptance = kTRUE;
+  Double_t Pi = TMath::Pi();
+  if(!((fPhiMin1==0.) && (fPhiMax1==0.) && (fPhiMin2==0.) && (fPhiMax2==0.)))
+  {
+   bUniformAcceptance = kFALSE;
+  }
   for(Int_t i=0;i<iNewMultiplicityOfRP;i++) {
     dTmpPt = fPtSpectra->GetRandom();
     dTmpEta = fMyTRandom3->Uniform(dEtaMin,dEtaMax);
@@ -198,15 +210,63 @@ AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly()
     // add the track to the event
     for(Int_t d=0;d<fNoOfLoops;d++) {
       AliFlowTrackSimple* pTrack = new AliFlowTrackSimple();
-      pTrack->SetPt(dTmpPt);
-      pTrack->SetEta(dTmpEta);
-      pTrack->SetPhi(dTmpPhi);
-      pTrack->SetForRPSelection(kTRUE);
-      iSelParticlesRP++;
-      pTrack->SetForPOISelection(kTRUE);
-      iSelParticlesPOI++;
-      pEvent->TrackCollection()->Add(pTrack);
-      iGoodTracks++;
+      
+      // uniform acceptance:
+      if(bUniformAcceptance)
+      {
+       pTrack->SetPt(dTmpPt); 
+       pTrack->SetEta(dTmpEta); 
+       pTrack->SetPhi(dTmpPhi); 
+       pTrack->SetForRPSelection(kTRUE); 
+       iSelParticlesRP++; 
+       pTrack->SetForPOISelection(kTRUE); 
+       iSelParticlesPOI++; 
+       pEvent->TrackCollection()->Add(pTrack); 
+       iGoodTracks++; 
+      }
+      // non-uniform acceptance, 1st sector:
+      else if ((dTmpPhi > fPhiMin1*Pi/180) && (dTmpPhi < fPhiMax1*Pi/180))
+      {
+       if(fMyTRandom3->Uniform(0,1) > 1 - fProbability1)
+       {
+        pTrack->SetPt(dTmpPt);
+        pTrack->SetEta(dTmpEta);
+        pTrack->SetPhi(dTmpPhi);
+        pTrack->SetForRPSelection(kTRUE);
+        iSelParticlesRP++;
+        pTrack->SetForPOISelection(kTRUE);
+        iSelParticlesPOI++;
+        pEvent->TrackCollection()->Add(pTrack);
+        iGoodTracks++;
+       }
+      } 
+      // non-uniform acceptance, 2nd sector:
+      else if ((dTmpPhi > fPhiMin2*Pi/180) && (dTmpPhi < fPhiMax2*Pi/180))
+      {
+       if(fMyTRandom3->Uniform(0,1) > 1 - fProbability2)
+       {
+        pTrack->SetPt(dTmpPt);
+        pTrack->SetEta(dTmpEta);
+        pTrack->SetPhi(dTmpPhi);
+        pTrack->SetForRPSelection(kTRUE);
+        iSelParticlesRP++;
+        pTrack->SetForPOISelection(kTRUE);
+        iSelParticlesPOI++;
+        pEvent->TrackCollection()->Add(pTrack);
+        iGoodTracks++;
+       }
+      } else 
+        {
+         pTrack->SetPt(dTmpPt);
+         pTrack->SetEta(dTmpEta);
+         pTrack->SetPhi(dTmpPhi);
+         pTrack->SetForRPSelection(kTRUE);
+         iSelParticlesRP++;
+         pTrack->SetForPOISelection(kTRUE);
+         iSelParticlesPOI++;
+         pEvent->TrackCollection()->Add(pTrack);
+         iGoodTracks++;
+        }
     }
   }
   
index 78121a4..9ae3aa9 100644 (file)
@@ -94,6 +94,24 @@ class AliFlowEventSimpleMakerOnTheFly {
   
   void SetPtCutOff(Double_t dPtCutOff) {this->fPtCutOff = dPtCutOff;}
   Double_t GetPtCutOff() const {return this->fPtCutOff;} 
+  
+  void SetFirstSectorPhiMin(Double_t dPhiMin1) {this->fPhiMin1 = dPhiMin1;}
+  Double_t GetFirstSectorPhiMin() const {return this->fPhiMin1;} 
+  
+  void SetFirstSectorPhiMax(Double_t dPhiMax1) {this->fPhiMax1 = dPhiMax1;}
+  Double_t GetFirstSectorPhiMax() const {return this->fPhiMax1;}
+  
+  void SetFirstSectorProbability(Double_t dProbability1) {this->fProbability1 = dProbability1;}
+  Double_t GetFirstProbability() const {return this->fProbability1;}  
+  
+  void SetSecondSectorPhiMin(Double_t dPhiMin2) {this->fPhiMin2 = dPhiMin2;}
+  Double_t GetSecondSectorPhiMin() const {return this->fPhiMin2;} 
+  
+  void SetSecondSectorPhiMax(Double_t dPhiMax2) {this->fPhiMax2 = dPhiMax2;}
+  Double_t GetSecondSectorPhiMax() const {return this->fPhiMax2;}
+  
+  void SetSecondSectorProbability(Double_t dProbability2) {this->fProbability2 = dProbability2;}
+  Double_t GetSecondProbability() const {return this->fProbability2;}  
   //................................................................................................
   
   void SetNoOfLoops(Int_t noofl) {this->fNoOfLoops = noofl;}
@@ -134,6 +152,13 @@ class AliFlowEventSimpleMakerOnTheFly {
   // (pt,eta) dependent harmonics:
   Double_t  fV2RPMax;                // elliptic flow of RPs
   Double_t  fPtCutOff;               // elliptic flow spread of RPs
+  // non-uniform acceptance:
+  Double_t  fPhiMin1;                // first non-uniform sector starts at azimuth fPhiMin1
+  Double_t  fPhiMax1;                // first non-uniform sector ends at azimuth fPhiMax1
+  Double_t  fProbability1;           // particles emitted in fPhiMin1 < phi < fPhiMax1 are taken with probability fProbability1 
+  Double_t  fPhiMin2;                // second non-uniform sector starts at azimuth fPhiMin2
+  Double_t  fPhiMax2;                // second non-uniform sector starts at azimuth fPhiMax2
+  Double_t  fProbability2;           // particles emitted in fPhiMin2 < phi < fPhiMax2 are taken with probability fProbability2
   //................................................................................................
   
   //................................................................................................
index eb9051c..d1964f3 100644 (file)
@@ -14,8 +14,8 @@ Bool_t LYZ2PROD = kFALSE;
 Bool_t LYZEP    = kFALSE;
 Bool_t GFC      = kTRUE;
 Bool_t QC       = kTRUE;
-Bool_t FQD      = kFALSE;
-Bool_t MCEP     = kFALSE;
+Bool_t FQD      = kTRUE;
+Bool_t MCEP     = kTRUE;
 //--------------------------------------------------------------------------------------
 
 // Weights 
@@ -25,7 +25,7 @@ Bool_t usePtWeights  = kFALSE; // pt weights
 Bool_t useEtaWeights = kFALSE; // eta weights
 
 // Parameters for the simulation of events 'on the fly': 
-Bool_t bSameSeed = kTRUE; // use always the same seed for random generators. 
+Bool_t bSameSeed = kFALSE; // use always the same seed for random generators. 
                            // usage of same seed (kTRUE) is relevant in two cases:
                            // 1.) If you want to use LYZ method to calcualte differential flow;
                            // 2.) If you want to use phi weights for GFC, QC and FQD
@@ -34,18 +34,26 @@ Bool_t bConstantHarmonics = kTRUE; // harmonics V1, V2, V4... are constant (kTRU
 
 Int_t iLoops = 1; // number of times to use each track (to simulate nonflow)
 
-Bool_t bMultDistrOfRPsIsGauss = kFALSE; // 1.) if kTRUE  = multiplicitiy of RPs is sampled e-b-e from Gaussian distribution with
+Bool_t bMultDistrOfRPsIsGauss = kTRUE; // 1.) if kTRUE  = multiplicitiy of RPs is sampled e-b-e from Gaussian distribution with
                                         //                 mean = iMultiplicityOfRP and spread = dMultiplicitySpreadOfRP
                                         // 2.) if kFALSE = multiplicitiy of RPs is sampled e-b-e uniformly from 
                                         //                 interval [iMinMultOfRP,iMaxMultOfRP]
                                         // 3.) for a fixed multiplicity use Gaussian with zero spread or use uniform with iMinMult=iMaxMult
                                         
-Bool_t bV2DistrOfRPsIsGauss = kFALSE; // 1.) if kTRUE  = elliptic flow of RPs is sampled e-b-e from Gaussian distribution with
+Bool_t bV2DistrOfRPsIsGauss = kTRUE; // 1.) if kTRUE  = elliptic flow of RPs is sampled e-b-e from Gaussian distribution with
                                       //                 mean = dV2RP and spread = dV2SpreadRP
                                       // 2.) if kFALSE = elliptic flow of RPs is sampled e-b-e uniformly from 
                                       //                 interval [dMinV2RP,dMaxV2RP]
                                       // 3.) for a fixed elliptic flow use Gaussian with zero spread or use uniform with dMinV2RP=dMaxV2RP
-                                                                                   
+
+Bool_t uniformAcceptance = kTRUE; // 1.) if kTRUE = detectors has uniform azimuthal acceptance
+                                  // 2.) if kFALSE = you will simulate detector with non-uniform acceptance in one or two sectors. 
+                                  //                 For each of two sectors you specify phi_min, phi_max and probability p. Then all particles 
+                                  //                 going in direction phi_min < phi < phi_max will be taken with probability p. If p = 0, that
+                                  //                 sector is blocked. Set bellow phimin1, phimax1, p1 for the first sector and phimin2, phimax2, p2 
+                                  //                 for the second sector. If you set phimin2 = phimax2 = p2 = 0, only first non-uniform sector is 
+                                  //                 simulated.
+                                                                                                                                                                                                                                                          
 Int_t iMultiplicityOfRP = 500;        // mean multiplicity of RPs (if sampled from Gaussian)
 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs (if sampled from Gaussian)
 Int_t iMinMultOfRP = 400;             // minimal multiplicity of RPs (if sampled uniformly)
@@ -73,12 +81,28 @@ Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed)
 Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
 //......................................................................................  
 
+//......................................................................................  
+// settings for non-uniform acceptance:
+// Remark: set the angles in degrees from interval [0,360] and probability from interval [0,1]
+
+// 1st non-uniform sector:
+Double_t phimin1 = 60;  // first non-uniform sector starts at this azimuth
+Double_t phimax1 = 120; // first non-uniform sector ends at this azimuth
+Double_t p1 = 0.33;        // e.g. if p1 = 0 all particles emitted in phimin1 < phi < phimax1 are blocked
+                        // e.g. if p1 = 0.5 half of the particles emitted in phimin1 < phi < phimax1 are blocked
+
+// 2nd non-uniform sector (Remark: if you do NOT want to simulate this sector, set phimin2 = phimax2 = p2 = 0):                 
+Double_t phimin2 = 0.0; // second non-uniform sector starts at this azimuth (make sure phimin2 > phimax1 !!!!)
+Double_t phimax2 = 0.0; // second non-uniform sector ends at this azimuth
+Double_t p2 = 0.0;
+//......................................................................................  
+
 enum anaModes {mLocal,mLocalSource,mLocalPAR};
 // mLocal: Analyze data on your computer using aliroot
 // mLocalPAR: Analyze data on your computer using root + PAR files
 // mLocalSource: Analyze data on your computer using root + source files
                                           
-int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=10)
+int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=440)
 {
  TStopwatch timer;
  timer.Start();
@@ -87,6 +111,31 @@ int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=10)
  if (LYZ1PROD && LYZ2PROD)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
  if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
+
+ if(!uniformAcceptance && phimin1 > phimax1)
+ {
+  cout<<"WARNING: you must have phimin1 < phimax1 !!!!"<<endl;
+  break;
+ }
+
+ if (!uniformAcceptance && !((phimin2 == 0.) && (phimax2 == 0.) && (p2 == 0.)) && (phimin2 < phimax1 || phimin2 > phimax2))
+ {
+  cout<<"WARNING: you must have phimin2 > phimax1 and phimin2 < phimax2 !!!!"<<endl;
+  break;
+ }
+ if((phimin1 < 0 || phimin1 > 360) || (phimax1 < 0 || phimax1 > 360) || 
+    (phimin2 < 0 || phimin2 > 360) || (phimax2 < 0 || phimax2 > 360) )
+ {
+  cout<<"WARNING: you must take azimuthal angles from interval [0,360] !!!!"<<endl;
+  break;
+ }
+ if((p1 < 0 || p1 > 1) || (p2 < 0 || p2 > 1))
+ {
+  cout<<"WARNING: you must take p1 and p2 from interval [0,1] !!!!"<<endl;
+  break;
+ }
  
  cout<<endl;
  cout<<endl;
@@ -114,6 +163,10 @@ int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=10)
   sseed = dt.GetNanoSec()/1000;
  }
  
+ cout<<endl;
+ cout<<"Seed for the random generators is "<<sseed<<endl;
+ cout<<endl;
  //---------------------------------------------------------------------------------------
  // If the weights are used: 
  TFile *fileWithWeights = NULL;
@@ -308,6 +361,7 @@ int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=10)
      eventMakerOnTheFly->SetMaxV2RP(dMaxV2RP);  
     }
  }
  // (pt,eta) dependent harmonic V2:
  if(!bConstantHarmonics)
  {
@@ -315,6 +369,17 @@ int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=10)
   eventMakerOnTheFly->SetV2RPMax(dV2RPMax);
   eventMakerOnTheFly->SetPtCutOff(dPtCutOff);  
  }
+ // non-uniform acceptance:
+ if(!uniformAcceptance)
+ {
+  eventMakerOnTheFly->SetFirstSectorPhiMin(phimin1);
+  eventMakerOnTheFly->SetFirstSectorPhiMax(phimax1);
+  eventMakerOnTheFly->SetFirstSectorProbability(p1);
+  eventMakerOnTheFly->SetSecondSectorPhiMin(phimin2);
+  eventMakerOnTheFly->SetSecondSectorPhiMax(phimax2);
+  eventMakerOnTheFly->SetSecondSectorProbability(p2);
+ }
        
  //---------------------------------------------------------------------------------------  
  // create and analyze events 'on the fly':