X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliRun.cxx;h=5ce56e1003559bec5a711ef999e71a899b904f02;hb=b9d0a01d7a0723a09071b0b56200d72f59a9c2b6;hp=c0e4580535d6fd96da3b68ec8ce8dc37254f35a3;hpb=7fb014804a171990cc00803e61957e4bbe893f38;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliRun.cxx b/STEER/AliRun.cxx index c0e4580535d..5ce56e10035 100644 --- a/STEER/AliRun.cxx +++ b/STEER/AliRun.cxx @@ -15,6 +15,208 @@ /* $Log$ +Revision 1.82.4.3 2002/10/14 09:45:57 hristov +Updating VirtualMC to v3-09-02 + +Revision 1.82.4.2 2002/06/10 17:54:06 hristov +Only one SetGenEventHeader function kept + +Revision 1.82.4.1 2002/06/10 14:43:06 hristov +Merged with v3-08-02 + +Revision 1.86 2002/05/28 14:24:57 hristov +Correct warning messages + +Revision 1.85 2002/05/24 13:29:58 hristov +AliTrackReference added, AliDisplay modified + +Revision 1.84 2002/05/21 16:26:07 hristov +Find correctly TreeK in case CONFIG_SPLIT_FILE is set (Y.Schutz) + +Revision 1.83 2002/04/04 13:16:17 jchudoba +add possibility to write sdigits, digits and rec. points into separate files + +Revision 1.82 2002/03/12 11:06:03 morsch +Add particle status code to argument list of SetTrack(..). + +Revision 1.81 2001/12/19 14:46:26 morsch +Add possibility to disable StepManager() for each module separately. + +Revision 1.80 2001/10/21 18:22:55 hristov +BranchOld replaced by Branch. It works correctly with Root 2.02.xx + +Revision 1.79 2001/10/09 18:00:35 hristov +Temporary fix to provide unique event number in the simulation (J.Chudoba) + +Revision 1.78 2001/10/04 15:30:56 hristov +Changes to accommodate the set of PHOS folders and tasks (Y.Schutz) + +Revision 1.77 2001/09/04 15:09:11 hristov +fTreeE->Branch replaced temporary by fTreeE->BranchOld to avoid data corruption in case of many events per file + +Revision 1.76 2001/08/03 14:38:35 morsch +Use original method to access TreeH. + +Revision 1.75 2001/07/28 10:52:27 hristov +Event number updated correctly (M.Ivanov) + +Revision 1.74 2001/07/28 10:39:16 morsch +GetEventsPerRun() method needed by afterburners added to AliRun.h +Corresponding setters and getters have been from AliGenerator. + +Revision 1.73 2001/07/27 12:34:30 jchudoba +remove the dummy argument in fStack->GetEvent call + +Revision 1.72 2001/07/03 08:10:57 hristov +J.Chudoba's changes merged correctly with the HEAD + +Revision 1.70 2001/06/29 08:01:36 morsch +Small correction to the previous. + +Revision 1.69 2001/06/28 16:27:50 morsch +AliReco() with user control of event range. + +Revision 1.68 2001/06/11 13:14:40 morsch +SetAliGenEventHeader() method added. + +Revision 1.67 2001/06/07 18:24:50 buncic +Removed compilation warning in AliConfig initialisation. + +Revision 1.66 2001/05/22 14:32:40 hristov +Weird inline removed + +Revision 1.65 2001/05/21 17:22:51 buncic +Fixed problem with missing AliConfig while reading galice.root + +Revision 1.64 2001/05/16 14:57:22 alibrary +New files for folders and Stack + +Revision 1.62 2001/04/06 11:12:33 morsch +Clear fParticles after each event. (Ivana Hrivnacova) + +Revision 1.61 2001/03/30 07:04:10 morsch +Call fGenerator->FinishRun() for final print-outs, cross-section and weight calculations. + +Revision 1.60 2001/03/21 18:22:30 hristov +fParticleFileMap fix (I.Hrivnacova) + +Revision 1.59 2001/03/12 17:47:03 hristov +Changes needed on Sun with CC 5.0 + +Revision 1.58 2001/03/09 14:27:26 morsch +Fix for multiple events per file: inhibit decrease of size of fParticleFileMap. + +Revision 1.57 2001/02/23 17:40:23 buncic +All trees needed for simulation created in RunMC(). TreeR and its branches +are now created in new RunReco() method. + +Revision 1.56 2001/02/14 15:45:20 hristov +Algorithmic way of getting entry index in fParticleMap. Protection of fParticleFileMap (I.Hrivnacova) + +Revision 1.55 2001/02/12 15:52:54 buncic +Removed OpenBaseFile(). + +Revision 1.54 2001/02/07 10:39:05 hristov +Remove default value for argument + +Revision 1.53 2001/02/06 11:02:26 hristov +New SetTrack interface added, added check for unfilled particles in FinishEvent (I.Hrivnacova) + +Revision 1.52 2001/02/05 16:22:25 buncic +Added TreeS to GetEvent(). + +Revision 1.51 2001/02/02 15:16:20 morsch +SetHighWaterMark method added to mark last particle in event. + +Revision 1.50 2001/01/27 10:32:00 hristov +Leave the loop when primaries are filled (I.Hrivnacova) + +Revision 1.49 2001/01/26 19:58:48 hristov +Major upgrade of AliRoot code + +Revision 1.48 2001/01/17 10:50:50 hristov +Corrections to destructors + +Revision 1.47 2000/12/18 10:44:01 morsch +Possibility to set field map by passing pointer to objet of type AliMagF via +SetField(). +Example: +gAlice->SetField(new AliMagFCM("Map2", "$(ALICE_ROOT)/data/field01.dat",2,1.,10.)); + +Revision 1.46 2000/12/14 19:29:27 fca +galice.cuts was not read any more + +Revision 1.45 2000/11/30 07:12:49 alibrary +Introducing new Rndm and QA classes + +Revision 1.44 2000/10/26 13:58:59 morsch +Add possibility to choose the lego generator (of type AliGeneratorLego or derived) when running +RunLego(). Default is the base class AliGeneratorLego. + +Revision 1.43 2000/10/09 09:43:17 fca +Special remapping of hits for TPC and TRD. End-of-primary action introduced + +Revision 1.42 2000/10/02 21:28:14 fca +Removal of useless dependecies via forward declarations + +Revision 1.41 2000/07/13 16:19:09 fca +Mainly coding conventions + some small bug fixes + +Revision 1.40 2000/07/12 08:56:25 fca +Coding convention correction and warning removal + +Revision 1.39 2000/07/11 18:24:59 fca +Coding convention corrections + few minor bug fixes + +Revision 1.38 2000/06/20 13:05:45 fca +Writing down the TREE headers before job starts + +Revision 1.37 2000/06/09 20:05:11 morsch +Introduce possibility to chose magnetic field version 3: AliMagFDM + field02.dat + +Revision 1.36 2000/06/08 14:03:58 hristov +Only one initializer for a default argument + +Revision 1.35 2000/06/07 10:13:14 hristov +Delete only existent objects. + +Revision 1.34 2000/05/18 10:45:38 fca +Delete Particle Factory properly + +Revision 1.33 2000/05/16 13:10:40 fca +New method IsNewTrack and fix for a problem in Father-Daughter relations + +Revision 1.32 2000/04/27 10:38:21 fca +Correct termination of Lego Run and introduce Lego getter in AliRun + +Revision 1.31 2000/04/26 10:17:32 fca +Changes in Lego for G4 compatibility + +Revision 1.30 2000/04/18 19:11:40 fca +Introduce variable Config.C function signature + +Revision 1.29 2000/04/07 11:12:34 fca +G4 compatibility changes + +Revision 1.28 2000/04/05 06:51:06 fca +Workaround for an HP compiler problem + +Revision 1.27 2000/03/22 18:08:07 fca +Rationalisation of the virtual MC interfaces + +Revision 1.26 2000/03/22 13:42:26 fca +SetGenerator does not replace an existing generator, ResetGenerator does + +Revision 1.25 2000/02/23 16:25:22 fca +AliVMC and AliGeant3 classes introduced +ReadEuclid moved from AliRun to AliModule + +Revision 1.24 2000/01/19 17:17:20 fca +Introducing a list of lists of hits -- more hits allowed for detector now + +Revision 1.23 1999/12/03 11:14:31 fca +Fixing previous wrong checking + Revision 1.21 1999/11/25 10:40:08 fca Fixing daughters information also in primary tracks @@ -51,66 +253,68 @@ Introduction of the Copyright and cvs Log // // /////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include + #include #include #include -#include #include #include #include - +#include +#include +#include +#include +#include +#include +#include #include "TParticle.h" #include "AliRun.h" #include "AliDisplay.h" +#include "AliMC.h" +#include "AliLego.h" +#include "AliMagFC.h" +#include "AliMagFCM.h" +#include "AliMagFDM.h" +#include "AliHit.h" +#include "TRandom3.h" +#include "AliMCQA.h" +#include "AliGenerator.h" +#include "AliLegoGenerator.h" +#include "AliConfig.h" +#include "AliStack.h" +#include "AliGenEventHeader.h" +#include "AliHeader.h" + +#include "AliDetector.h" -#include "AliCallf77.h" - -#include -#include -#include - AliRun *gAlice; -static AliHeader *header; - -#ifndef WIN32 - -# define rxgtrak rxgtrak_ -# define rxstrak rxstrak_ -# define rxkeep rxkeep_ -# define rxouth rxouth_ -#else - -# define rxgtrak RXGTRAK -# define rxstrak RXSTRAK -# define rxkeep RXKEEP -# define rxouth RXOUTH -#endif - -static TArrayF sEventEnergy; -static TArrayF sSummEnergy; -static TArrayF sSum2Energy; - ClassImp(AliRun) //_____________________________________________________________________________ AliRun::AliRun() + : TVirtualMCApplication() { // // Default constructor for AliRun // - header=&fHeader; + fHeader = 0; fRun = 0; fEvent = 0; - fCurrent = -1; - fModules = 0; + fEventNrInRun = 0; + fStack = 0; + fModules = 0; fGenerator = 0; fTreeD = 0; - fTreeK = 0; fTreeH = 0; + fTreeTR = 0; fTreeE = 0; fTreeR = 0; - fParticles = 0; + fTreeS = 0; fGeometry = 0; fDisplay = 0; fField = 0; @@ -122,11 +326,24 @@ AliRun::AliRun() fInitDone = kFALSE; fLego = 0; fPDGDB = 0; //Particle factory object! + fHitLists = 0; + fConfigFunction = "\0"; + fRandom = 0; + fMCQA = 0; + fTransParName = "\0"; + fBaseFileName = ".\0"; + fDebug = 0; + fTreeDFile = 0; + fTreeSFileName = ""; + fTreeSFile = 0; + fTreeSFileName = ""; + fTreeRFile = 0; + fTreeRFileName = ""; } //_____________________________________________________________________________ AliRun::AliRun(const char *name, const char *title) - : TNamed(name,title) + : TVirtualMCApplication(name,title) { // // Constructor for the main processor. @@ -138,19 +355,36 @@ AliRun::AliRun(const char *name, const char *title) gAlice = this; fTreeD = 0; - fTreeK = 0; fTreeH = 0; + fTreeTR = 0; fTreeE = 0; fTreeR = 0; + fTreeS = 0; fTrRmax = 1.e10; fTrZmax = 1.e10; fGenerator = 0; fInitDone = kFALSE; fLego = 0; fField = 0; + fConfigFunction = "Config();"; + fTreeDFile = 0; + fTreeSFileName = ""; + fTreeSFile = 0; + fTreeSFileName = ""; + fTreeRFile = 0; + fTreeRFileName = ""; + + // Set random number generator + gRandom = fRandom = new TRandom3(); + + if (gSystem->Getenv("CONFIG_SEED")) { + gRandom->SetSeed((UInt_t)atoi(gSystem->Getenv("CONFIG_SEED"))); + } gROOT->GetListOfBrowsables()->Add(this,name); // + // Particle stack + fStack = new AliStack(10000); // create the support list for the various Detectors fModules = new TObjArray(77); // @@ -158,17 +392,11 @@ AliRun::AliRun(const char *name, const char *title) BuildSimpleGeometry(); - - fNtrack=0; - fHgwmk=0; - fCurrent=-1; - header=&fHeader; + fHeader = new AliHeader(); fRun = 0; fEvent = 0; + fEventNrInRun = 0; // - // Create the particle stack - fParticles = new TClonesArray("TParticle",100); - fDisplay = 0; // // Create default mag field @@ -182,14 +410,27 @@ AliRun::AliRun(const char *name, const char *title) // // Make particles fPDGDB = TDatabasePDG::Instance(); //Particle factory object! + + AliConfig::Instance()->Add(fPDGDB); + // + // Create HitLists list + fHitLists = new TList(); + // + SetTransPar(); + fBaseFileName = ".\0"; + // + fDebug = 0; } + //_____________________________________________________________________________ AliRun::~AliRun() { // - // Defaullt AliRun destructor + // Default AliRun destructor // + TFile *curfil =0; + if(fTreeE)curfil=fTreeE->GetCurrentFile(); delete fImedia; delete fField; delete fMC; @@ -198,18 +439,44 @@ AliRun::~AliRun() delete fGenerator; delete fLego; delete fTreeD; - delete fTreeK; delete fTreeH; + delete fTreeTR; delete fTreeE; delete fTreeR; + delete fTreeS; if (fModules) { fModules->Delete(); delete fModules; } - if (fParticles) { - fParticles->Delete(); - delete fParticles; + delete fStack; + delete fHitLists; + delete fPDGDB; + delete fMCQA; + delete fHeader; + // avoid to delete TFile objects not owned by this object + // avoid multiple deletions + if(curfil == fTreeDFile) fTreeDFile=0; + if(curfil == fTreeSFile) fTreeSFile=0; + if(curfil == fTreeRFile) fTreeRFile=0; + if(fTreeSFile == fTreeDFile) fTreeSFile=0; + if(fTreeRFile == fTreeDFile) fTreeRFile=0; + if(fTreeRFile == fTreeSFile) fTreeRFile=0; + if(fTreeDFile){ + if(fTreeDFile->IsOpen())fTreeDFile->Close(); + delete fTreeDFile; } + if(fTreeSFile){ + if(fTreeSFile->IsOpen())fTreeSFile->Close(); + delete fTreeSFile; + } + if(fTreeRFile){ + if(fTreeRFile->IsOpen())fTreeRFile->Close(); + delete fTreeRFile; + } + if (gROOT->GetListOfBrowsables()) + gROOT->GetListOfBrowsables()->Remove(this); + + gAlice=0; } //_____________________________________________________________________________ @@ -240,17 +507,23 @@ void AliRun::Browse(TBrowser *b) // of the Root browser. // It displays the Root Trees and all detectors. // - if (fTreeK) b->Add(fTreeK,fTreeK->GetName()); + if(!fStack) fStack=fHeader->Stack(); + TTree* pTreeK = fStack->TreeK(); + + if (pTreeK) b->Add(pTreeK,pTreeK->GetName()); if (fTreeH) b->Add(fTreeH,fTreeH->GetName()); + if (fTreeTR) b->Add(fTreeTR,fTreeH->GetName()); if (fTreeD) b->Add(fTreeD,fTreeD->GetName()); if (fTreeE) b->Add(fTreeE,fTreeE->GetName()); if (fTreeR) b->Add(fTreeR,fTreeR->GetName()); + if (fTreeS) b->Add(fTreeS,fTreeS->GetName()); TIter next(fModules); AliModule *detector; while((detector = (AliModule*)next())) { b->Add(detector,detector->GetName()); } + b->Add(fMCQA,"AliMCQA"); } //_____________________________________________________________________________ @@ -290,25 +563,6 @@ void AliRun::CleanDetectors() } } -//_____________________________________________________________________________ -void AliRun::CleanParents() -{ - // - // Clean Particles stack. - // Set parent/daughter relations - // - TClonesArray &particles = *(gAlice->Particles()); - TParticle *part; - int i; - for(i=0; iTestBit(Daughters_Bit)) { - part->SetFirstDaughter(-1); - part->SetLastDaughter(-1); - } - } -} - //_____________________________________________________________________________ Int_t AliRun::DistancetoPrimitive(Int_t, Int_t) { @@ -320,31 +574,29 @@ Int_t AliRun::DistancetoPrimitive(Int_t, Int_t) } //_____________________________________________________________________________ -void AliRun::DumpPart (Int_t i) +void AliRun::DumpPart (Int_t i) const { // // Dumps particle i in the stack // - TClonesArray &particles = *fParticles; - ((TParticle*) particles[i])->Print(); + fStack->DumpPart(i); } //_____________________________________________________________________________ -void AliRun::DumpPStack () +void AliRun::DumpPStack () const { // // Dumps the particle stack // - TClonesArray &particles = *fParticles; - printf( - "\n\n=======================================================================\n"); - for (Int_t i=0;i %d ",i); ((TParticle*) particles[i])->Print(); - printf("--------------------------------------------------------------\n"); - } - printf( - "\n=======================================================================\n\n"); + fStack->DumpPStack(); +} + +//_____________________________________________________________________________ +void AliRun::SetField(AliMagF* magField) +{ + // Set Magnetic Field Map + fField = magField; + fField->ReadField(); } //_____________________________________________________________________________ @@ -360,113 +612,20 @@ void AliRun::SetField(Int_t type, Int_t version, Float_t scale, // // --- Sanity check on mag field flags - if(type<0 || type > 2) { - Warning("SetField", - "Invalid magnetic field flag: %5d; Helix tracking chosen instead\n" - ,type); - type=2; - } if(fField) delete fField; if(version==1) { - fField = new AliMagFC("Map1"," ",type,version,scale,maxField); - } else if(version<=3) { - fField = new AliMagFCM("Map2-3",filename,type,version,scale,maxField); + fField = new AliMagFC("Map1"," ",type,scale,maxField); + } else if(version<=2) { + fField = new AliMagFCM("Map2-3",filename,type,scale,maxField); + fField->ReadField(); + } else if(version==3) { + fField = new AliMagFDM("Map4",filename,type,scale,maxField); fField->ReadField(); } else { Warning("SetField","Invalid map %d\n",version); } } - -//_____________________________________________________________________________ -void AliRun::FillTree() -{ - // - // Fills all AliRun TTrees - // - if (fTreeK) fTreeK->Fill(); - if (fTreeH) fTreeH->Fill(); - if (fTreeD) fTreeD->Fill(); - if (fTreeR) fTreeR->Fill(); -} -//_____________________________________________________________________________ -void AliRun::FinishPrimary() -{ - // - // Called at the end of each primary track - // - - // static Int_t count=0; - // const Int_t times=10; - // This primary is finished, purify stack - gAlice->PurifyKine(); - - // Write out hits if any - if (gAlice->TreeH()) { - gAlice->TreeH()->Fill(); - } - - // Reset Hits info - gAlice->ResetHits(); - - // - // if(++count%times==1) gObjectTable->Print(); -} - -//_____________________________________________________________________________ -void AliRun::FinishEvent() -{ - // - // Called at the end of the event. - // - - //Update the energy deposit tables - Int_t i; - for(i=0;iFill(); - } - - // Write out the digits - if (fTreeD) { - fTreeD->Fill(); - ResetDigits(); - } - - // Write out reconstructed clusters - if (fTreeR) { - fTreeR->Fill(); - } - - // Write out the event Header information - if (fTreeE) fTreeE->Fill(); - - // Reset stack info - ResetStack(); - - // Write Tree headers - // Int_t ievent = fHeader.GetEvent(); - // char hname[30]; - // sprintf(hname,"TreeK%d",ievent); - if (fTreeK) fTreeK->Write(); - // sprintf(hname,"TreeH%d",ievent); - if (fTreeH) fTreeH->Write(); - // sprintf(hname,"TreeD%d",ievent); - if (fTreeD) fTreeD->Write(); - // sprintf(hname,"TreeR%d",ievent); - if (fTreeR) fTreeR->Write(); -} - //_____________________________________________________________________________ void AliRun::FinishRun() { @@ -474,6 +633,9 @@ void AliRun::FinishRun() // Called at the end of the run. // + // + if(fLego) fLego->FinishRun(); + // Clean detector information TIter next(fModules); AliModule *detector; @@ -483,57 +645,50 @@ void AliRun::FinishRun() //Output energy summary tables EnergySummary(); + + TFile *file = fTreeE->GetCurrentFile(); + + file->cd(); - // file is retrieved from whatever tree - TFile *File = 0; - if (fTreeK) File = fTreeK->GetCurrentFile(); - if ((!File) && (fTreeH)) File = fTreeH->GetCurrentFile(); - if ((!File) && (fTreeD)) File = fTreeD->GetCurrentFile(); - if ((!File) && (fTreeE)) File = fTreeE->GetCurrentFile(); - if( NULL==File ) { - Error("FinishRun","There isn't root file!"); - exit(1); - } - File->cd(); - fTreeE->Write(); + fTreeE->Write(0,TObject::kOverwrite); + // Write AliRun info and all detectors parameters + Write(0,TObject::kOverwrite); + // Clean tree information - delete fTreeK; fTreeK = 0; - delete fTreeH; fTreeH = 0; - delete fTreeD; fTreeD = 0; - delete fTreeR; fTreeR = 0; - delete fTreeE; fTreeE = 0; + + fStack->FinishRun(); - // Write AliRun info and all detectors parameters - Write(); + if (fTreeH) { + delete fTreeH; fTreeH = 0; + } + if (fTreeTR) { + delete fTreeTR; fTreeTR = 0; + } + if (fTreeD) { + delete fTreeD; fTreeD = 0; + } + if (fTreeR) { + delete fTreeR; fTreeR = 0; + } +// if (fTreeE) { +// delete fTreeE; fTreeE = 0; +// } + if (fTreeS) { + delete fTreeS; fTreeS = 0; + } + fGenerator->FinishRun(); // Close output file - File->Write(); - File->Close(); + file->Write(); } //_____________________________________________________________________________ void AliRun::FlagTrack(Int_t track) { + // Delegate to stack // - // Flags a track and all its family tree to be kept - // - int curr; - TParticle *particle; - - curr=track; - while(1) { - particle=(TParticle*)fParticles->UncheckedAt(curr); - - // If the particle is flagged the three from here upward is saved already - if(particle->TestBit(Keep_Bit)) return; - - // Save this particle - particle->SetBit(Keep_Bit); - - // Move to father if any - if((curr=particle->GetFirstMother())==-1) return; - } + fStack->FlagTrack(track); } //_____________________________________________________________________________ @@ -547,25 +702,25 @@ void AliRun::EnergySummary() Float_t edtot=0; Float_t ed, ed2; Int_t kn, i, left, j, id; - const Float_t zero=0; - Int_t ievent=fHeader.GetEvent()+1; + const Float_t kzero=0; + Int_t ievent=fHeader->GetEvent()+1; // // Energy loss information if(ievent) { printf("***************** Energy Loss Information per event (GEV) *****************\n"); - for(kn=1;kn0) { - sEventEnergy[ndep]=kn; + fEventEnergy[ndep]=kn; if(ievent>1) { ed=ed/ievent; - ed2=sSum2Energy[kn]; + ed2=fSum2Energy[kn]; ed2=ed2/ievent; - ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,zero))/ed; + ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed; } else ed2=99; - sSummEnergy[ndep]=ed; - sSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,zero)); + fSummEnergy[ndep]=ed; + fSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,kzero)); edtot+=ed; ndep++; } @@ -574,8 +729,8 @@ void AliRun::EnergySummary() left=ndep-kn*3; for(i=0;i<(3VolName(id),sSummEnergy[j],sSum2Energy[j]); + id=Int_t (fEventEnergy[j]+0.1); + printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]); } printf("\n"); } @@ -587,8 +742,8 @@ void AliRun::EnergySummary() left=ndep-kn*5; for(i=0;i<(5VolName(id),100*sSummEnergy[j]/edtot); + id=Int_t (fEventEnergy[j]+0.1); + printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot); } printf("\n"); } @@ -597,13 +752,13 @@ void AliRun::EnergySummary() } // // Reset the TArray's - sEventEnergy.Set(0); - sSummEnergy.Set(0); - sSum2Energy.Set(0); + // fEventEnergy.Set(0); + // fSummEnergy.Set(0); + // fSum2Energy.Set(0); } //_____________________________________________________________________________ -AliModule *AliRun::GetModule(const char *name) +AliModule *AliRun::GetModule(const char *name) const { // // Return pointer to detector from name @@ -612,7 +767,7 @@ AliModule *AliRun::GetModule(const char *name) } //_____________________________________________________________________________ -AliDetector *AliRun::GetDetector(const char *name) +AliDetector *AliRun::GetDetector(const char *name) const { // // Return pointer to detector from name @@ -621,7 +776,7 @@ AliDetector *AliRun::GetDetector(const char *name) } //_____________________________________________________________________________ -Int_t AliRun::GetModuleID(const char *name) +Int_t AliRun::GetModuleID(const char *name) const { // // Return galice internal detector identifier from name @@ -641,48 +796,135 @@ Int_t AliRun::GetEvent(Int_t event) // // Reset existing structures - ResetStack(); ResetHits(); + ResetTrackReferences(); ResetDigits(); + ResetSDigits(); // Delete Trees already connected - if (fTreeK) delete fTreeK; - if (fTreeH) delete fTreeH; - if (fTreeD) delete fTreeD; - if (fTreeR) delete fTreeR; - - // Get header from file - if(fTreeE) fTreeE->GetEntry(event); - else Error("GetEvent","Cannot file Header Tree\n"); + if (fTreeH) { delete fTreeH; fTreeH = 0;} + if (fTreeTR) { delete fTreeTR; fTreeTR = 0;} + if (fTreeD) { delete fTreeD; fTreeD = 0;} + if (fTreeR) { delete fTreeR; fTreeR = 0;} + if (fTreeS) { delete fTreeS; fTreeS = 0;} + + // Create the particle stack + if (fHeader) delete fHeader; + fHeader = 0; - // Get Kine Tree from file + // Get header from file + if(fTreeE) { + fTreeE->SetBranchAddress("Header", &fHeader); + + if (!fTreeE->GetEntry(event)) { + Error("GetEvent","Cannot find event:%d\n",event); + return -1; + } + } + else { + Error("GetEvent","Cannot find Header Tree (TE)\n"); + return -1; + } + + // Get the stack from the header, set fStack to 0 if it + // fails to get event + // + TFile *file = fTreeE->GetCurrentFile(); char treeName[20]; - sprintf(treeName,"TreeK%d",event); - fTreeK = (TTree*)gDirectory->Get(treeName); - if (fTreeK) fTreeK->SetBranchAddress("Particles", &fParticles); - else Error("GetEvent","cannot find Kine Tree for event:%d\n",event); + file->cd(); + + if (fStack) delete fStack; + fStack = fHeader->Stack(); + if (fStack) { + if (!fStack->GetEvent(event)) fStack = 0; + } + // Get Hits Tree header from file sprintf(treeName,"TreeH%d",event); fTreeH = (TTree*)gDirectory->Get(treeName); if (!fTreeH) { - Error("GetEvent","cannot find Hits Tree for event:%d\n",event); + Warning("GetEvent","cannot find Hits Tree for event:%d\n",event); } - + + // Get TracReferences Tree header from file + sprintf(treeName,"TreeTR%d",event); + fTreeTR = (TTree*)gDirectory->Get(treeName); + if (!fTreeTR) { + Warning("GetEvent","cannot find TrackReferences Tree for event:%d\n",event); + } + + // get current file name and compare with names containing trees S,D,R + TString curfilname=(TString)fTreeE->GetCurrentFile()->GetName(); + if(fTreeDFileName==curfilname)fTreeDFileName=""; + if(fTreeSFileName==curfilname)fTreeSFileName=""; + if(fTreeRFileName==curfilname)fTreeRFileName=""; + // Get Digits Tree header from file sprintf(treeName,"TreeD%d",event); - fTreeD = (TTree*)gDirectory->Get(treeName); + + if (!fTreeDFile && fTreeDFileName != "") { + InitTreeFile("D",fTreeDFileName); + } + if (fTreeDFile) { + fTreeD = (TTree*)fTreeDFile->Get(treeName); + } else { + fTreeD = (TTree*)file->Get(treeName); + } if (!fTreeD) { - Warning("GetEvent","cannot find Digits Tree for event:%d\n",event); + // Warning("GetEvent","cannot find Digits Tree for event:%d\n",event); } - + if(fTreeDFileName != ""){ + if(fTreeDFileName==fTreeSFileName) { + fTreeSFileName = ""; + fTreeSFile = fTreeDFile; + } + if(fTreeDFileName==fTreeRFileName) { + fTreeRFileName = ""; + fTreeRFile = fTreeDFile; + } + } + + file->cd(); + + // Get SDigits Tree header from file + sprintf(treeName,"TreeS%d",event); + if (!fTreeSFile && fTreeSFileName != "") { + InitTreeFile("S",fTreeSFileName); + } + if (fTreeSFile) { + fTreeS = (TTree*)fTreeSFile->Get(treeName); + } else { + fTreeS = (TTree*)gDirectory->Get(treeName); + } + if (!fTreeS) { + // Warning("GetEvent","cannot find SDigits Tree for event:%d\n",event); + } + + if(fTreeSFileName != ""){ + if(fTreeSFileName==fTreeRFileName){ + fTreeRFileName = ""; + fTreeRFile = fTreeSFile; + } + } + + file->cd(); // Get Reconstruct Tree header from file sprintf(treeName,"TreeR%d",event); - fTreeR = (TTree*)gDirectory->Get(treeName); + if (!fTreeRFile && fTreeRFileName != "") { + InitTreeFile("R",fTreeRFileName); + } + if(fTreeRFile) { + fTreeR = (TTree*)fTreeRFile->Get(treeName); + } else { + fTreeR = (TTree*)gDirectory->Get(treeName); + } if (!fTreeR) { // printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event); } + + file->cd(); // Set Trees branch addresses TIter next(fModules); @@ -690,10 +932,10 @@ Int_t AliRun::GetEvent(Int_t event) while((detector = (AliModule*)next())) { detector->SetTreeAddress(); } - - if (fTreeK) fTreeK->GetEvent(0); - fNtrack = Int_t (fParticles->GetEntries()); - return fNtrack; + + fEvent=event; //MI change + + return fHeader->GetNtrack(); } //_____________________________________________________________________________ @@ -712,7 +954,6 @@ TGeometry *AliRun::GetGeometry() TIter next(fModules); AliModule *detector; while((detector = (AliModule*)next())) { - detector->SetTreeAddress(); TList *dnodes=detector->Nodes(); Int_t j; TNode *node, *node1; @@ -727,143 +968,25 @@ TGeometry *AliRun::GetGeometry() } //_____________________________________________________________________________ -void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom, - Float_t &e, Float_t *vpos, Float_t *polar, - Float_t &tof) -{ - // - // Return next track from stack of particles - // - TVector3 pol; - fCurrent=-1; - TParticle *track; - for(Int_t i=fNtrack-1; i>=0; i--) { - track=(TParticle*) fParticles->UncheckedAt(i); - if(!track->TestBit(Done_Bit)) { - // - // The track has not yet been processed - fCurrent=i; - ipart=track->GetPdgCode(); - pmom[0]=track->Px(); - pmom[1]=track->Py(); - pmom[2]=track->Pz(); - e =track->Energy(); - vpos[0]=track->Vx(); - vpos[1]=track->Vy(); - vpos[2]=track->Vz(); - track->GetPolarisation(pol); - polar[0]=pol.X(); - polar[1]=pol.Y(); - polar[2]=pol.Z(); - tof=track->T(); - track->SetBit(Done_Bit); - break; - } - } - mtrack=fCurrent; - // - // stop and start timer when we start a primary track - Int_t nprimaries = fHeader.GetNprimary(); - if (fCurrent >= nprimaries) return; - if (fCurrent < nprimaries-1) { - fTimer.Stop(); - track=(TParticle*) fParticles->UncheckedAt(fCurrent+1); - // track->SetProcessTime(fTimer.CpuTime()); - } - fTimer.Start(); -} - -//_____________________________________________________________________________ -Int_t AliRun::GetPrimary(Int_t track) +Int_t AliRun::GetPrimary(Int_t track) const { // // return number of primary that has generated track // - int current, parent; - TParticle *part; - // - parent=track; - while (1) { - current=parent; - part = (TParticle *)fParticles->UncheckedAt(current); - parent=part->GetFirstMother(); - if(parent<0) return current; - } + return fStack->GetPrimary(track); } //_____________________________________________________________________________ -void AliRun::Init(const char *setup) +void AliRun::MediaTable() { // - // Initialize the Alice setup + // Built media table to get from the media number to + // the detector id // - - gROOT->LoadMacro(setup); - gInterpreter->ProcessLine("Config();"); - - gMC->DefineParticles(); //Create standard MC particles - - TObject *objfirst, *objlast; - - fNdets = fModules->GetLast()+1; - - // - //=================Create Materials, geometry, histograms, etc - TIter next(fModules); - AliModule *detector; - while((detector = (AliModule*)next())) { - detector->SetTreeAddress(); - objlast = gDirectory->GetList()->Last(); - - // Initialise detector materials, geometry, histograms,etc - detector->CreateMaterials(); - detector->CreateGeometry(); - detector->BuildGeometry(); - detector->Init(); - - // Add Detector histograms in Detector list of histograms - if (objlast) objfirst = gDirectory->GetList()->After(objlast); - else objfirst = gDirectory->GetList()->First(); - while (objfirst) { - detector->Histograms()->Add(objfirst); - objfirst = gDirectory->GetList()->After(objfirst); - } - } - SetTransPar(); //Read the cuts for all materials - - MediaTable(); //Build the special IMEDIA table - - //Close the geometry structure - gMC->Ggclos(); - - //Initialise geometry deposition table - sEventEnergy.Set(gMC->NofVolumes()+1); - sSummEnergy.Set(gMC->NofVolumes()+1); - sSum2Energy.Set(gMC->NofVolumes()+1); - - //Create the color table - gMC->SetColors(); - - //Compute cross-sections - gMC->Gphysi(); - - //Write Geometry object to current file. - fGeometry->Write(); - - fInitDone = kTRUE; -} - -//_____________________________________________________________________________ -void AliRun::MediaTable() -{ - // - // Built media table to get from the media number to - // the detector id - // - Int_t kz, nz, idt, lz, i, k, ind; - // Int_t ibeg; - TObjArray &dets = *gAlice->Detectors(); - AliModule *det; + Int_t kz, nz, idt, lz, i, k, ind; + // Int_t ibeg; + TObjArray &dets = *gAlice->Detectors(); + AliModule *det; // // For all detectors for (kz=0;kzGetName(),generator->GetName()); + else + Warning("ResetGenerator","Replacing generator %s with NULL\n", + fGenerator->GetName()); + fGenerator = generator; +} + +//____________________________________________________________________________ +void AliRun::SetTransPar(char *filename) +{ + fTransParName = filename; +} + +//____________________________________________________________________________ +void AliRun::SetBaseFile(char *filename) +{ + fBaseFileName = filename; +} + +//____________________________________________________________________________ +void AliRun::ReadTransPar() { // // Read filename to set the transport parameters // - const Int_t ncuts=10; - const Int_t nflags=11; - const Int_t npars=ncuts+nflags; - const char pars[npars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO", + const Int_t kncuts=10; + const Int_t knflags=11; + const Int_t knpars=kncuts+knflags; + const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO", "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI", "BREM","COMP","DCAY","DRAY","HADR","LOSS", "MULS","PAIR","PHOT","RAYL"}; char line[256]; char detName[7]; char* filtmp; - Float_t cut[ncuts]; - Int_t flag[nflags]; + Float_t cut[kncuts]; + Int_t flag[knflags]; Int_t i, itmed, iret, ktmed, kz; FILE *lun; // // See whether the file is there - filtmp=gSystem->ExpandPathName(filename); + filtmp=gSystem->ExpandPathName(fTransParName.Data()); lun=fopen(filtmp,"r"); delete [] filtmp; if(!lun) { - Warning("SetTransPar","File %s does not exist!\n",filename); + Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data()); return; } // - printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n"); - printf(" *%59s\n","*"); - printf(" * Please check carefully what you are doing!%10s\n","*"); - printf(" *%59s\n","*"); + if(fDebug) { + printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n"); + printf(" *%59s\n","*"); + printf(" * Please check carefully what you are doing!%10s\n","*"); + printf(" *%59s\n","*"); + } // while(1) { // Initialise cuts and flags - for(i=0;iGetName()); + Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName()); continue; } // Set energy thresholds - for(kz=0;kz=0) { - printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n", - pars[kz],cut[kz],itmed,mod->GetName()); - gMC->Gstpar(ktmed,pars[kz],cut[kz]); + if(fDebug) printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n", + kpars[kz],cut[kz],itmed,mod->GetName()); + gMC->Gstpar(ktmed,kpars[kz],cut[kz]); } } // Set transport mechanisms - for(kz=0;kz=0) { - printf(" * %-6s set to %10d for tracking medium code %4d for %s\n", - pars[ncuts+kz],flag[kz],itmed,mod->GetName()); - gMC->Gstpar(ktmed,pars[ncuts+kz],Float_t(flag[kz])); + if(fDebug) printf(" * %-6s set to %10d for tracking medium code %4d for %s\n", + kpars[kncuts+kz],flag[kz],itmed,mod->GetName()); + gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz])); } } } else { - Warning("SetTransPar","Invalid medium code %d *\n",itmed); + Warning("ReadTransPar","Invalid medium code %d *\n",itmed); continue; } } else { - Warning("SetTransPar","Module %s not present\n",detName); + if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName); continue; } } } + //_____________________________________________________________________________ -void AliRun::MakeTree(Option_t *option) +void AliRun::MakeTree(Option_t *option, const char *file) { // // Create the ROOT trees // Loop on all detectors to create the Root branch (if any) // + char hname[30]; // // Analyse options - char *K = strstr(option,"K"); - char *H = strstr(option,"H"); - char *E = strstr(option,"E"); - char *D = strstr(option,"D"); - char *R = strstr(option,"R"); - // - if (K && !fTreeK) fTreeK = new TTree("TreeK0","Kinematics"); - if (H && !fTreeH) fTreeH = new TTree("TreeH0","Hits"); - if (D && !fTreeD) fTreeD = new TTree("TreeD0","Digits"); - if (E && !fTreeE) fTreeE = new TTree("TE","Header"); - if (R && !fTreeR) fTreeR = new TTree("TreeR0","Reconstruction"); - if (fTreeH) fTreeH->SetAutoSave(1000000000); //no autosave + const char *oK = strstr(option,"K"); + const char *oH = strstr(option,"H"); + const char *oTR = strstr(option,"T"); + const char *oE = strstr(option,"E"); + const char *oD = strstr(option,"D"); + const char *oR = strstr(option,"R"); + const char *oS = strstr(option,"S"); + // + + TDirectory *cwd = gDirectory; + + TBranch *branch = 0; + + if (oK) fStack->MakeTree(fEvent, file); + + if (oE && !fTreeE) { + fTreeE = new TTree("TE","Header"); + // branch = fTreeE->Branch("Header", "AliHeader", &fHeader, 4000, 0); + branch = fTreeE->Branch("Header", "AliHeader", &fHeader, 4000, 0); + branch->SetAutoDelete(kFALSE); + TFolder *folder = (TFolder *)gROOT->FindObjectAny("/Folders/RunMC/Event/Header"); + if (folder) folder->Add(fHeader); +// branch = fTreeE->Branch("Stack","AliStack", &fStack, 4000, 0); +// branch->SetAutoDelete(kFALSE); +// if (folder) folder->Add(fStack); + fTreeE->Write(0,TObject::kOverwrite); + } + + if (file && branch) { + char * outFile = new char[strlen(gAlice->GetBaseFile())+strlen(file)+2]; + sprintf(outFile,"%s/%s",GetBaseFile(),file); + branch->SetFile(outFile); + TIter next( branch->GetListOfBranches()); + while ((branch=(TBranch*)next())) { + branch->SetFile(outFile); + } + if (GetDebug()>1) + printf("* MakeBranch * Diverting Branch %s to file %s\n", branch->GetName(),file); + cwd->cd(); + delete outFile; + } + + if (oH && !fTreeH) { + sprintf(hname,"TreeH%d",fEvent); + fTreeH = new TTree(hname,"Hits"); + fTreeH->SetAutoSave(1000000000); //no autosave + fTreeH->Write(0,TObject::kOverwrite); + } + + if (oTR && !fTreeTR) { + sprintf(hname,"TreeTR%d",fEvent); + fTreeTR = new TTree(hname,"TrackReferences"); + fTreeTR->SetAutoSave(1000000000); //no autosave + fTreeTR->Write(0,TObject::kOverwrite); + } + + if (oD && !fTreeD) { + sprintf(hname,"TreeD%d",fEvent); + fTreeD = new TTree(hname,"Digits"); + fTreeD->Write(0,TObject::kOverwrite); + } + if (oS && !fTreeS) { + sprintf(hname,"TreeS%d",fEvent); + fTreeS = new TTree(hname,"SDigits"); + fTreeS->Write(0,TObject::kOverwrite); + } + if (oR && !fTreeR) { + sprintf(hname,"TreeR%d",fEvent); + fTreeR = new TTree(hname,"Reconstruction"); + fTreeR->Write(0,TObject::kOverwrite); + } + // // Create a branch for hits/digits for each detector // Each branch is a TClonesArray. Each data member of the Hits classes @@ -1054,162 +1271,58 @@ void AliRun::MakeTree(Option_t *option) TIter next(fModules); AliModule *detector; while((detector = (AliModule*)next())) { - if (H || D || R) detector->MakeBranch(option); + if (oH) detector->MakeBranch(option,file); + if (oTR) detector->MakeBranchTR(option,file); } - // Create a branch for particles - if (fTreeK && K) fTreeK->Branch("Particles",&fParticles,4000); - - // Create a branch for Header - if (fTreeE && E) fTreeE->Branch("Header","AliHeader",&header,4000); } //_____________________________________________________________________________ -Int_t AliRun::PurifyKine(Int_t lastSavedTrack, Int_t nofTracks) +TParticle* AliRun::Particle(Int_t i) { - // - // PurifyKine with external parameters - // - fHgwmk = lastSavedTrack; - fNtrack = nofTracks; - PurifyKine(); - return fHgwmk; + return fStack->Particle(i); } //_____________________________________________________________________________ -void AliRun::PurifyKine() +void AliRun::ResetDigits() { // - // Compress kinematic tree keeping only flagged particles - // and renaming the particle id's in all the hits - // - TClonesArray &particles = *fParticles; - int nkeep=fHgwmk+1, parent, i; - TParticle *part, *partnew, *father; - AliHit *OneHit; - int *map = new int[particles.GetEntries()]; - - // Save in Header total number of tracks before compression - fHeader.SetNtrack(fHeader.GetNtrack()+fNtrack-fHgwmk); - - // Preset map, to be removed later - for(i=0; iTestBit(Keep_Bit)) { - - // This particle has to be kept - map[i]=nkeep; - if(i!=nkeep) { - - // Old and new are different, have to copy - partnew = (TParticle *)particles.UncheckedAt(nkeep); - *partnew = *part; - } else partnew = part; - - // as the parent is always *before*, it must be already - // in place. This is what we are checking anyway! - if((parent=partnew->GetFirstMother())>fHgwmk) { - if(map[parent]==-99) printf("map[%d] = -99!\n",parent); - partnew->SetFirstMother(map[parent]); - } - nkeep++; - } - } - fNtrack=nkeep; - - // Fix daughters information - for (i=0; iGetFirstMother(); - if(parent>=0) { - father = (TParticle *)particles.UncheckedAt(parent); - if(father->TestBit(Daughters_Bit)) { - - if(iGetFirstDaughter()) father->SetFirstDaughter(i); - if(i>father->GetLastDaughter()) father->SetLastDaughter(i); - } else { - // Iitialise daughters info for first pass - father->SetFirstDaughter(i); - father->SetLastDaughter(i); - father->SetBit(Daughters_Bit); - } - } - } - - // Now loop on all detectors and reset the hits + // Reset all Detectors digits + // TIter next(fModules); AliModule *detector; while((detector = (AliModule*)next())) { - if (!detector->Hits()) continue; - TClonesArray &vHits=*(detector->Hits()); - if(vHits.GetEntries() != detector->GetNhits()) - printf("vHits.GetEntries()!=detector->GetNhits(): %d != %d\n", - vHits.GetEntries(),detector->GetNhits()); - for (i=0; iGetNhits(); i++) { - OneHit = (AliHit *)vHits.UncheckedAt(i); - OneHit->SetTrack(map[OneHit->GetTrack()]); - } + detector->ResetDigits(); } - - fHgwmk=nkeep-1; - particles.SetLast(fHgwmk); - delete [] map; } //_____________________________________________________________________________ -void AliRun::Reset(Int_t run, Int_t idevent) +void AliRun::ResetSDigits() { // - // Reset all Detectors & kinematics & trees - // - char hname[30]; + // Reset all Detectors digits // - ResetStack(); - ResetHits(); - ResetDigits(); - - // Initialise event header - fHeader.Reset(run,idevent); - - if(fTreeK) { - fTreeK->Reset(); - sprintf(hname,"TreeK%d",idevent); - fTreeK->SetName(hname); - } - if(fTreeH) { - fTreeH->Reset(); - sprintf(hname,"TreeH%d",idevent); - fTreeH->SetName(hname); - } - if(fTreeD) { - fTreeD->Reset(); - sprintf(hname,"TreeD%d",idevent); - fTreeD->SetName(hname); - } - if(fTreeR) { - fTreeR->Reset(); - sprintf(hname,"TreeR%d",idevent); - fTreeR->SetName(hname); + TIter next(fModules); + AliModule *detector; + while((detector = (AliModule*)next())) { + detector->ResetSDigits(); } } //_____________________________________________________________________________ -void AliRun::ResetDigits() +void AliRun::ResetHits() { // - // Reset all Detectors digits + // Reset all Detectors hits // TIter next(fModules); AliModule *detector; while((detector = (AliModule*)next())) { - detector->ResetDigits(); + detector->ResetHits(); } } //_____________________________________________________________________________ -void AliRun::ResetHits() +void AliRun::ResetTrackReferences() { // // Reset all Detectors hits @@ -1217,7 +1330,7 @@ void AliRun::ResetHits() TIter next(fModules); AliModule *detector; while((detector = (AliModule*)next())) { - detector->ResetHits(); + detector->ResetTrackReferences(); } } @@ -1235,7 +1348,78 @@ void AliRun::ResetPoints() } //_____________________________________________________________________________ -void AliRun::Run(Int_t nevent, const char *setup) +void AliRun::InitMC(const char *setup) +{ + // + // Initialize the Alice setup + // + + if(fInitDone) { + Warning("Init","Cannot initialise AliRun twice!\n"); + return; + } + + gROOT->LoadMacro(setup); + gInterpreter->ProcessLine(fConfigFunction.Data()); + + // Register MC in configuration + AliConfig::Instance()->Add(gMC); + gMC->SetStack(fStack); + + gMC->DefineParticles(); //Create standard MC particles + + TObject *objfirst, *objlast; + + fNdets = fModules->GetLast()+1; + + // + //=================Create Materials and geometry + gMC->Init(); + + // Added also after in case of interactive initialisation of modules + fNdets = fModules->GetLast()+1; + + TIter next(fModules); + AliModule *detector; + while((detector = (AliModule*)next())) { + detector->SetTreeAddress(); + objlast = gDirectory->GetList()->Last(); + + // Add Detector histograms in Detector list of histograms + if (objlast) objfirst = gDirectory->GetList()->After(objlast); + else objfirst = gDirectory->GetList()->First(); + while (objfirst) { + detector->Histograms()->Add(objfirst); + objfirst = gDirectory->GetList()->After(objfirst); + } + } + ReadTransPar(); //Read the cuts for all materials + + MediaTable(); //Build the special IMEDIA table + + //Initialise geometry deposition table + fEventEnergy.Set(gMC->NofVolumes()+1); + fSummEnergy.Set(gMC->NofVolumes()+1); + fSum2Energy.Set(gMC->NofVolumes()+1); + + //Compute cross-sections + gMC->BuildPhysics(); + + //Write Geometry object to current file. + fGeometry->Write(); + + fInitDone = kTRUE; + + fMCQA = new AliMCQA(fNdets); + + AliConfig::Instance(); + // + // Save stuff at the beginning of the file to avoid file corruption + Write(); +} + +//_____________________________________________________________________________ +void AliRun::RunMC(Int_t nevent, const char *setup) { // // Main function to be called to process a galice run @@ -1244,35 +1428,137 @@ void AliRun::Run(Int_t nevent, const char *setup) // a positive number of events will cause the finish routine // to be called // - - Int_t i, todo; + fEventsPerRun = nevent; // check if initialisation has been done - if (!fInitDone) Init(setup); + if (!fInitDone) InitMC(setup); // Create the Root Tree with one branch per detector - if(!fEvent) { - gAlice->MakeTree("KHDER"); - } - todo = TMath::Abs(nevent); - for (i=0; iReset(fRun, fEvent); - gMC->Gtrigi(); - gMC->Gtrigc(); - gMC->Gtrig(); - gAlice->FinishEvent(); - fEvent++; + MakeTree("ESDRT"); + + if (gSystem->Getenv("CONFIG_SPLIT_FILE")) { + MakeTree("K","Kine.root"); + MakeTree("H","Hits.root"); + } else { + MakeTree("KH"); } - + + gMC->ProcessRun(nevent); + // End of this run, close files - if(nevent>0) gAlice->FinishRun(); + if(nevent>0) FinishRun(); } //_____________________________________________________________________________ -void AliRun::RunLego(const char *setup,Int_t ntheta,Float_t themin, - Float_t themax,Int_t nphi,Float_t phimin,Float_t phimax, - Float_t rmin,Float_t rmax,Float_t zmax) +void AliRun::RunReco(const char *selected, Int_t first, Int_t last) +{ + // + // Main function to be called to reconstruct Alice event + // + cout << "Found "<< gAlice->TreeE()->GetEntries() << "events" << endl; + Int_t nFirst = first; + Int_t nLast = (last < 0)? (Int_t) gAlice->TreeE()->GetEntries() : last; + + for (Int_t nevent = nFirst; nevent <= nLast; nevent++) { + cout << "Processing event "<< nevent << endl; + GetEvent(nevent); + // MakeTree("R"); + Digits2Reco(selected); + } +} + +//_____________________________________________________________________________ + +void AliRun::Hits2Digits(const char *selected) +{ + + // Convert Hits to sumable digits + // + for (Int_t nevent=0; neventTreeE()->GetEntries(); nevent++) { + GetEvent(nevent); + // MakeTree("D"); + Hits2SDigits(selected); + SDigits2Digits(selected); + } +} + + +//_____________________________________________________________________________ + +void AliRun::Tree2Tree(Option_t *option, const char *selected) +{ + // + // Function to transform the content of + // + // - TreeH to TreeS (option "S") + // - TreeS to TreeD (option "D") + // - TreeD to TreeR (option "R") + // + // If multiple options are specified ("SDR"), transformation will be done in sequence for + // selected detector and for all detectors if none is selected (detector string + // can contain blank separated list of detector names). + + + const char *oS = strstr(option,"S"); + const char *oD = strstr(option,"D"); + const char *oR = strstr(option,"R"); + + TObjArray *detectors = Detectors(); + + TIter next(detectors); + + AliDetector *detector = 0; + + TDirectory *cwd = gDirectory; + + char outFile[32]; + + while((detector = (AliDetector*)next())) { + if (selected) + if (strcmp(detector->GetName(),selected)) continue; + if (detector->IsActive()){ + if (gSystem->Getenv("CONFIG_SPLIT_FILE")) { + if (oS) { + sprintf(outFile,"SDigits.%s.root",detector->GetName()); + detector->MakeBranch("S",outFile); + } + if (oD) { + sprintf(outFile,"Digits.%s.root",detector->GetName()); + detector->MakeBranch("D",outFile); + } + if (oR) { + sprintf(outFile,"Reco.%s.root",detector->GetName()); + detector->MakeBranch("R",outFile); + } + } else { + detector->MakeBranch(option); + } + + cwd->cd(); + + if (oS) { + cout << "Hits2SDigits: Processing " << detector->GetName() << "..." << endl; + detector->Hits2SDigits(); + } + if (oD) { + cout << "SDigits2Digits: Processing " << detector->GetName() << "..." << endl; + detector->SDigits2Digits(); + } + if (oR) { + cout << "Digits2Reco: Processing " << detector->GetName() << "..." << endl; + detector->Digits2Reco(); + } + + cwd->cd(); + } + } +} + + +//_____________________________________________________________________________ +void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min, + Float_t c1max,Int_t nc2,Float_t c2min,Float_t c2max, + Float_t rmin,Float_t rmax,Float_t zmax, AliLegoGenerator* gener) { // // Generates lego plots of: @@ -1313,17 +1599,51 @@ void AliRun::RunLego(const char *setup,Int_t ntheta,Float_t themin, // // check if initialisation has been done - if (!fInitDone) Init(setup); + if (!fInitDone) InitMC(setup); + //Save current generator + AliGenerator *gen=Generator(); + + // Set new generator + if (!gener) gener = new AliLegoGenerator(); + ResetGenerator(gener); + // + // Configure Generator + gener->SetRadiusRange(rmin, rmax); + gener->SetZMax(zmax); + gener->SetCoor1Range(nc1, c1min, c1max); + gener->SetCoor2Range(nc2, c2min, c2max); - fLego = new AliLego("lego","lego"); - fLego->Init(ntheta,themin,themax,nphi,phimin,phimax,rmin,rmax,zmax); - fLego->Run(); + + //Create Lego object + fLego = new AliLego("lego",gener); + + //Prepare MC for Lego Run + gMC->InitLego(); + + //Run Lego Object + + //gMC->ProcessRun(nc1*nc2+1); + gMC->ProcessRun(nc1*nc2); // Create only the Root event Tree - gAlice->MakeTree("E"); + MakeTree("E"); // End of this run, close files - gAlice->FinishRun(); + FinishRun(); + // Restore current generator + ResetGenerator(gen); + // Delete Lego Object + delete fLego; fLego=0; +} + +//_____________________________________________________________________________ +void AliRun::SetConfigFunction(const char * config) +{ + // + // Set the signature of the function contained in Config.C to configure + // the run + // + fConfigFunction=config; } //_____________________________________________________________________________ @@ -1332,516 +1652,657 @@ void AliRun::SetCurrentTrack(Int_t track) // // Set current track number // - fCurrent = track; + fStack->SetCurrentTrack(track); } //_____________________________________________________________________________ void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, Float_t *vpos, Float_t *polar, Float_t tof, - const char *mecha, Int_t &ntr, Float_t weight) + AliMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) { +// Delegate to stack +// + + fStack->SetTrack(done, parent, pdg, pmom, vpos, polar, tof, + mech, ntr, weight, is); +} + +//_____________________________________________________________________________ +void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg, + Double_t px, Double_t py, Double_t pz, Double_t e, + Double_t vx, Double_t vy, Double_t vz, Double_t tof, + Double_t polx, Double_t poly, Double_t polz, + AliMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) +{ + // Delegate to stack // - // Load a track on the stack - // - // done 0 if the track has to be transported - // 1 if not - // parent identifier of the parent track. -1 for a primary - // pdg particle code - // pmom momentum GeV/c - // vpos position - // polar polarisation - // tof time of flight in seconds - // mecha production mechanism - // ntr on output the number of the track stored - // - TClonesArray &particles = *fParticles; - TParticle *particle; - Float_t mass; - const Int_t firstdaughter=-1; - const Int_t lastdaughter=-1; - const Int_t KS=0; - // const Float_t tlife=0; - - // - // Here we get the static mass - // For MC is ok, but a more sophisticated method could be necessary - // if the calculated mass is required - // also, this method is potentially dangerous if the mass - // used in the MC is not the same of the PDG database - // - mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); - Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+ - pmom[1]*pmom[1]+pmom[2]*pmom[2]); - - //printf("Loading particle %s mass %f ene %f No %d ip %d pos %f %f %f mom %f %f %f KS %d m %s\n", - //pname,mass,e,fNtrack,pdg,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],KS,mecha); - - particle=new(particles[fNtrack]) TParticle(pdg,KS,parent,-1,firstdaughter, - lastdaughter,pmom[0],pmom[1],pmom[2], - e,vpos[0],vpos[1],vpos[2],tof); - // polar[0],polar[1],polar[2],tof, - // mecha,weight); - ((TParticle*)particles[fNtrack])->SetPolarisation(TVector3(polar[0],polar[1],polar[2])); - ((TParticle*)particles[fNtrack])->SetWeight(weight); - if(!done) particle->SetBit(Done_Bit); - - if(parent>=0) { - particle=(TParticle*) fParticles->UncheckedAt(parent); - particle->SetLastDaughter(fNtrack); - if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); - } else { - // - // This is a primary track. Set high water mark for this event - fHgwmk=fNtrack; + fStack->SetTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof, + polx, poly, polz, mech, ntr, weight, is); + +} + +//_____________________________________________________________________________ +void AliRun::SetHighWaterMark(const Int_t nt) +{ // - // Set also number if primary tracks - fHeader.SetNprimary(fHgwmk+1); - fHeader.SetNtrack(fHgwmk+1); - } - ntr = fNtrack++; + // Set high water mark for last track in event + fStack->SetHighWaterMark(nt); } //_____________________________________________________________________________ void AliRun::KeepTrack(const Int_t track) { // - // flags a track to be kept + // Delegate to stack // - TClonesArray &particles = *fParticles; - ((TParticle*)particles[track])->SetBit(Keep_Bit); + fStack->KeepTrack(track); } +// +// MC Application +// + //_____________________________________________________________________________ -void AliRun::StepManager(Int_t id) const +void AliRun::ConstructGeometry() +{ + // + // Create modules, materials, geometry + // + + TStopwatch stw; + TIter next(fModules); + AliModule *detector; + printf("Geometry creation:\n"); + while((detector = (AliModule*)next())) { + stw.Start(); + // Initialise detector materials and geometry + detector->CreateMaterials(); + detector->CreateGeometry(); + printf("%10s R:%.2fs C:%.2fs\n", + detector->GetName(),stw.RealTime(),stw.CpuTime()); + } +} + +//_____________________________________________________________________________ +void AliRun::InitGeometry() +{ + // + // Initialize detectors and display geometry + // + + printf("Initialisation:\n"); + TStopwatch stw; + TIter next(fModules); + AliModule *detector; + while((detector = (AliModule*)next())) { + stw.Start(); + // Initialise detector and display geometry + detector->Init(); + detector->BuildGeometry(); + printf("%10s R:%.2fs C:%.2fs\n", + detector->GetName(),stw.RealTime(),stw.CpuTime()); + } + +} + +//_____________________________________________________________________________ +void AliRun::GeneratePrimaries() +{ + // + // Generate primary particles and fill them in the stack. + // + + Generator()->Generate(); +} + +//_____________________________________________________________________________ +void AliRun::BeginEvent() +{ + // Clean-up previous event + // Energy scores + fEventEnergy.Reset(); + // Clean detector information + CleanDetectors(); + // Reset stack info + fStack->Reset(); + + + // + // Reset all Detectors & kinematics & trees + // + char hname[30]; + // + // Initialise event header + fHeader->Reset(fRun,fEvent,fEventNrInRun); + // + fStack->BeginEvent(fEvent); + + // + if(fLego) { + fLego->BeginEvent(); + return; + } + + // + + ResetHits(); + ResetTrackReferences(); + ResetDigits(); + ResetSDigits(); + + + if(fTreeH) { + fTreeH->Reset(); + sprintf(hname,"TreeH%d",fEvent); + fTreeH->SetName(hname); + } + + if(fTreeTR) { + fTreeTR->Reset(); + sprintf(hname,"TreeTR%d",fEvent); + fTreeTR->SetName(hname); + } + + if(fTreeD) { + fTreeD->Reset(); + sprintf(hname,"TreeD%d",fEvent); + fTreeD->SetName(hname); + fTreeD->Write(0,TObject::kOverwrite); + } + if(fTreeS) { + fTreeS->Reset(); + sprintf(hname,"TreeS%d",fEvent); + fTreeS->SetName(hname); + fTreeS->Write(0,TObject::kOverwrite); + } + if(fTreeR) { + fTreeR->Reset(); + sprintf(hname,"TreeR%d",fEvent); + fTreeR->SetName(hname); + fTreeR->Write(0,TObject::kOverwrite); + } +} + +//_____________________________________________________________________________ +void AliRun::BeginPrimary() +{ + // + // Called at the beginning of each primary track + // + + // Reset Hits info + gAlice->ResetHits(); + gAlice->ResetTrackReferences(); + +} + +//_____________________________________________________________________________ +void AliRun::PreTrack() +{ + TObjArray &dets = *fModules; + AliModule *module; + + for(Int_t i=0; i<=fNdets; i++) + if((module = (AliModule*)dets[i])) + module->PreTrack(); + + fMCQA->PreTrack(); +} + +//_____________________________________________________________________________ +void AliRun::Stepping() { // // Called at every step during transport // - Int_t copy; + Int_t id = DetFromMate(gMC->GetMedium()); + if (id < 0) return; + // // --- If lego option, do it and leave - if (fLego) { + if (fLego) fLego->StepManager(); - return; - } - //Update energy deposition tables - sEventEnergy[gMC->CurrentVolID(copy)]+=gMC->Edep(); + else { + Int_t copy; + //Update energy deposition tables + AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep()); - //Call the appropriate stepping routine; - AliModule *det = (AliModule*)fModules->At(id); - if(det) det->StepManager(); + //Call the appropriate stepping routine; + AliModule *det = (AliModule*)fModules->At(id); + if(det && det->StepManagerIsEnabled()) { + fMCQA->StepManager(id); + det->StepManager(); + } + } } //_____________________________________________________________________________ -void AliRun::ReadEuclid(const char* filnam, const AliModule *det, char* topvol) +void AliRun::PostTrack() +{ + TObjArray &dets = *fModules; + AliModule *module; + + for(Int_t i=0; i<=fNdets; i++) + if((module = (AliModule*)dets[i])) + module->PostTrack(); +} + +//_____________________________________________________________________________ +void AliRun::FinishPrimary() { - // - // read in the geometry of the detector in euclid file format - // - // id_det : the detector identification (2=its,...) - // topvol : return parameter describing the name of the top - // volume of geometry. - // - // author : m. maire - // - // 28.07.98 - // several changes have been made by miroslav helbich - // subroutine is rewrited to follow the new established way of memory - // booking for tracking medias and rotation matrices. - // all used tracking media have to be defined first, for this you can use - // subroutine greutmed. - // top volume is searched as only volume not positioned into another - // - - Int_t i, nvol, iret, itmed, irot, numed, npar, ndiv, iaxe; - Int_t ndvmx, nr, flag; - char key[5], card[77], natmed[21]; - char name[5], mother[5], shape[5], konly[5], volst[7000][5]; - char *filtmp; - Float_t par[50]; - Float_t teta1, phi1, teta2, phi2, teta3, phi3, orig, step; - Float_t xo, yo, zo; - const Int_t maxrot=5000; - Int_t idrot[maxrot],istop[7000]; - FILE *lun; // - // *** The input filnam name will be with extension '.euc' - filtmp=gSystem->ExpandPathName(filnam); - lun=fopen(filtmp,"r"); - delete [] filtmp; - if(!lun) { - Error("ReadEuclid","Could not open file %s\n",filnam); - return; + // Called at the end of each primary track + // + + // static Int_t count=0; + // const Int_t times=10; + // This primary is finished, purify stack + fStack->PurifyKine(); + + TIter next(fModules); + AliModule *detector; + while((detector = (AliModule*)next())) { + detector->FinishPrimary(); } - //* --- definition of rotation matrix 0 --- - TArrayI &idtmed = *(det->GetIdtmed()); - for(i=1; i=100 ) { - Error("ReadEuclid","TMED illegal medium number %d for %s\n",itmed,natmed); - exit(1); - } - //Pad the string with blanks - i=-1; - while(natmed[++i]); - while(i<20) natmed[i++]=' '; - natmed[i]='\0'; - // - if( idtmed[itmed]<=0 ) { - Error("ReadEuclid","TMED undefined medium number %d for %s\n",itmed,natmed); - exit(1); - } - gMC->Gckmat(idtmed[itmed],natmed); - //* - } else if (!strcmp(key,"ROTM")) { - sscanf(&card[4],"%d %f %f %f %f %f %f",&irot,&teta1,&phi1,&teta2,&phi2,&teta3,&phi3); - if( irot<=0 || irot>=maxrot ) { - Error("ReadEuclid","ROTM rotation matrix number %d illegal\n",irot); - exit(1); - } - det->AliMatrix(idrot[irot],teta1,phi1,teta2,phi2,teta3,phi3); - //* - } else if (!strcmp(key,"VOLU")) { - sscanf(&card[5],"'%[^']' '%[^']' %d %d", name, shape, &numed, &npar); - if (npar>0) { - for(i=0;iGsvolu( name, shape, idtmed[numed], par, npar); - //* save the defined volumes - strcpy(volst[++nvol],name); - istop[nvol]=1; - //* - } else if (!strcmp(key,"DIVN")) { - sscanf(&card[5],"'%[^']' '%[^']' %d %d", name, mother, &ndiv, &iaxe); - gMC->Gsdvn ( name, mother, ndiv, iaxe ); - //* - } else if (!strcmp(key,"DVN2")) { - sscanf(&card[5],"'%[^']' '%[^']' %d %d %f %d",name, mother, &ndiv, &iaxe, &orig, &numed); - gMC->Gsdvn2( name, mother, ndiv, iaxe, orig,idtmed[numed]); - //* - } else if (!strcmp(key,"DIVT")) { - sscanf(&card[5],"'%[^']' '%[^']' %f %d %d %d", name, mother, &step, &iaxe, &numed, &ndvmx); - gMC->Gsdvt ( name, mother, step, iaxe, idtmed[numed], ndvmx); - //* - } else if (!strcmp(key,"DVT2")) { - sscanf(&card[5],"'%[^']' '%[^']' %f %d %f %d %d", name, mother, &step, &iaxe, &orig, &numed, &ndvmx); - gMC->Gsdvt2 ( name, mother, step, iaxe, orig, idtmed[numed], ndvmx ); - //* - } else if (!strcmp(key,"POSI")) { - sscanf(&card[5],"'%[^']' %d '%[^']' %f %f %f %d '%[^']'", name, &nr, mother, &xo, &yo, &zo, &irot, konly); - if( irot<0 || irot>=maxrot ) { - Error("ReadEuclid","POSI %s#%d rotation matrix number %d illegal\n",name,nr,irot); - exit(1); - } - if( idrot[irot] == -99) { - Error("ReadEuclid","POSI %s#%d undefined matrix number %d\n",name,nr,irot); - exit(1); - } - //*** volume name cannot be the top volume - for(i=1;i<=nvol;i++) { - if (!strcmp(volst[i],name)) istop[i]=0; - } - //* - gMC->Gspos ( name, nr, mother, xo, yo, zo, idrot[irot], konly ); - //* - } else if (!strcmp(key,"POSP")) { - sscanf(&card[5],"'%[^']' %d '%[^']' %f %f %f %d '%[^']' %d", name, &nr, mother, &xo, &yo, &zo, &irot, konly, &npar); - if( irot<0 || irot>=maxrot ) { - Error("ReadEuclid","POSP %s#%d rotation matrix number %d illegal\n",name,nr,irot); - exit(1); - } - if( idrot[irot] == -99) { - Error("ReadEuclid","POSP %s#%d undefined matrix number %d\n",name,nr,irot); - exit(1); - } - if (npar > 0) { - for(i=0;iGsposp ( name, nr, mother, xo,yo,zo, idrot[irot], konly, par, npar); - } - //* - if (strcmp(key,"END")) goto L10; - //* find top volume in the geometry - flag=0; - for(i=1;i<=nvol;i++) { - if (istop[i] && flag) { - Warning("ReadEuclid"," %s is another possible top volume\n",volst[i]); - } - if (istop[i] && !flag) { - strcpy(topvol,volst[i]); - printf(" *** GREUCL *** volume %s taken as a top volume\n",topvol); - flag=1; - } + + // Write out hits if any + if (gAlice->TreeH()) { + gAlice->TreeH()->Fill(); } - if (!flag) { - Warning("ReadEuclid","top volume not found\n"); + + // Write out hits if any + if (gAlice->TreeTR()) { + gAlice->TreeTR()->Fill(); } - fclose (lun); - //* - //* commented out only for the not cernlib version - printf(" *** GREUCL *** file: %s is now read in\n",filnam); + // - return; - //* - L20: - Error("ReadEuclid","reading error or premature end of file\n"); + // if(++count%times==1) gObjectTable->Print(); } //_____________________________________________________________________________ -void AliRun::ReadEuclidMedia(const char* filnam, const AliModule *det) +void AliRun::FinishEvent() { - // - // read in the materials and tracking media for the detector - // in euclid file format - // - // filnam: name of the input file - // id_det: id_det is the detector identification (2=its,...) - // - // author : miroslav helbich - // - Float_t sxmgmx = gAlice->Field()->Max(); - Int_t isxfld = gAlice->Field()->Integ(); - Int_t end, i, iret, itmed; - char key[5], card[130], natmed[21], namate[21]; - Float_t ubuf[50]; - char* filtmp; - FILE *lun; - Int_t imate; - Int_t nwbuf, isvol, ifield, nmat; - Float_t a, z, dens, radl, absl, fieldm, tmaxfd, stemax, deemax, epsil, stmin; // - end=strlen(filnam); - for(i=0;iExpandPathName(filnam); - lun=fopen(filtmp,"r"); - delete [] filtmp; - if(!lun) { - Warning("ReadEuclidMedia","Could not open file %s\n",filnam); - return; - } + // - // Retrieve Mag Field parameters - Int_t ISXFLD=gAlice->Field()->Integ(); - Float_t SXMGMX=gAlice->Field()->Max(); - // TArrayI &idtmed = *(det->GetIdtmed()); - // - L10: - for(i=0;i<130;i++) card[i]=0; - iret=fscanf(lun,"%4s %[^\n]",key,card); - if(iret<=0) goto L20; - fscanf(lun,"%*c"); - //* - //* read material - if (!strcmp(key,"MATE")) { - sscanf(card,"%d '%[^']' %f %f %f %f %f %d",&imate,namate,&a,&z,&dens,&radl,&absl,&nwbuf); - if (nwbuf>0) for(i=0;iAliMaterial(imate,namate,a,z,dens,radl,absl,ubuf,nwbuf); - //* read tracking medium - } else if (!strcmp(key,"TMED")) { - sscanf(card,"%d '%[^']' %d %d %d %f %f %f %f %f %f %d", - &itmed,natmed,&nmat,&isvol,&ifield,&fieldm,&tmaxfd, - &stemax,&deemax,&epsil,&stmin,&nwbuf); - if (nwbuf>0) for(i=0;iAliMedium(itmed,natmed,nmat,isvol,ISXFLD,SXMGMX,tmaxfd, - stemax,deemax,epsil,stmin,ubuf,nwbuf); - // (*fImedia)[idtmed[itmed]-1]=id_det; - //* - } - //* - if (strcmp(key,"END")) goto L10; - fclose (lun); - //* - //* commented out only for the not cernlib version - Warning("ReadEuclidMedia","file: %s is now read in\n",filnam); - //* - return; - //* - L20: - Warning("ReadEuclidMedia","reading error or premature end of file\n"); -} - + if(fLego) fLego->FinishEvent(); + + //Update the energy deposit tables + Int_t i; + for(i=0;iSetNprimary(fStack->GetNprimary()); + fHeader->SetNtrack(fStack->GetNtrack()); + + + // Write out the kinematics + fStack->FinishEvent(); + + // Write out the event Header information + if (fTreeE) { + fHeader->SetStack(fStack); + fTreeE->Fill(); + } + + + // Write Tree headers + TTree* pTreeK = fStack->TreeK(); + if (pTreeK) pTreeK->Write(0,TObject::kOverwrite); + if (fTreeH) fTreeH->Write(0,TObject::kOverwrite); + if (fTreeTR) fTreeTR->Write(0,TObject::kOverwrite); + + ++fEvent; + ++fEventNrInRun; +} + +//_____________________________________________________________________________ +void AliRun::Field(const Double_t* x, Double_t *b) const +{ + Float_t xfloat[3]; + for (Int_t i=0; i<3; i++) xfloat[i] = x[i]; + + if (Field()) { + Float_t bfloat[3]; + Field()->Field(xfloat,bfloat); + for (Int_t j=0; j<3; j++) b[j] = bfloat[j]; + } + else { + printf("No mag field defined!\n"); + b[0]=b[1]=b[2]=0.; + } + +} + +// +// End of MC Application +// + //_____________________________________________________________________________ void AliRun::Streamer(TBuffer &R__b) { - // // Stream an object of class AliRun. - // + if (R__b.IsReading()) { - Version_t R__v = R__b.ReadVersion(); if (R__v) { } - TNamed::Streamer(R__b); if (!gAlice) gAlice = this; + + AliRun::Class()->ReadBuffer(R__b, this); + // gROOT->GetListOfBrowsables()->Add(this,"Run"); + fTreeE = (TTree*)gDirectory->Get("TE"); - if (fTreeE) fTreeE->SetBranchAddress("Header", &header); + if (fTreeE) { + fTreeE->SetBranchAddress("Header", &fHeader); + } else Error("Streamer","cannot find Header Tree\n"); - R__b >> fNtrack; - R__b >> fHgwmk; - R__b >> fDebug; - fHeader.Streamer(R__b); - R__b >> fModules; - R__b >> fParticles; - R__b >> fField; - // R__b >> fMC; - R__b >> fNdets; - R__b >> fTrRmax; - R__b >> fTrZmax; - R__b >> fGenerator; - if(R__v>1) { - R__b >> fPDGDB; //Particle factory object! - fTreeE->GetEntry(0); - } else { - fHeader.SetEvent(0); - fPDGDB = TDatabasePDG::Instance(); //Particle factory object! - } + + fTreeE->GetEntry(0); + gRandom = fRandom; } else { - R__b.WriteVersion(AliRun::IsA()); - TNamed::Streamer(R__b); - R__b << fNtrack; - R__b << fHgwmk; - R__b << fDebug; - fHeader.Streamer(R__b); - R__b << fModules; - R__b << fParticles; - R__b << fField; - // R__b << fMC; - R__b << fNdets; - R__b << fTrRmax; - R__b << fTrZmax; - R__b << fGenerator; - R__b << fPDGDB; //Particle factory object! - } -} + AliRun::Class()->WriteBuffer(R__b, this); + } +} -//_____________________________________________________________________________ -// -// Interfaces to Fortran -// -//_____________________________________________________________________________ +//___________________________________________________________________________ +Int_t AliRun::CurrentTrack() const { + // + // Returns current track + // + return fStack->CurrentTrack(); +} -extern "C" void type_of_call rxgtrak (Int_t &mtrack, Int_t &ipart, Float_t *pmom, - Float_t &e, Float_t *vpos, Float_t *polar, - Float_t &tof) -{ +//___________________________________________________________________________ +Int_t AliRun::GetNtrack() const { + // + // Returns number of tracks in stack + // + return fStack->GetNtrack(); +} + +//___________________________________________________________________________ +TObjArray* AliRun::Particles() { // - // Fetches next track from the ROOT stack for transport. Called by the - // modified version of GTREVE. + // Returns pointer to Particles array + // + return fStack->Particles(); +} + +//___________________________________________________________________________ +TTree* AliRun::TreeK() { // - // Track number in the ROOT stack. If MTRACK=0 no - // mtrack more tracks are left in the stack to be - // transported. - // ipart Particle code in the GEANT conventions. - // pmom[3] Particle momentum in GeV/c - // e Particle energy in GeV - // vpos[3] Particle position - // tof Particle time of flight in seconds + // Returns pointer to the TreeK array // - Int_t pdg; - gAlice->GetNextTrack(mtrack, pdg, pmom, e, vpos, polar, tof); - ipart = gMC->IdFromPDG(pdg); - mtrack++; + return fStack->TreeK(); } -//_____________________________________________________________________________ -extern "C" void type_of_call -#ifndef WIN32 -rxstrak (Int_t &keep, Int_t &parent, Int_t &ipart, Float_t *pmom, - Float_t *vpos, Float_t &tof, const char* cmech, Int_t &ntr, const int cmlen) -#else -rxstrak (Int_t &keep, Int_t &parent, Int_t &ipart, Float_t *pmom, - Float_t *vpos, Float_t &tof, const char* cmech, const int cmlen, - Int_t &ntr) -#endif + +//___________________________________________________________________________ +void AliRun::SetGenEventHeader(AliGenEventHeader* header) +{ + fHeader->SetGenEventHeader(header); +} + +//___________________________________________________________________________ +TFile* AliRun::InitFile(TString fileName) +{ +// +// create the file where the whole tree will be saved +// + TDirectory *wd = gDirectory; + TFile* file = TFile::Open(fileName,"update"); + gDirectory = wd; + if (!file->IsOpen()) { + Error("Cannot open file, %s\n",fileName); + return 0; + } + return file; +} + +//___________________________________________________________________________ +TFile* AliRun::InitTreeFile(Option_t *option, TString fileName) { // - // Fetches next track from the ROOT stack for transport. Called by GUKINE - // and GUSTEP. - // - // Status of the track. If keep=0 the track is put - // keep on the ROOT stack but it is not fetched for - // transport. - // parent Parent track. If parent=0 the track is a primary. - // In GUSTEP the routine is normally called to store - // secondaries generated by the current track whose - // ROOT stack number is MTRACK (common SCKINE. - // ipart Particle code in the GEANT conventions. - // pmom[3] Particle momentum in GeV/c - // vpos[3] Particle position - // tof Particle time of flight in seconds - // - // cmech (CHARACTER*10) Particle origin. This field is user - // defined and it is not used inside the GALICE code. - // ntr Number assigned to the particle in the ROOT stack. - // - char mecha[11]; - Float_t polar[3]={0.,0.,0.}; - for(int i=0; i<10 && iPDGFromId(ipart); - gAlice->SetTrack(keep, parent-1, pdg, pmom, vpos, polar, tof, mecha, ntr); - ntr++; + // create the file where one of the following trees will be saved + // trees: S,D,R + // WARNING: by default these trees are saved on the file on which + // hits are stored. If you divert one of these trees, you cannot restore + // it to the original file (usually galice.root) in the same aliroot session + Bool_t oS = (strstr(option,"S")!=0); + Bool_t oR = (strstr(option,"R")!=0); + Bool_t oD = (strstr(option,"D")!=0); + Int_t choice[3]; + for (Int_t i=0; i<3; i++) choice[i] = 0; + if(oS)choice[0] = 1; + if(oD)choice[1] = 1; + if(oR)choice[2] = 1; + + TFile *ptr=0; + + if(!(oS || oR || oD))return ptr; + + Int_t active[3]; + for (Int_t i=0; i<3; i++) active[i] = 0; + if(fTreeSFileName != "") active[0] = 1; + if(fTreeDFileName != "") active[1] = 1; + if(fTreeDFileName != "") active[2] = 1; + + Bool_t alreadyopen1 = kFALSE; + Bool_t alreadyopen2 = kFALSE; + + if(oS){ + // if already active and same name with non-null ptr + if(active[0]==1 && fileName == fTreeSFileName && fTreeSFile){ + Warning("InitTreeFile","File %s already opened",fTreeSFileName.Data()); + ptr = fTreeSFile; + } + else { + // if already active with different name with non-null ptr + if(active[0]==1 && fileName != fTreeSFileName && fTreeSFile){ + // close the active files and also the other possible files in option + CloseTreeFile(option); + } + fTreeSFileName = fileName; + alreadyopen1 = + (active[1] == 1 && fTreeDFileName == fTreeSFileName && fTreeDFile); + alreadyopen2 = + (active[2] == 1 && fTreeRFileName == fTreeSFileName && fTreeRFile); + if(!(alreadyopen1 || alreadyopen2)){ + ptr = InitFile(fileName); + fTreeSFile = ptr; + } + else { + if(alreadyopen1){fTreeSFile = fTreeDFile; ptr = fTreeSFile;} + if(alreadyopen2){fTreeSFile = fTreeRFile; ptr = fTreeSFile;} + } + if(choice[1] == 1) { fTreeDFileName = fileName; fTreeDFile = ptr;} + if(choice[2] == 1) { fTreeRFileName = fileName; fTreeRFile = ptr;} + } + return ptr; + } + + if(oD){ + // if already active and same name with non-null ptr + if(active[1]==1 && fileName == fTreeDFileName && fTreeDFile){ + Warning("InitTreeFile","File %s already opened",fTreeDFileName.Data()); + ptr = fTreeDFile; + } + else { + // if already active with different name with non-null ptr + if(active[1]==1 && fileName != fTreeDFileName && fTreeDFile){ + // close the active files and also the other possible files in option + CloseTreeFile(option); + } + fTreeDFileName = fileName; + alreadyopen1 = + (active[0] == 1 && fTreeSFileName == fTreeDFileName && fTreeSFile); + alreadyopen2 = + (active[2] == 1 && fTreeRFileName == fTreeDFileName && fTreeRFile); + if(!(alreadyopen1 || alreadyopen2)){ + ptr = InitFile(fileName); + fTreeDFile = ptr; + } + else { + if(alreadyopen1){fTreeDFile = fTreeSFile; ptr = fTreeDFile;} + if(alreadyopen2){fTreeDFile = fTreeRFile; ptr = fTreeDFile;} + } + if(choice[2] == 1) { fTreeRFileName = fileName; fTreeRFile = ptr;} + } + return ptr; + } + + if(oR){ + // if already active and same name with non-null ptr + if(active[2]==1 && fileName == fTreeRFileName && fTreeRFile){ + Warning("InitTreeFile","File %s already opened",fTreeRFileName.Data()); + ptr = fTreeRFile; + } + else { + // if already active with different name with non-null ptr + if(active[2]==1 && fileName != fTreeRFileName && fTreeRFile){ + // close the active files and also the other possible files in option + CloseTreeFile(option); + } + fTreeRFileName = fileName; + alreadyopen1 = + (active[1] == 1 && fTreeDFileName == fTreeRFileName && fTreeDFile); + alreadyopen2 = + (active[0]== 1 && fTreeSFileName == fTreeRFileName && fTreeSFile); + if(!(alreadyopen1 || alreadyopen2)){ + ptr = InitFile(fileName); + fTreeRFile = ptr; + } + else { + if(alreadyopen1){fTreeRFile = fTreeDFile; ptr = fTreeRFile;} + if(alreadyopen2){fTreeRFile = fTreeSFile; ptr = fTreeRFile;} + } + } + return ptr; + } + } -//_____________________________________________________________________________ -extern "C" void type_of_call rxkeep(const Int_t &n) +//___________________________________________________________________________ +void AliRun::PrintTreeFile() { - if( NULL==gAlice ) exit(1); - - if( n<=0 || n>gAlice->Particles()->GetEntries() ) - { - printf(" Bad index n=%d must be 0Particles()->GetEntries()); - exit(1); + // + // prints the file names and pointer associated to S,D,R trees + // + cout<<"===================================================\n"; + TFile *file = fTreeE->GetCurrentFile(); + TString curfilname=""; + if(file)curfilname=(TString)file->GetName(); + cout<<" Current tree file name: "<IsOpen()){ + fTreeSFile->Close(); + delete fTreeSFile; + } } - - ((TParticle*)(gAlice->Particles()->UncheckedAt(n-1)))->SetBit(Keep_Bit); + fTreeSFile = 0; + } + if(oD){ + fTreeDFileName = ""; + if(fTreeDFile){ + if(!((fTreeDFile == fTreeRFile) || (fTreeDFile == fTreeSFile)) && + fTreeDFile->IsOpen()){ + fTreeDFile->Close(); + delete fTreeDFile; + } + } + fTreeDFile = 0; + } + if(oR){ + fTreeRFileName = ""; + if(fTreeRFile){ + if(!((fTreeRFile == fTreeSFile) || (fTreeRFile == fTreeDFile)) && + fTreeRFile->IsOpen()){ + fTreeRFile->Close(); + delete fTreeRFile; + } + } + fTreeRFile = 0; + } } -//_____________________________________________________________________________ -extern "C" void type_of_call rxouth () +//___________________________________________________________________________ +void AliRun::MakeTree(Option_t *option, TFile *file) { // - // Called by Gtreve at the end of each primary track + // Create some trees in the separate file // - gAlice->FinishPrimary(); -} + const char *oD = strstr(option,"D"); + const char *oR = strstr(option,"R"); + const char *oS = strstr(option,"S"); + + TDirectory *cwd = gDirectory; + char hname[30]; + if (oD) { + delete fTreeD; + sprintf(hname,"TreeD%d",fEvent); + file->cd(); + fTreeD = static_cast(file->Get("hname")); + if (!fTreeD) { + fTreeD = new TTree(hname,"Digits"); + fTreeD->Write(0,TObject::kOverwrite); + } + cwd->cd(); + } + if (oS) { + delete fTreeS; + sprintf(hname,"TreeS%d",fEvent); + file->cd(); + fTreeS = static_cast(file->Get("hname")); + if (!fTreeS) { + fTreeS = new TTree(hname,"SDigits"); + fTreeS->Write(0,TObject::kOverwrite); + } + cwd->cd(); + } + + if (oR) { + delete fTreeR; + sprintf(hname,"TreeR%d",fEvent); + file->cd(); + fTreeR = static_cast(file->Get("hname")); + if (!fTreeR) { + fTreeR = new TTree(hname,"RecPoint"); + fTreeR->Write(0,TObject::kOverwrite); + } + cwd->cd(); + } +}