remove old function
[u/mrichter/AliRoot.git] / STARLIGHT / AliGenStarLight.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id:  $ */
17 #include <Riostream.h>
18 #include "AliLog.h"
19 #include "TStarLight.h"
20 #include "starlight.h"
21 #include "upcevent.h"
22 #include "AliRunLoader.h"
23 #include "AliGenEventHeader.h"
24 #include "AliGenStarLight.h"
25
26 ClassImp(AliGenStarLight);
27
28 AliGenStarLight::AliGenStarLight()
29   : AliGenMC()
30   , fRapidityMotherMin( 1) // Max < Min: no cut
31   , fRapidityMotherMax(-1)
32   , fEtaChildMin( 1)       // Max < Min: no cut
33   , fEtaChildMax(-1)
34   , fSLgenerator(NULL)
35 {
36 }
37 //----------------------------------------------------------------------
38 AliGenStarLight::AliGenStarLight(Int_t npart)
39   : AliGenMC(npart)
40   , fRapidityMotherMin( 1) // Max < Min: no cut
41   , fRapidityMotherMax(-1)
42   , fEtaChildMin( 1)       // Max < Min: no cut
43   , fEtaChildMax(-1)
44   , fSLgenerator(new TStarLight("TStarLight",
45                                 "StarLight UPC Event Generator",
46                                 ""))  // no config file name  
47 {
48   //
49 }
50 //----------------------------------------------------------------------
51 AliGenStarLight::~AliGenStarLight() {
52   if (NULL != fSLgenerator) delete fSLgenerator;
53   fSLgenerator = NULL;
54 }
55 void AliGenStarLight::ImportConfigurationFromFile(const char* filename) {
56   if (NULL == fSLgenerator) {
57     AliFatal("AliGenStarLight class not constructed properly. ");
58     return;
59   }
60   fSLgenerator->ImportConfigurationFromFile(filename);
61 }
62 void AliGenStarLight::SetParameter(const char* line) {
63   if (NULL == fSLgenerator) {
64     AliFatal("AliGenStarLight class not constructed properly. ");
65     return;
66   }
67   fSLgenerator->SetParameter(line);
68 }
69 //----------------------------------------------------------------------
70 void AliGenStarLight::Init() {
71   if (NULL == fSLgenerator) {
72     AliFatal("AliGenStarLight class not constructed properly. ");
73     return;
74   }
75   fSLgenerator->InitStarLight();
76 }
77 //----------------------------------------------------------------------
78 void AliGenStarLight::Generate() {
79   Float_t vpos[4] = { 0, 0, 0, 0 };
80   if (fVertexSmear == kPerEvent) {
81     Vertex(); // get vertex
82     for (Int_t i=0; i<3; ++i)
83       vpos[i] = fVertex[i];
84     vpos[3] = fTime;
85   }
86
87   Int_t   nt(0);     // number of tracks
88   Bool_t  genOK(kFALSE);  
89   // generate events until all constraints are fulfilled
90   for (Int_t trials=0; !genOK && trials < 100*1000; ++trials) {
91     fSLgenerator->GenerateEvent();
92     fSLgenerator->BoostEvent();
93     fSLgenerator->ImportParticles(&fParticles, "ALL");
94     
95     TLorentzVector vSum;
96     genOK = kTRUE;
97     const Long64_t n(fParticles.GetEntries());
98     if (n == 0) {
99       AliFatal("no particles generated");
100       return;
101     }
102     for (Long64_t i(0); i<n; ++i) {
103       TParticle *part(dynamic_cast<TParticle*>(fParticles.At(i)));
104       if (NULL == part) {
105         AliFatal("NULL == part");
106         return;
107       }
108       genOK = 
109         (part->Theta() >= fThetaMin) && (part->Theta() <  fThetaMax) &&
110         (part->Phi()   >= fPhiMin)   && (part->Phi()   <  fPhiMax)   &&
111         (part->Y()     >= fYMin)     && (part->Y()     <  fYMax)     &&
112         (part->P()     >= fPMin)     && (part->P()     <  fPMax)     &&
113         (part->Pt()    >= fPtMin)    && (part->Pt()    <  fPtMax);
114       if (fEtaChildMin <= fEtaChildMax) // no cut if Max < Min
115         genOK = genOK && (part->Eta() >= fEtaChildMin &&
116                           part->Eta() <  fEtaChildMax);
117       if (kFALSE == genOK)
118         break;
119
120       TLorentzVector v;
121       part->Momentum(v);
122       vSum += v;
123     }
124
125     if (fRapidityMotherMin <= fRapidityMotherMax) // no cut if Max < Min
126       genOK = (genOK
127                ? (vSum.Rapidity() > fRapidityMotherMin &&
128                   vSum.Rapidity() < fRapidityMotherMax)
129                : kFALSE);
130
131     if (kFALSE == genOK) continue;
132     
133     fNprimaries = 0;
134     for (Long64_t i(0), n(fParticles.GetEntries()); i<n; ++i) {
135       const TParticle *part(dynamic_cast<TParticle*>(fParticles.At(i)));        
136       if (NULL == part) {
137         AliFatal("NULL == part");
138         return;
139       }
140       const Int_t   iparent(-1);
141       const Float_t polar[3] = { 0, 0, 0 };
142       const Float_t weight(trials+1);
143       PushTrack(fTrackIt, iparent, part->GetPdgCode(), 
144                 part->Px(), part->Py(), part->Pz(), part->Energy(),
145                 vpos[0],    vpos[1],    vpos[2],    vpos[3],
146                 polar[0],   polar[1],   polar[2], 
147                 kPPrimary, nt, weight, part->GetStatusCode());
148       AliInfo(Form("weight=%.0f nt=%d fTrackIt=%d statusCode=%d",
149                    weight, nt, fTrackIt, part->GetStatusCode()));
150       part->Print();
151       KeepTrack(nt);
152       ++fNprimaries;
153     }
154     fParticles.Clear();
155   }
156   if (kFALSE == genOK)
157     AliFatal("Maximum number of trials reached");
158
159   AliGenEventHeader* header = new AliGenEventHeader("SL");
160   const TArrayF vertexPosition(3, vpos);
161   header->SetPrimaryVertex(vertexPosition);
162   header->SetInteractionTime(vpos[3]);
163   header->SetNProduced(nt);
164   AddHeader(header);
165   SetHighWaterMark(nt);
166   AliRunLoader::Instance()->CdGAFile();
167 }