[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderESD.cxx

#include "AliHBTReaderESD.h"

#include <TPDGCode.h>
#include <TString.h>
#include <TObjString.h>
#include <TTree.h>
#include <TFile.h>
#include <TParticle.h>

#include <AliRunLoader.h>
#include <AliStack.h>
#include <AliESDtrack.h>
#include <AliESD.h>

#include "AliHBTRun.h"
#include "AliHBTEvent.h"
#include "AliHBTParticle.h"
#include "AliHBTParticleCut.h"

ClassImp(AliHBTReaderESD)

AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
 fESDFileName(esdfilename),
 fGAlFileName(galfilename),
 fFile(0x0),
 fRunLoader(0x0),
 fReadParticles(kFALSE)
{
  //cosntructor
  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
    Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/
  
AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
 AliHBTReader(dirs), 
 fESDFileName(esdfilename),
 fGAlFileName(galfilename),
 fFile(0x0),
 fRunLoader(0x0),
 fReadParticles(kFALSE)
{
  //cosntructor
  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
    Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/

AliHBTReaderESD::~AliHBTReaderESD()
{
 //desctructor
  delete fFile;
  delete fRunLoader;
}
/**********************************************************/
Int_t AliHBTReaderESD::ReadNext()
{
//reads next event from fFile
  
  //****** Tentative particle type "concentrations"
  static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
  
  Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
  Double_t w[kNSpecies];
  Double_t mom[3];//momentum
  Double_t pos[3];//position
  
  if (AliHBTParticle::fgDebug)
    Info("ReadNext","Entered");
    
  TDatabasePDG* pdgdb = TDatabasePDG::Instance();
  if (pdgdb == 0x0)
   {
     Error("Next","Can not get PDG Database Instance.");
     return 1;
   }
  
  if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
  if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
  
  fParticlesEvent->Reset();
  fTracksEvent->Reset();
        
  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
   {
     if (fFile == 0x0)
      {
       fFile = OpenFile(fCurrentDir);//rl is opened here
       if (fFile == 0x0)
         {
           Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
           fCurrentDir++;
           continue;
         }
       fCurrentEvent = 0;
      } 
      
     //try to read
     TString esdname;
     esdname+=fCurrentEvent;
     AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
     if (esd == 0x0)
      {
        if (AliHBTParticle::fgDebug > 2 )
          {
            Info("ReadNext","Can not find ESD object named %s.",esdname.Data());
          }
        fCurrentDir++;
        delete fFile;//we have to assume there is no more ESD objects in the fFile
        delete fRunLoader;
        fFile = 0x0;
        continue;
      }
     AliStack* stack = 0x0;
     if (fReadParticles && fRunLoader)
      {
        fRunLoader->GetEvent(fCurrentEvent);
        stack = fRunLoader->Stack();
      }
     Info("ReadNext","Reading Event %d",fCurrentEvent);
  
     Int_t ntr = esd->GetNumberOfTracks();
     Info("ReadNext","Found %d tracks.",ntr);
     for (Int_t i = 0;i<ntr; i++)
      {
        AliESDtrack *esdtrack = esd->GetTrack(i);
        if (esdtrack == 0x0)
         {
           Error("Next","Can not get track %d", i);
           continue;
         }
        if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) 
         {
           if (AliHBTParticle::fgDebug > 2) 
             Info("ReadNext","Particle skipped: PID BIT is not set.");
           continue;
         }

        esdtrack->GetESDpid(pidtable);
        esdtrack->GetPxPyPz(mom);
        esdtrack->GetXYZ(pos);
        
        Double_t extx;
        Double_t extp[5];
        esdtrack->GetExternalParameters(extx,extp);
        
        Int_t charge = (extp[4] > 0)?-1:1;//if curvature=-charg/Pt is positive charge is negative
        
        //Particle from kinematics
        AliHBTParticle* particle = 0;
        Bool_t keeppart = kFALSE;
        if ( fReadParticles && stack  )
         {
           if (esdtrack->GetLabel() < 0) continue;//this is fake -  we are not able to match any track 
           TParticle *p = stack->Particle(esdtrack->GetLabel());
           if (p==0x0) 
            {
              Error("ReadNext","Can not find track with such label.");
              continue;
            }
//           if(p->GetPdgCode()<0) charge = -1;
           particle = new AliHBTParticle(*p,i);
           
         }
         
        //Here we apply Bayes' formula
        Double_t rc=0.;
        for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
        if (rc==0.0) 
         {
           if (AliHBTParticle::fgDebug > 2) 
             Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
           continue;
         }

        for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;

        if (AliHBTParticle::fgDebug > 4)
         { 
           Info("ReadNext","###########################################################################");
           Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
           Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
           TString msg("Pid list got from track:");
           for (Int_t s = 0;s<kNSpecies;s++)
            {
              msg+="\n    ";
              msg+=s;
              msg+="(";
              msg+=charge*GetSpeciesPdgCode((ESpecies)s);
              msg+="): ";
              msg+=w[s];
              msg+=" (";
              msg+=pidtable[s];
              msg+=")";
            }
           Info("ReadNext","%s",msg.Data());
         }//if (AliHBTParticle::fgDebug>4)

        for (Int_t s = 0; s<kNSpecies; s++)
         {
           Float_t pp = w[s];
           if (pp == 0.0) continue;

           Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
           if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type 

           Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
           Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track

           AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, 
                                                      mom[0], mom[1], mom[2], tEtot,
                                                      pos[0], pos[1], pos[2], 0.);
           //copy probabilitis of other species (if not zero)
           for (Int_t k = 0; k<kNSpecies; k++)
            {
              if (k == s) continue;
              if (w[k] == 0.0) continue;
              track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
            }
           if(Pass(track))//check if meets all criteria of any of our cuts
                          //if it does not delete it and take next good track
            { 
              delete track;
              continue;
            }
           
           fTracksEvent->AddParticle(track);
           if (particle) fParticlesEvent->AddParticle(particle);
           keeppart = kTRUE;
               
           if (AliHBTParticle::fgDebug > 4 )
            {
              Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
              track->Print();
            }
         }//for (Int_t s = 0; s<kNSpecies; s++)
         
        if (keeppart == kFALSE) delete particle;//particle was not stored in event
        
      }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
     
     Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
            fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(),
            fNEventsRead,fCurrentEvent,fCurrentDir);
      
     fCurrentEvent++;
     fNEventsRead++;
     delete esd;
     return 0;//success -> read one event
   }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
   
  return 1; //no more directories to read
}
/**********************************************************/
void AliHBTReaderESD::Rewind()
{
  delete fFile;
  fFile = 0;
  delete fRunLoader;
  fCurrentDir = 0;
  fNEventsRead = 0;
  fCurrentEvent++;
  
}
/**********************************************************/

TFile* AliHBTReaderESD::OpenFile(Int_t n)
{
//opens fFile with kine tree

 const TString& dirname = GetDirName(n);
 if (dirname == "")
  {
   Error("OpenFiles","Can not get directory name");
   return 0x0;
  }
 TString filename = dirname +"/"+ fESDFileName;
 TFile *ret = TFile::Open(filename.Data()); 

 if ( ret == 0x0)
  {
    Error("OpenFiles","Can't open fFile %s",filename.Data());
    return 0x0;
  }
 if (!ret->IsOpen())
  {
    Error("OpenFiles","Can't open fFile  %s",filename.Data());
    return 0x0;
  }
 
 if (fReadParticles)
  {
   fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
   if (fRunLoader == 0x0)
    {
      Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
      delete ret;
      return 0x0;
    }
   fRunLoader->LoadHeader();
   if (fRunLoader->LoadKinematics())
    {
      Error("Next","Error occured while loading kinematics.");
      delete fRunLoader;
      delete ret;
      return 0x0;
    }
  }
 return ret;
}
/**********************************************************/

Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
{
  //returns pdg code from the PID index
  //ask jura about charge
  switch (spec)
   {
     case kESDElectron:
       return kPositron;
       break;
     case kESDMuon:
       return kMuonPlus;
       break;
     case kESDPion:
       return kPiPlus;
       break;
     case kESDKaon:
       return kKPlus;
       break;
     case kESDProton:
       return kProton;
       break;
     default:
       ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
       break;
   }
  return 0;
}
Commit	Line	Data
7be13e68	1	#include "AliHBTReaderESD.h"
	2
	3	#include <TPDGCode.h>
	4	#include <TString.h>
	5	#include <TObjString.h>
	6	#include <TTree.h>
	7	#include <TFile.h>
	8	#include <TParticle.h>
	9
bed069a4	10	#include <AliRunLoader.h>
bed069a4	11	#include <AliStack.h>
7be13e68	12	#include <AliESDtrack.h>
	13	#include <AliESD.h>
	14
	15	#include "AliHBTRun.h"
	16	#include "AliHBTEvent.h"
	17	#include "AliHBTParticle.h"
	18	#include "AliHBTParticleCut.h"
	19
	20	ClassImp(AliHBTReaderESD)
	21
bed069a4	22	AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
7be13e68	23	fESDFileName(esdfilename),
bed069a4	24	fGAlFileName(galfilename),
	25	fFile(0x0),
	26	fRunLoader(0x0),
	27	fReadParticles(kFALSE)
7be13e68	28	{
	29	//cosntructor
	30	if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
	31	Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
	32	}
	33	/********************************************************************/
	34
bed069a4	35	AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
7be13e68	36	AliHBTReader(dirs),
7be13e68	37	fESDFileName(esdfilename),
bed069a4	38	fGAlFileName(galfilename),
	39	fFile(0x0),
	40	fRunLoader(0x0),
	41	fReadParticles(kFALSE)
7be13e68	42	{
	43	//cosntructor
	44	if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
	45	Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
	46	}
	47	/********************************************************************/
	48
	49	AliHBTReaderESD::~AliHBTReaderESD()
	50	{
	51	//desctructor
bed069a4	52	delete fFile;
bed069a4	53	delete fRunLoader;
7be13e68	54	}
bed069a4	55	/**********************************************************/
bed069a4	56	Int_t AliHBTReaderESD::ReadNext()
7be13e68	57	{
bed069a4	58	//reads next event from fFile
	59
	60	//****** Tentative particle type "concentrations"
	61	static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
7be13e68	62
	63	Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
	64	Double_t w[kNSpecies];
	65	Double_t mom[3];//momentum
	66	Double_t pos[3];//position
7be13e68	67
bed069a4	68	if (AliHBTParticle::fgDebug)
	69	Info("ReadNext","Entered");
	70
7be13e68	71	TDatabasePDG* pdgdb = TDatabasePDG::Instance();
	72	if (pdgdb == 0x0)
	73	{
bed069a4	74	Error("Next","Can not get PDG Database Instance.");
7be13e68	75	return 1;
7be13e68	76	}
bed069a4	77
	78	if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
	79	if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
	80
	81	fParticlesEvent->Reset();
	82	fTracksEvent->Reset();
	83
7be13e68	84	do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
7be13e68	85	{
bed069a4	86	if (fFile == 0x0)
7be13e68	87	{
bed069a4	88	fFile = OpenFile(fCurrentDir);//rl is opened here
bed069a4	89	if (fFile == 0x0)
7be13e68	90	{
bed069a4	91	Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
	92	fCurrentDir++;
	93	continue;
	94	}
	95	fCurrentEvent = 0;
	96	}
	97
	98	//try to read
	99	TString esdname;
	100	esdname+=fCurrentEvent;
	101	AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
	102	if (esd == 0x0)
	103	{
	104	if (AliHBTParticle::fgDebug > 2 )
	105	{
	106	Info("ReadNext","Can not find ESD object named %s.",esdname.Data());
	107	}
	108	fCurrentDir++;
	109	delete fFile;//we have to assume there is no more ESD objects in the fFile
	110	delete fRunLoader;
	111	fFile = 0x0;
	112	continue;
	113	}
	114	AliStack* stack = 0x0;
	115	if (fReadParticles && fRunLoader)
	116	{
	117	fRunLoader->GetEvent(fCurrentEvent);
	118	stack = fRunLoader->Stack();
	119	}
	120	Info("ReadNext","Reading Event %d",fCurrentEvent);
7be13e68	121
bed069a4	122	Int_t ntr = esd->GetNumberOfTracks();
	123	Info("ReadNext","Found %d tracks.",ntr);
	124	for (Int_t i = 0;i<ntr; i++)
	125	{
	126	AliESDtrack *esdtrack = esd->GetTrack(i);
	127	if (esdtrack == 0x0)
	128	{
	129	Error("Next","Can not get track %d", i);
	130	continue;
	131	}
	132	if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE)
	133	{
	134	if (AliHBTParticle::fgDebug > 2)
	135	Info("ReadNext","Particle skipped: PID BIT is not set.");
	136	continue;
	137	}
7be13e68	138
bed069a4	139	esdtrack->GetESDpid(pidtable);
	140	esdtrack->GetPxPyPz(mom);
	141	esdtrack->GetXYZ(pos);
	142
3c69895b	143	Double_t extx;
	144	Double_t extp[5];
	145	esdtrack->GetExternalParameters(extx,extp);
	146
	147	Int_t charge = (extp[4] > 0)?-1:1;//if curvature=-charg/Pt is positive charge is negative
	148
bed069a4	149	//Particle from kinematics
	150	AliHBTParticle* particle = 0;
	151	Bool_t keeppart = kFALSE;
	152	if ( fReadParticles && stack )
	153	{
	154	if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track
	155	TParticle *p = stack->Particle(esdtrack->GetLabel());
	156	if (p==0x0)
	157	{
	158	Error("ReadNext","Can not find track with such label.");
	159	continue;
	160	}
3c69895b	161	// if(p->GetPdgCode()<0) charge = -1;
bed069a4	162	particle = new AliHBTParticle(*p,i);
3c69895b	163
bed069a4	164	}
	165
	166	//Here we apply Bayes' formula
	167	Double_t rc=0.;
	168	for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
	169	if (rc==0.0)
	170	{
	171	if (AliHBTParticle::fgDebug > 2)
	172	Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
	173	continue;
	174	}
7be13e68	175
bed069a4	176	for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
7be13e68	177
bed069a4	178	if (AliHBTParticle::fgDebug > 4)
	179	{
	180	Info("ReadNext","###########################################################################");
	181	Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
	182	Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
	183	TString msg("Pid list got from track:");
	184	for (Int_t s = 0;s<kNSpecies;s++)
	185	{
	186	msg+="\n ";
	187	msg+=s;
	188	msg+="(";
3c69895b	189	msg+=charge*GetSpeciesPdgCode((ESpecies)s);
bed069a4	190	msg+="): ";
	191	msg+=w[s];
	192	msg+=" (";
	193	msg+=pidtable[s];
	194	msg+=")";
	195	}
	196	Info("ReadNext","%s",msg.Data());
	197	}//if (AliHBTParticle::fgDebug>4)
7be13e68	198
bed069a4	199	for (Int_t s = 0; s<kNSpecies; s++)
bed069a4	200	{
dd82cadc	201	Float_t pp = w[s];
dd82cadc	202	if (pp == 0.0) continue;
7be13e68	203
3c69895b	204	Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
bed069a4	205	if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type
7be13e68	206
bed069a4	207	Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
bed069a4	208	Double_t tEtot = TMath::Sqrt( mom[0]mom[0] + mom[1]mom[1] + mom[2]mom[2] + massmass);//total energy of the track
7be13e68	209
bed069a4	210	AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i,
	211	mom[0], mom[1], mom[2], tEtot,
	212	pos[0], pos[1], pos[2], 0.);
	213	//copy probabilitis of other species (if not zero)
	214	for (Int_t k = 0; k<kNSpecies; k++)
	215	{
	216	if (k == s) continue;
	217	if (w[k] == 0.0) continue;
3c69895b	218	track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
bed069a4	219	}
	220	if(Pass(track))//check if meets all criteria of any of our cuts
	221	//if it does not delete it and take next good track
	222	{
	223	delete track;
	224	continue;
	225	}
	226
	227	fTracksEvent->AddParticle(track);
	228	if (particle) fParticlesEvent->AddParticle(particle);
	229	keeppart = kTRUE;
	230
	231	if (AliHBTParticle::fgDebug > 4 )
	232	{
	233	Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
	234	track->Print();
	235	}
	236	}//for (Int_t s = 0; s<kNSpecies; s++)
	237
	238	if (keeppart == kFALSE) delete particle;//particle was not stored in event
	239
	240	}//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
	241
	242	Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
dd82cadc	243	fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(),
bed069a4	244	fNEventsRead,fCurrentEvent,fCurrentDir);
7be13e68	245
bed069a4	246	fCurrentEvent++;
	247	fNEventsRead++;
	248	delete esd;
	249	return 0;//success -> read one event
	250	}while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
	251
	252	return 1; //no more directories to read
	253	}
	254	/**********************************************************/
	255	void AliHBTReaderESD::Rewind()
	256	{
	257	delete fFile;
	258	fFile = 0;
	259	delete fRunLoader;
	260	fCurrentDir = 0;
	261	fNEventsRead = 0;
	262	fCurrentEvent++;
7be13e68	263
bed069a4	264	}
7be13e68	265	/**********************************************************/
bed069a4	266
7be13e68	267	TFile* AliHBTReaderESD::OpenFile(Int_t n)
7be13e68	268	{
bed069a4	269	//opens fFile with kine tree
7be13e68	270
	271	const TString& dirname = GetDirName(n);
	272	if (dirname == "")
	273	{
	274	Error("OpenFiles","Can not get directory name");
	275	return 0x0;
	276	}
	277	TString filename = dirname +"/"+ fESDFileName;
	278	TFile *ret = TFile::Open(filename.Data());
	279
	280	if ( ret == 0x0)
	281	{
bed069a4	282	Error("OpenFiles","Can't open fFile %s",filename.Data());
7be13e68	283	return 0x0;
	284	}
	285	if (!ret->IsOpen())
	286	{
bed069a4	287	Error("OpenFiles","Can't open fFile %s",filename.Data());
7be13e68	288	return 0x0;
	289	}
	290
bed069a4	291	if (fReadParticles)
	292	{
	293	fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
	294	if (fRunLoader == 0x0)
	295	{
	296	Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
	297	delete ret;
	298	return 0x0;
	299	}
	300	fRunLoader->LoadHeader();
	301	if (fRunLoader->LoadKinematics())
	302	{
	303	Error("Next","Error occured while loading kinematics.");
	304	delete fRunLoader;
	305	delete ret;
	306	return 0x0;
	307	}
	308	}
7be13e68	309	return ret;
	310	}
	311	/**********************************************************/
	312
	313	Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
	314	{
	315	//returns pdg code from the PID index
	316	//ask jura about charge
	317	switch (spec)
	318	{
	319	case kESDElectron:
	320	return kPositron;
	321	break;
	322	case kESDMuon:
	323	return kMuonPlus;
	324	break;
	325	case kESDPion:
	326	return kPiPlus;
	327	break;
	328	case kESDKaon:
	329	return kKPlus;
	330	break;
	331	case kESDProton:
	332	return kProton;
	333	break;
	334	default:
	335	::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
	336	break;
	337	}
	338	return 0;
	339	}